Parallel Molecular Dynamics Simulation on Elastic Properties of Solid Argon
Futoshi Shimizu, Hajime Kimizuka, Hideo Kaburaki
Center for Promotion of Computational Science and Engineering
Japan Atomic Energy Research InstituteNakameguro 2-2-54, Japan 153-0061
shimizu@koma.jaeri.go.jp kimizuka@koma.jaeri.go.jp kaburaki@sugar.tokai.jaeri.go.jp
Ju Li, Sidney Yip
Department of Nuclear Engineering
Massachusetts Institute of Technology
Cambridge, Massachusetts 02139
liju99@mit.edu syip@mit.edu

Abstract
Parallel Molecular Dynamics Stencil has been developed to execute eﬀectively large-scaleparallel molecular dynamics simulations. The Stencil is adapted to varieties of molecular dynamics simulations without special attention to parallelization techniques. As an example oflarge-scale simulation using this Stencil, the adiabatic elastic constants of solid argon in crystalline and amorphous states, have been evaluated over the temperature range of 15–75K. An
MD model of solid argon, consisting of 500–1,000,188 atoms interacting via the Lennard-Jones
potential, has been used throughout the simulations. It is found that the cutoﬀ-length is criticalin evaluating elastic constants and that the calculated values for the temperature dependenceof crystalline argon agree well with the measured ones.

1 Introduction
Molecular dynamics (MD) method is widely employed to predict macroscopic properties
from atomistic motions in diverse applications such as physics, chemistry, and material science.
In particular, equilibrium molecular dynamics is used to predict physical properties of matter,such as macroscopic thermodynamics quantities, atomistic structural properties, and dynamicalquantities of transport coeﬃcients [1].
In today’s computing environment, the distributed computing is the way to extend the system to a large scale. Parallelization is the basis towards this future distributed computing.
Many parallel MD simulation programs, such as AL CMD [2], IMD [3], Moldy [4], NAMD [5],and SPaSM [6], have been developed and published. Almost all of these programs introduce
one speciﬁc parallel algorithm while there are many kinds of parallel algorithms for MD. In thispaper, we take up two parallel algorithms, that is, particle decomposition and spatial decomposition methods, and make a comparison between them by our newly developed MD program,
Parallel Molecular Dynamics Stencil [7]. We also report an MD simulation of solid argon as anapplication of the Stencil.
This paper is published in SNA 2000, Proceedings of the Fourth International Conference on Supercomputing
in Nuclear Applications, September 4-7, 2000, Tokyo, Japan.

2 Parallel MD simulation
In the method of equilibrium molecular dynamics, the size eﬀect would become more important when we derive the statistically averaged results, in particular, in the lower temperature
region. The introduction of the massively parallel supercomputers leads us to the possibilityof attacking larger systems, but also to the involvement in the complexities of parallel calculations along with choosing varieties of potentials and statistical ensemble methods for the MD
method. In this section, we describe an empirical Lennard-Jones (LJ) potential, cutoﬀ schemes,
and parallel algorithms.
2.1 Lennard-Jones potential
As a model system for the method of MD simulation, a system interacting with the LennardJones pair potential has long been used to elucidate the fundamental properties of many-bodyinteracting particles. The standard Lennard-Jones 12-6 potential for argon below is suitable
for measuring all the computational eﬃciency of the parallel MD simulation.
σ
σφ(r) = 4ε[( )12 − ( )6 ],
r
r

(1)

where parameters σ(= 0.34[nm]) and ε(= 120[K] × kB ) are the distance to the zero in φ(r) andthe energy at the minimum in φ(r), respectively. Also, the prediction of elastic constants has
been found to be a good benchmark to test the accuracy of the empirical potential and thenumerical eﬃciency of the MD method [8, 9].
2.2 Cutoﬀ schemes
In a system of N particles interacting via pair potential, there are N(N − 1)/2 calculations in
total. Thus, the cost of force calculation increases with the square of N, and almost all execution
time of MD simulation is taken up in the calculation of force as the system size becomes larger.
The cutoﬀ distance rc is introduced in the LJ system to save the computational cost. The
shifted-force potential is used throughout this paper where pair interactions beyond rc areneglected and the discontinuity involved in calculating the force is removed at the cutoﬀ point.


dφ
dφ 

