MOLOCH computer code for molecular-dynamics simulation of processes in condensed matter

. Theoretical and experimental investigation into properties of condensed matter is one of the mainstreams in RFNC-VNIITF scientiﬁc activity. The method of molecular dynamics (MD) is an innovative method of theoretical materials science. Modern supercomputers allow the direct simulation of collective e ﬀ ects in multibillion atom sample, making it possible to model physical processes on the atomistic level, including material response to dynamic load, radiation damage, inﬂuence of defects and alloying additions upon material mechanical properties, or aging of actinides. During past ten years, the computer code MOLOCH has been developed at RFNC-VNIITF. It is a parallel code suitable for massive parallel computing. Modern programming techniques were used to make the code almost 100% e ﬃ cient. Practically all instruments required for modelling were implemented in the code: a potential builder for di ﬀ erent materials, simulation of physical processes in arbitrary 3D geometry, and calculated data processing. A set of tests was developed to analyse algorithms e ﬃ ciency. It can be used to compare codes with di ﬀ erent MD implementation between each other.


Short-range molecular dynamics challenges
There are several challenges to deal with while implementing short-range MD, including neighbours search within a predefined radius and decomposition algorithms to parallelize between computer nodes. The better they are implemented, the more flexible the code is.

Neighbours search within predefined radius
For the fast calculation of forces on atoms, a cutoff radius is used in MD, i.e., an atom is assumed to interact only with atoms within the cutoff radius.
One of the approaches to MD computing is the creation of neighbours lists for each atom. For this end, atoms within the cutoff radius need to be found as fast as possible.
The Link Cell Method (LCM), which is easy-to-implement and quick, is commonly used for this purpose. MOLOCH implements an original algorithm called the Ordered Space Method (OSM).
The LCM divides space ( Fig. 1(b)) into cubes with a constant rib length that is equal to or smaller than the cutoff radius, and links all atoms into these cubes ( Fig. 1(a)). Atom's neighbours are located in 27 nearest cubes. The LCM is simple in coding and quick in use, but its main disadvantage is memory failure when the sample is expanding (Fig. 1(d)).
This LCM disadvantage originates from its incontestable advantage -the constant rib length. In our method, the length is not constant, but this makes atom neighbours finding much more difficult. a e-mail: f.a.sapozhnikov@vniitf.ru This is an Open Access article distributed under the terms of the Creative Commons Attribution-Noncommercial License 3.0, which permits unrestricted use, distribution, and reproduction in any noncommercial medium, provided the original work is properly cited.  The OSM divides space into boxes with approximately equal atomic concentrations ( Fig. 1(c)) to get rid of the dependence on sample size. It is slower than the LCM by 10% in LJ like potentials (a 10 Å cutoff radius and 140 neighbours) and by 20% in MEAM like potentials (a 5 Å cutoff and 14 neighbours), but works without memory failure.

Parallelization method
Parallel MD techniques use force, atom and spatial decompositions. As shown in [1], the spatial decomposition is most effective when the numbers of atoms and processors are large. MOLOCH uses a hierarchical spatial decomposition (Fig. 2) which allows processor responsibility zones to be excellently determined with respect to local atomic concentrations. It also allows effective dynamic balancing with processor responsibility zones being locally changed simultaneously with the exchange of atoms either at each time step, or when imbalance threshold is exceeded.
The algorithm of zoning includes the following steps.
1. Cut a sample along the X axis into slices with approximately equal number of atoms ( Fig. 2(b)); 2. Cut each slice, independently of the others, along the Y axis into bars with approximately equal number of atoms ( Fig. 2(c)); 3. Cut each bar, independently of the others, along the Z axis into zones with approximately equal number of atoms ( Fig. 2(d)); 4. Each processor is assigned a zone of space with atoms (or without them if balancing simply does not have time to change the zones at a high rate of atomic flow).

MOLOCH code
An MD calculation implies three phases: sample generation, MD simulation and results analysis. Reliability of the entire calculation is questionable if at least one of the phases gives non-physical results. MOLOCH implements effective algorithms that ensure reliable results.

00017-p.2
New Models and Hydrocodes for Shock Wave Processes in Condensed Matter