dΦ
r̂
(r < rc )− r̂ +r̂ =f (r) =
(2)drdr rcdr0(r ≥ rc ), φ(r) − φ(r ) − (r − r ) dφ cc
Φ(r) =
dr rc0(r < rc )(3)(r ≥ rc ).

The pair list is introduced for the eﬀective calculation of the shifted-force potential system.
The one method is to use the typical Verlet neighbor list, where each (i-th) particle in a system
has a list of neighboring (j-th) particles whose positions are within a distance of rlist (Figure1(a)). The distance rlist is taken as the sum of rc and some buﬀer margin. The updating intervalof the neighbor list is determined by this buﬀer margin. With this method, the computational
cost is proportional to the number of total particles N and the numbers of particles in the list
which is independent of N.

r celliijr list(a) Verlet neighbor list methodj(b) Cell partitioning method

Figure 1: Cutoﬀ schemes for short range interactions.

The other method is cell partitioning. The simulation space is divided into rectangular cells
whose side length rcell is taken larger than rc (Figure 1(b)). In this method, the pair list for an
i-th particle in a particular cell can be constructed by searching particles in the cell containing
the i-th particle and the contacting neighboring cells. The computational cost for this method
becomes also proportional to the number of total particles N and the numbers of particles incells which is independent of N.
2.3 Parallel algorithms
Here we consider the parallelization of the MD simulation program, in which the data structures such as positions and velocities of particles and the computing procedures such as calculations of forces and updates of particle positions are partitioned and allotted to each processor.
One of the key factor to be taken into account for the parallel computation is the communication between processors. The communication time, consisting of the data transfer time andthe start-up time, always tends to reduce the parallel eﬃciency. Also, in the case that communications among processors are synchronous, load imbalance, in which one processor with thejob done has to wait other processors to ﬁnish their jobs, deteriorates the parallel eﬃciency.
Since the method of partitioning the data and computational procedures is closely related, agood parallel MD program should be made in such a way to lessen the communication data
and frequencies and to realize a good load balance among processors. Next we describe thetypical partitioning methods, particle and spatial decomposition method [10].
In the particle decomposition method, particles are globally numbered, partitioned, andallotted to processors, as shown in Fig.2(a). Basically, each processor takes care of the allotted
particles from the beginning to the end of the calculation. Here, for simplicity of the program,the replicated data method is used where all processors store and share positions of particles
for the whole system. It is predicted to be disadvantageous that scalability in this systemdeteriorates as the number of processors increases, since the synchronous global communicationof broadcasting is required. On the other hand, load balance among processors can be easilyadjusted to be uniform by allotting the same number of particles to each processor.
In the spatial decomposition method, communication cost can be reduced using the localityof particles, that is, by allotting particles belonging to the partitioned space to each processor(see Fig.2(b)). In this method, communication between neighboring processors is required to
calculate the force, because particles near the boundary region interact with those in neigh-simulation cell
10

0
21
13

20

26
4

33

29

7

23
34

9
10

#0

28

19
9

8
39

0

1

3

19
20

#1

29
30

#2

8
1

2

4
7

6

9 10
1

0

2
3

3

#3

6

4

7

(a) Particle decomposition method

7

5
8

9

processors
0
8
0

#0

11
0

#1

7
0

#2

10

#3

8

11

5
4

39

3
5

5
7
0

2

1

0
4

6

35
27

simulation cell
2

24

5
17

38

2

14

16

37

11

31 32

22
12

3

36

18

1

30

25
15

6

processors
0610(b) Spatial decomposition method

Figure 2: Two decomposition methods for parallel MD; each forty particles are allotted
to 4 processors.

boring regions. Also, sending and receiving particles are necessary among processors whenparticles move across the boundary. Scalability for the program execution time as a functionof the number of processors maintains well since the communication is restricted to local areas.
Load balance among processors is generally good except that the distribution of particles is
extremely imbalance.
2.4 Development of Parallel Molecular Dynamics Stencil
It has not been straightforward to implement the parallel algorithms methods, into the MD
program, writing necessary procedures for data communication explicitly. Combination withthe cutoﬀ schemes for eﬃcient computation of short range particle interactions makes the programming furthermore diﬃcult and cumbersome. We have been developing Parallel Molecular
Dynamics Stencil [7], which consists of parallel model programs for executing MD simulations
by incorporating various statistical ensembles methods to realize constant temperature andconstant pressure conditions, and the integration algorithms such as velocity Verlet, Beeman,and ﬁve-value Gear methods. The Stencil is designed in such a way to separate and hide partsof the programs for cutoﬀ schemes and parallel algorithms, so that users can concentrate on
the numerical simulation itself. Although the whole program is written in C using MPI [11],users do not need to know how to use each MPI function.
In the present version of the Stencil, the pair list for cutoﬀ schemes is taken as the interface
to calculate force. This part of the program can be revised according to other physical models.
With the Stencil’s framework, parallel programming can be done in the same way as sequential
programming for a single processor. This stage of programming is designed to be independent ofother programming stages of cutoﬀ schemes, integration algorithms, and statistical ensembles.
Any of these routines can be easily linked to construct a parallelized executable program.
2.5 Performance of Parallel MD Stencil
The rare-gas elements are good examples for the calculation of MD properties since they
have comparatively simple structures, and weak (van der Waals) forces can be represented bywell-deﬁned short-range interactions. We have evaluated performance of Parallel MD Stencil
on Intel Paragon XP/S75 with argon atoms interacting via the shifted-force LJ potential ofElapsed time per update: t [sec]

(a)
256
64

100

SD/CP
PD/NL
2,048,000
256,000
32,000
4,000
500

SD/CP
PD/NL
2,048,000
256,000
32,000
4,000
500

(b)
10
Updates: S [/sec]

1024

16
4

1

0.1

1
0.01
0.25
0.001

0.0625
1

4
16
64
256
Number of processors: p

1024

1

10

100

1000

10000

100000

Number of atoms per processor: nFigure 3: Performance of Parallel MD Stencil; (a) time versus number of processorsand (b) speed versus number of atoms per processor are shown by two methods, spatial
decomposition with cell partitioning method (solid line) and particle decomposition
with Verlet neighbor list method (dashed line).

Eq.3. Here, Nosé-Hoover [12] thermostat and Parrinello-Rahman [13] constant-stress method
are incorporated to realize the NT P ensemble. We chose the cutoﬀ length to be 2.5σ(0.85[nm])and measured the performance over many cases of simulations using atoms ranging from 500
to 2,048,000.
Figure 3(a) shows elapsed time per update of parallel execution by two methods; spatial
decomposition with cell partitioning method (solid line) and particle decomposition with Verlet
neighbor list method (dashed line). Because of the limitation of memories of each processor,
the maximum number of atoms that can be treated by the particle decomposition method is
32,000 in the present computing environment. As the total number of atoms decreases, there
is an upper limit to the number of processors that can be used in the spatial decomposition
method. This limitation comes from the condition of interaction range, in which the side lengthof processor cell has to be taken larger than rc . On the other hand, in the particle decomposition
method, it is noted that, in principle, even one atom can be assigned to each processor without
regard to its eﬃciency. Figure 3(b) represents the speed S as a function of the number of atoms
per processor n. In the spatial decomposition method, scalability is conﬁrmed by the fact thatdata for diﬀerent total number atoms collapse onto a single diagonal line in the region ofn > 1, 000. This is because the region for search to make the pair list based on cutoﬀ algorithm
is restricted to local areas. By ignoring inter-processor communications, the execution time per
update is estimated to be tSD = C(An + Bn), where A, B, C are the constants. Here, An and
Bn stand for costs for force calculation and pair list making, respectively. The logarithm of
the number of updates per second SSD = 1/tSD is written aslog SSD = − log n + D,

(4)

where D = − log C − log(A + B).
In the particle decomposition method, on the other hand, there are three diﬀerent lines forcases of 32,000 atoms, 4,000 atoms, and 500 atoms. Since all atoms are objects of search for apair list construction, time per update is estimated to be tP D = C (An+qB n), where A, B , C 
are the constants. A parameter q = N/n is nearly equal to the number of processors p. Anand qB n stand for costs for force calculation and pair list making, respectively. In this case,data for diﬀerent total number atoms do not conform to a single line, because the logarithm of
the number of updates per second SP D = 1/tP D depends on q:
log SP D = − log n + D (q).

(5)