Sample generation
For getting physically true samples, MOLOCH implements a generator of single-and poly-crystals by the Voronoy cell method ( Fig. 3(a)). It is possible to add dopes and defects. There is a special library of pre-calculated defects, for example, a helium bubble in plutonium lattice ( Fig. 3(b)) or a loop dislocation (Fig. 3(c)). A superposition of geometrical constraints can be used: half plane, sphere, cylinder, box ( Fig. 3(d)).
We developed and implemented in the code the following potentials: GEAM (Generalized Embedded Atom Model) for materials with complicated electronic structures and Modified REBO ABCD (a long-range modification of REBO) for a model explosive containing atoms of 4 types (parameterization for TATB) [2].
Parameters for the existing and new potentials are adjusted with a parallel genetic algorithm which minimizes error in the description of a given set of experimental or calculated data by a particular potential.
MOLOCH is written in C++ which allows new potentials to be easily implemented at minimal changes in the code.

Simulation
The NVE, NVT, and NPT ensembles can be modeled with the simultaneous use of several potentials and a variable time step. Complicated geometry constraints are implemented (e.g., a moving collapsing spherical piston).

Results analysis
MOLOCH creates 3D meshes for macroscopic characteristics (temperature, pressure and density), detects bonds and calculates chemical compositions (e.g., during carbon nanocluster evaporation modeling).

Automatic structure recognition
We developed a method called Adaptive Template Analysis (ATA) to recognize crystal structures and defects [3]. ATA recognizes BCC, FCC, HCP, Diamond and Graphen lattices (other ones can be easily added), and reconstructs vacancies, identifies interstitials and split interstitials, dislocations and stacking faults. ATA doesn't have any adjustable parameters, has low sensitivity to temperature disordering and local density fluctuations.

Test set for MD codes algorithms
There are many tests that can be used to verify if MD results are physically correct but we also need tests that could help measure the efficiency of each particular algorithm and the entire MD implementation. The test setups provided below are devised to test such important algorithms as dynamic balancing, multi-potential modelling and neighbours list building.

00017-p.4
New Models and Hydrocodes for Shock Wave Processes in Condensed Matter a) b) c) MOLOCH uses sum of neighbours list building time and energy and force calculation time to balance processors loading. Let's calculate imbalance (D) using following formula:

"Death Star" test
The test aim is to measure the efficiency of dynamic balancing (if implemented) which must be quick (to balance on every MD step if high imbalance) and independent of the axial direction of atomic flows. The problem implements a fully 3D deformation of a thin spherical layer (Fig. 8).
Test sample: A zero-centred spherical layer is cut from single crystal FCC copper (EAM [4]) with outer and inner radii R big and R small = 0.95R big , respectively ( Fig. 7(a)). Velocities of all atoms are pointed to the centre (Fig. 7(b)) and defined as: v = V max e − arccos(cos θ cos θ 0 +sin θ sin θ 0 cos(φ−φ 0 )) , here V max = 3000m/s is maximal velocity. (1) The spherical coordinates of a specified direction of loading are φ 0 = π/6 and θ 0 = π/3. The spherical coordinates of an atom are φ = atan2(y/R big , x/R big ) and θ = acos(z/R big ) where x, y, z are its Cartesian coordinates.
Simulation: Free boundary conditions; timestep 1.5 fs, 30'000 timesteps. Outcome: The ratio μ of test sample calculation time to the control sample calculation time is to be close to 1. If it is much greater, then the balancing efficiency is low.

"Cocktail B-52" test
The test shows the efficiency of MD simulation with several potentials with strongly different cutoff radii, specifically MEAM (3.5 Å) and LJ (10 Å). Physical analogue: Modelling of interface instability. Problem: Time of MEAM triple interactions is extremely sensitive to the cutoff value, so it may slow down several times on looking for interaction within 'useless' for this potential large cutoff.
Control samples: Independent calculations of HCP Be and FCC Ar cubes ( Fig. 9(b) and 9(c)). Simulation: Free boundary conditions, timestep 1.5fs, 1000 timesteps. Outcome: The ratio μ of test sample calculation time to the sum of control samples calculation times is to be close to 1. If it is much greater, then the efficiency of multi-potential modelling is low.

"Supernova" test
The test demonstrates feasibility of MD simulation for expanding samples.