Since D (q) = − log C  −log(A +qB ) is a function of q which depends on N, there appear three
separate lines for each N. In addition, owing to the collective communications which involve
all processors, the deviation from the diagonal line is observed for each N. In the case of 500
atoms, it is noted that the updates for the particle decomposition method are larger than thosefor the spatial decomposition method.
With all the performance described above, it becomes straightforward to extract the best
performance from your computing environment when the parameters for simulation, such as
number of atoms, cut-oﬀ length, etc., are ﬁxed. The Parallel MD Stencil is designed in such a
way to adapt easily to a wide range of parameter space.

3 Elastic Properties of Solid Argon
The spherically symmetric rare-gases atoms crystallize with face-centered cubic (fcc) structures, and their solids have been the subject of many theoretical investigations. Also, a quantitative agreement has been obtained for several of the experimentally determined properties. In
this study, we focused on the elastic properties of crystalline and amorphous argon solid, andinvestigated these system-size eﬀect to show the eﬀectiveness of parallel computing.
3.1 Elastic property of crystalline argon
We have performed MD simulations of the crystalline argon using several techniques employed in the Parallel MD Stencil. We started from the model system in the fcc structure(500–1,000,188 atoms in a cubic cell with periodic boundary conditions, interacting via theshifted-force LJ potential of Eq.3 with the cutoﬀ-length rc =2.5σ or 4.1σ), and heated the system
continuously to the desired temperature. The structures of the reference states at the desired
temperature and stress were determined in the NT P ensemble using the Nosé-Hoover [12] thermostat and Parrinello-Rahman [13] constant-stress methods. The adiabatic elastic constants Cij
at each temperature were calculated within the constant-volume and constant-energy (NV E)
ensemble. The ﬂuctuation formula for the internal stress tensor [14] was used to obtain acomplete set of elastic constants over the temperature range of 15–75K.
Figure 4 shows the temperature dependence of the Cij ’s of crystalline argon. The calculated
C11 has a smaller temperature dependence than that of the experiment whereas the calculated
values of C12 and C44 tend to have a larger temperature dependence [15]. However, it is
conﬁrmed that the overall temperature dependence is well reproduced. In particular, the model
agrees with experiment in predicting that C44 decreases larger with temperature than C12 . It
is generally known that the Cauchy relation C12 = C44 holds in the cubic symmetry if theatoms interact with central forces such as the LJ potential. (Also, the Poisson ratio will be
estimated to be 1/4.) In Fig. 4, the calculated C12 and C44 values are distinct from each other,

5.0
C11

4.5
Adiabatic elastic constants (GPa)

N=4,000
4.0
3.5

N=1,000,000

3.0
2.5

C11
C12

2.0
1.5

C12

C44

C44

1.0
0.5
0.0
0

10

20

30

40

50

60

70

80

90

Temperature (K)Figure 4: Temperature evolution of the adiabatic elastic stiﬀness constants of crystalline
argon. The calculated elastic constants for rc = 4.1σ and 2.5σ are shown as solid andopen plots, respectively; C11 (circle), C12 (square), and C44 (triangle). The experimental
values for crystalline argon [15] are shown by curves.

10.0

Mean square displacement (Å2)

9.0
N=500

8.0

N=1372

7.0
6.0

N=500

N=4000

5.0

45K

4.0
30K

3.0
2.0
1.0

15K

0.0
0

5000

10000 15000

20000

25000

30000 35000

40000

Simulation time step

Figure 5: Time evolution of the mean-square displacements of argon atoms.

and this result agrees with experiment. It indicates that the present results reproduce the ﬁnite
temperature eﬀects in an exact manner. However, the Cauchy relation is shown to hold in theextremely low-temperature region, since the elastic constants are expressed by the ﬂuctuationof a certain ensemble and this ﬂuctuation approaches zero in the limit of T =0K.
Comparing the two cutoﬀ-lengths, the Cij values for rc =2.5σ are signiﬁcantly smaller thanthose for rc =4.1σ, and the temperature dependence is underestimated in this temperature
range. It depends on the number of interactions taken into account to evaluate the elastic
constants. The cutoﬀ-length eﬀect is important for the elastic properties and the relatively
long-range interactions are recommended in this method.
3.2 Elastic property of amorphous argon
An MD model of amorphous argon has been obtained from a model of a supercooled liquid (T =100K) by relaxing in the NT P ensemble at constant zero pressure using the abovementioned methods.
Figure 5 shows the time evolution of the mean-square displacements (MSD) of argon atoms.
Overall, the MSD values increase linearly with time, and the slope of the MSD curve, which isequivalent to the diﬀusion coeﬃcient, becomes steep as the temperature increases. It suggeststhat the argon atoms be loosely jointed, and drift gradually keeping the amorphous states.
Rapid changes in the MSD curves indicate the transformation to the more stable phase due tocrystallization. Moreover, it reveals that the duration (lifetime) of amorphous state depends
on the system size.
Figure 6 shows the temperature dependence of the Cij ’s of amorphous argon. It is found that
the Cij values of the model of amorphous argon satisfy adequately the condition for an isotropic
material: C44 = 12 (C11 − C12 ). It is noted that C44 is signiﬁcantly small compared with the5.0

Adiabatic elastic constants (GPa)

4.5
4.0
3.5
C11

3.0
2.5

C11

C12

2.0

C12

1.5

C44

1.0
C44

0.5
0.0
0

10

20

30

40

50

60

70

80

90

Temperature (K)Figure 6: Temperature evolution of the adiabatic elastic stiﬀness constants of solidargon. The calculated elastic constants of crystalline and amorphous argon are shown
as solid and open plots, respectively; C11 (circle), C12 (square), and C44 (triangle).

other components, and a Poisson ratio exhibits the value close to the isotropic upper limit of0.5, as found in rubbery materials. The Cij values tend to decrease with temperature, and C44
particularly disappears (approaches zero) at around 60K. It indicates that the obtained solid
phase becomes mechanically unstable and transforms to a liquid phase at this temperature. It
contrasts with the fact that argon solids remain in the fcc structure until about 84K (melting
temperature).

4 Conclusions
We have implemented atom and spatial decomposition methods into Parallel Molecular Dynamics Stencil. By performing the evaluation analysis of the Stencil using argon atoms, a guiding principle in parallel MD simulation has been obtained; the spatial decomposition method
achieves a good parallel eﬃciency for system with a large number of atoms, while the particle
decomposition method is eﬃcient for systems with a small number of atoms.
The Parallel MD Stencil was found to be useful in simulations for searching parameters in awide range, in particular, changing the system size. It is found that the cutoﬀ-length of at leastrc = 4.1σ is needed in evaluating the elastic properties of argon solid. The simulation reproduces
well the absolute values and the temperature dependence of adiabatic elastic constants ofcrystalline argon, and, with the same method, predicts elastic constants of amorphous argon.

References[1] J. M. Haile, Molecular Dynamics Simulation: Elementary Methods, A Wiley-Interscience
Publication (1992).
[2] Ames Lab Classical Molecular Dynamics Source Code,
http://cmp.ameslab.gov/cmp/CMP Theory/cmd/alcmd source.html.
[3] The ITAP Molecular Dynamics Package Homepage,
http://www.itap.physik.uni-stuttgart.de/˜imd/.
[4] A GENERAL-PURPOSE MD CODE, http://www.earth.ox.ac.uk/˜keith/moldy.html.
[5] NAMD – Scalable Molecular Dynamics, http://www.ks.uiuc.edu/Research/namd/.
[6] SPaSM’s home, http://bifrost.lanl.gov/MD/MD.html.
[7] F. Shimizu, H. Kimizuka and H. Kaburaki, Proc. of the Conf. on Comp. Eng. and Sci. 5,361 (2000), (in Japanese).
[8] M. Tang and S. Yip, J. Appl. Phys. 76, 2719(1994).
[9] H. Kimizuka, H. Kaburaki, and Y .Kogure, Phys. Rev. Lett. 84, 5548 (2000).
[10] D. M. Beazley, et.al, Parallel Algorithm for Short-Range Molecular Dynamics, World Scientiﬁc’s Annual Reviews in Computational Physics, 3, 119 (1995).
[11] MPI - The Message Passing Interface Standard, http://www-unix.mcs.anl.gov/mpi/.
[12] S. Nosé, Mol. Phys. 52, 255 (1984); W. G. Hoover, Phys. Rev. A 31, 1695 (1985).
[13] M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).
[14] J. R. Ray and A. Rahman, J. Chem. Phys. 80, 4423 (1984); J. R. Ray, Comput. Phys.
Rep. 8, 109 (1988).
[15] G. L. Keeler and D. N. Batchelder, J. Phys. C: Solid St. Phys. 3, 510 (1970).