**Large-Scale Molecular Dynamics Simulations of Materials on Parallel
Computers**

Aiichiro Nakano,^a^ Timothy J. Campbell,^a,b^ Rajiv K. Kalia,^a^ Sanjay
Kodiyalam,^a^ Shuji Ogata,^a,c^ Fuyuki Shimojo,^a,d^ Priya Vashishta,^a^
and Phillip Walsh^a^

^a^Concurrent Computing Laboratory for Materials Simulations

Department of Computer Science, Department of Physics & Astronomy

Louisiana State University, Baton Rouge, LA 70803, USA

^b^Logicon Inc., Naval Oceanographic Office Major Shared Resource
Center, Stennis Space Center, MS 39529, USA

^c^Department of Applied Sciences, Yamaguchi University, Ube 755-8611,
Japan

^d^Faculty of Integrated Arts and Sciences, Hiroshima University,
Higashi-Hiroshima739-8521, Japan

**Abstract.** Scalable space-time multiresolution algorithms implemented
on massively parallel computers enable large-scale molecular dynamics
simulations involving up to a billion atoms, which are applied to the
study of nanosystems of great technological importance. These include
sintering, structure, and mechanical properties of nanostructured
ceramics and nanocomposites, structural transformation in semiconductor
nanocrystals, nanoindentation, and oxidation of metallic nanoparticles.

# 

# Introduction {#introduction .INTRODUCTION}

Advanced materials and devices with nanometer grain/feature sizes are
being developed to achieve higher strength and toughness in ceramic
materials and greater speeds in electronic devices. Below 100 nm,
however, continuum description of materials and devices must be
supplemented by atomistic descriptions \[1\]. Current state-of-the-art
atomistic simulations involve 1 million to 1 billion atoms. Scalable
multiresolution algorithms that enable these large-scale simulations are
described. We discuss the application of these algorithms to molecular
dynamics simulations of various nanosystems of great technological
importance.

# Scalable Atomistic Simulation Algorithms

We have developed a suite of scalable simulation, parallel-computing,
and data-management algorithms \[2,3\]. Billion-atom molecular dynamics
(MD) simulations have been demonstrated \[4\] on massively parallel
computers using: i) space-time multiresolution algorithms; ii)
wavelet-based adaptive load balancing in curvilinear computational
space; and iii) spacefilling-curve-based adaptive data compression for
scalable I/O.

We have also developed \[4\] scalable algorithms for quantum-mechanical
(QM) calculations based on the density functional theory (DFT) \[5\],
thereby demonstrating that 20,000-atom DFT-based QM calculations are
feasible on massively parallel computers.

## Space-time Multiresolution Molecular Dynamics

The MD approach obtains the phase-space trajectories (positions and
velocities of all atoms at all time) by numerically integrating coupled
ordinary differential equations. The dynamics is encoded in the
interatomic potential energy, *E*~MD~(**r***^N^*), which is a function
of the positions of all *N* atoms, **r***^N^* = {**r**~1~, **r**~2~,
\..., **r***~N~*}, in the system. In our many-body interatomic potential
scheme, *E*~MD~(**r***^N^*) is an analytic function of relative
positions of atomic pairs and triplets \[6\].

We have developed highly efficient, multiresolution algorithms to carry
out large-scale MD simulations on parallel computers. The most
compute-intensive problem in an MD simulation is the *O*(*N*^2^)
computation of the electrostatic energy for *N* charged atoms. We use
the Fast Multipole Method (FMM) by Greengard and Rokhlin to reduce the
complexity to *O*(*N*) by computing the electrostatic field recursively
on an octree \[7\]. To extend the simulated time scales, we combine the
FMM with a multiple time-scale (MTS) method that applies different
force-updating schedules for different force components \[8\].

We have implemented this multiresolution molecular-dynamics (MRMD)
algorithm on a number of parallel computers (1,280-processor IBM SP3,
1,088-processor Cray T3E, and 512-processor SGI Origin supercomputers as
well as on our 166-node PC and 64-node Digital Alpha clusters) using
spatial decomposition. The MRMD algorithm is highly scalable \[4\]: For
a 1.02-billion-atom silica system, one MD step takes only 26.4 seconds
on 1,024 Cray T3E processors. For this system, the parallel efficiency
is 0.97.

## Wavelet-based Adaptive Load Balancing

To overcome the load-imbalance problem associated with parallel
computation of irregular atomic distribution, we have developed an
adaptive load-balancing algorithm \[9\]. The
"computational-space-decomposition" scheme is based on a novel idea that
computational space shrinks where the workload density is high. The
algorithm introduces a curvilinear coordinate transformation to
partition workloads according to a uniform 3D-mesh topology in the
curved computational space. This leads to curved partition boundaries in
the physical Euclidean space. Subsequently simulated annealing is used
to minimize load-imbalance and communication costs as a functional of
the coordinate transformation. The multiresolution analysis based on
wavelets allows compact representation of partition boundaries, leading
to fast convergence of the minimization. The wavelet load balancer (WLB)
achieves nearly perfect speed-up with negligible computational overhead.

## Spacefilling-curve-based Adaptive Data Compression for Scalable I/O

For scalable input/output (I/O), we have designed a data compression
algorithm, which uses octree indexing and sorts atoms accordingly on the
resulting spacefilling curve \[10\]. By storing differences between
successive atomic coordinates, the I/O requirement for the same error
tolerance level reduces from *O*(*N*log*N*) to *O*(*N*). An adaptive,
variable-length encoding scheme is used to make the algorithm tolerant
to outliers and optimized dynamically. The spacefilling-curve data
compression (SCDC) algorithm achieves an order-of-magnitude improvement
in the I/O performance for actual MD data with user-controlled error
bounds.

## Linear-scaling Quantum-Mechanical Algorithms

An atom consists of a nucleus and surrounding electrons, and quantum
mechanics explicitly treats the electronic degrees-of-freedom. The DFT
reduces the exponentially complex quantum problem to a self-consistent
matrix eigenvalue problem, which can be solved with *O*(*N*~wf~^3^)
operations (*N*~wf~ is the number of independent wave functions, or
electronic bands). The DFT can be formulated as a minimization of the
energy, *E*~QM~(**r***^N^*, ψ^*N*wf^), with respect to electron wave
functions, ψ^*N*wf^(**r**) = {ψ~1~(**r**), ψ~2~(**r**), \...,
ψ~*N*wf~(**r**)}, subject to orthonormalization constraints. Efficient
parallel implementation of DFT is possible with real-space approaches
based on higher-order finite differencing \[11\] and multigrid
acceleration \[12\]. We include electron-ion interactions using
norm-conserving pseudopotentials \[13\] and the exchange-correlation
energy in a generalized gradient approximation \[14\]. For large systems
(*N*~wf~ \> 1,000), however, the *O*(*N*~wf~^3^) orthonormalization
becomes the bottleneck.

For scalable DFT calculations, linear-scaling algorithms are essential
\[15\]. We have implemented \[4\] an *O*(*N*~wf~) algorithm based on
unconstrained minimization of a modified energy functional and a
localized-basis approximation \[16\]. In the parallel linear-scaling DFT
algorithm, the computation time scales as *O*(*N*~wf~/*P*) on *P*
processors, whereas the communication time scales as
*O*((*N*~wf~*/P*)^2/3^). This is in contrast to the
*O*(*N*~wf~(*N*~wf~*/P*)^2/3^) communication in the conventional
parallel real-space DFT algorithm. Global communication for calculating
overlap integrals of wave functions (which scales as *N*~wf~^2^log*P* in
the conventional DFT algorithm) is unnecessary.

Our parallel linear-scaling DFT (LSDFT) algorithm is highly scalable:
For a 22,528-atom GaAs system on 1,024 Cray T3E processors, the parallel
efficiency is 0.96 \[4\]. Figure 1 shows the performance of our scalable
MD and QM algorithms on 1,024 T3E processors.

![](media/image1.png){width="2.4722222222222223in"
height="2.272222222222222in"}

**Figure 1.** Computation time of MD and QM algorithms as a function of
the number of atoms on 1,024 Cray T3E processors: (circles) classical MD
based on many-body interatomic potentials; (squares)
environment-dependant, variable-charge MD; (triangles) QM calculation
based on the tight-binding method; and, (diamonds) self-consistent
DFT-QM.

# Large-Scale Atomistic Simulations of Nanostructured

The scalable simulation algorithms described in the previous section
have been used to perform large-scale atomistic simulations of various
nanosystems. In this section, we summarize some of our simulation
results.

## Sintering, Structure, and Mechanical Properties of Nanostructured Ceramics and Nanocomposites

Advanced structural ceramics are highly desirable materials for
applications in extreme operating conditions. Light weight, elevated
melting temperatures, high strengths, and wear- and corrosion-resistance
make them very attractive for high-temperature and high-stress
applications. The only serious drawback of ceramics is that they are
brittle at low to moderately high temperatures. In recent years, a great
deal of progress has been made in the synthesis of ceramics that are
much more ductile than conventional coarse-grained materials \[17\].
These so called nanostructured materials are fabricated by *in-situ*
consolidation of nanometer size clusters.

We have performed multimillion-atom MD simulations to investigate
sintering, structure, and mechanical behavior of nanostructured
Si~3~N~4~ \[18-20\], SiO~2~ \[21\], and SiC \[22\].

Nanostructured Si~3~N~4~ is generated by consolidating a random cluster
configuration \[18-20\]. Pair distribution functions and bond angle
distributions reveal that interfacial regions in the consolidated
nanophase Si~3~N~4~ are amorphous. Systems sintered at low pressures
(1GPa) have percolating pores whose surface roughness exponents (0.46
and 0.86) are in excellent agreement with experiments. We also find that
the dependence of elastic moduli on porosity and grain size in
nanostructured Si~3~N~4~ can be understood in terms of a three-phase
model for heterogeneous materials.

Nanocomposite materials often exhibit superior properties compared with
conventional materials \[23\]. In recent years, large-scale MD
simulations have played a central role in the investigation of dynamic
fracture in such complex materials. For example, we are performing
1.5-billion-atom MD simulations to investigate fracture in a Si~3~N~4~
matrix reinforced with SiC fibers of linear dimension 0.3 μm. In order
to simulate the effect of a glassy phase that lubricates the
fiber-matrix interfaces, SiC fibers are coated with amorphous silica
layer (Fig. 2). We are investigating the effect of interphase structure
and residual stresses on the fracture toughness. Immersive and
interactive visualization reveals a rich diversity of atomistic
processes including fiber rupture, frictional pullout, and emission of
molecular fragments, which must all be taken into account in the design
of tough composites.

![](media/image2.png){width="2.808333333333333in" height="1.4in"}

**Figure 2.** Atomistic model of fractured Si~3~N~4~ matrix reinforced
with SiC fibers coated with amorphous silica. Small spheres represent
silicon atoms, and large spheres represent nitrogen (green), carbon
(magenta), and oxygen (cyan) atoms.

## Structural Transformation in Semiconductor Nanocrystals

Despite numerous experimental studies, structural transformations in SiC
and GaAs at high pressures are not well understood at the atomistic
level. We have investigated the mechanisms of these transformations
using MD simulations. In SiC, a reversible transformation between the
four-fold coordinated zinc-blende structure and the six-fold coordinated
rocksalt structure is found at a pressure of 100 GPa, in good agreement
with experimental data \[24\]. We have found that the atomistic
mechanism for the structural transformation is a cubic-to-monoclinic
unit-cell transformation and a relative shift of Si and C sublattices in
the \[100\] direction. This new transition path does not involve any
bond breaking and it has a significantly lower activation energy
compared with a previously proposed transformation mechanism.

We have also investigated pressure-induced structural transformations in
GaAs nanocrystals of different sizes using MD simulations \[25\].
Semiconductor nanocrystals have numerous applications as optical devices
and new synthetic paths to novel materials. It is found that the
transformation from four-fold (zinc blende) to six-fold (rocksalt)
coordination starts at the surfaces of nanocrystals and proceeds inwards
with increasing pressure. Inequivalent nucleation of the rocksalt phase
at different sites leads to an inhomogeneous deformation of the
nanocrystal. For sufficiently large spherical nanocrystals, this gives
rise to rocksalt structures of different orientations separated by grain
boundaries (Fig. 3). The absence of such grain boundaries in a faceted
nanocrystal of moderate size indicates sensitivity of the transformation
to the initial nanocrystal shape. The pressure corresponding to the
complete transformation increases with the nanocrystal radius and it
approaches the bulk value for a spherical nanocrystal of \~ 5,000 atoms.

![](media/image3.png){width="2.8993055555555554in" height="2.14375in"}

**Figure 3.** (Bottom left) Grain structure in a GaAs nanoparticle at a
pressure of 22.5 GPa is color coded. Insets show grain boundary types.

## Nanoindentation of Silicon Nitride

Nanoindentation testing is a unique probe of mechanical properties of
materials. Typically, an atomic force microscope tip is modified to
indent the surface of a very thin film. The resulting damage is used to
rank the ability of the material to withstand plastic damage against
that of other materials.

We have performed 10-million-atom MD simulations of nanoindentation of
Si~3~N~4~ thin films \[26\]. The films were indented with a square-based
pyramidal indentor to maximum depths of 8-9 nm. The hardness value
calculated from the load-displacement curve at 300K (approximately 50
GPa) is reduced to approximately 90% of its value at 2000K. Figure 4
shows the local pressure distribution directly under the indenter for
the fully loaded and fully unloaded configurations. These images have
been used in conjunction with local bond-angle calculations to
characterize a process of local amorphization under the indenter, which
is arrested by either piling up of material along the indenter edges or
by cracking under the indenter corners.

![](media/image4.wmf){width="3.111111111111111in"
height="1.5763888888888888in"}

**Figure 4.** (a) A half-slice view of pressure in Si~3~N~4~ during
nanoindentation. (b) An atomic view showing amorphization under the
indentor and material pileup at the edges of the indentor. Red and
yellow are silicon and nitrogen atoms, respectively.

## Oxidation of Aluminum Nanoclusters

Oxidation plays a critical role in the performance and durability of
various nanosystems. Oxidation of metallic nanoparticles offers an
interesting possibility of synthesizing nanocomposites with both
metallic and ceramic properties. We have performed the first successful
MD simulation of oxidation of an Al nanoparticle (diameter 200Å, see
Fig. 5) \[27\]. The MD simulations are based on an interaction scheme
developed by Streitz and Mintmire, which can successfully describe a
wide range of physical properties of both metallic and ceramic systems
\[28\]. This scheme is capable of treating bond formation and bond
breakage and changes in charge transfer as the atoms move and their
local environments are altered \[29,30\].

The MD simulations provide detailed picture of the rapid evolution and
culmination of the surface oxide thickness, local stresses, and atomic
diffusivities. Clusters of OAl~4~ coalesce to form a neutral,
percolating tetrahedral network that impedes further intrusion of oxygen
atoms into and of Al atoms out of the nanoparticle. As a result, a
stable oxide scale is formed. Structural analysis reveals a 40Å thick
amorphous oxide scale on the Al nanoparticle. The thickness and
structure of the oxide scale are in accordance with experimental
results.

![](media/image5.png){width="2.1083333333333334in"
height="2.4902777777777776in"}

**Figure 5.** Snapshot of the Al nanocluster after 0.5 ns of simulation
time. (A quarter of the system is cut out to show the
aluminum/aluminum-oxide interface.) The larger spheres correspond to
oxygen and smaller spheres to aluminum; color represents the charge on
an atom.

# Conclusion

Modern MD simulations of materials started in 1964, when Aneesur Rahman
simulated 864 argon atoms on a CDC 3600 computer \[31\]. Assuming a
simple exponential growth, the number of atoms that can be simulated in
a classical MD simulation has doubled every 21 months to reach 1.02
billion atoms \[4\]. Similarly, the number of atoms in a DFT-based MD
simulation (started by Roberto Car and Michele Parrinello in 1985 for 8
Si atoms \[32\]) has doubled every 16 months to reach 22,500 atoms
\[4\]. Petaflop computers anticipated to be built in the next ten years
will help maintain the growth rates in these "MD Moore's Laws" \[4\],
and enable trillion-atom MD and 10-million-atom DFT calculations, see
Fig. 6. (Our latest benchmark on a 1,280-processor IBM SP3 including
8.1-billion-atom MD on silica and 140,000-atom DFT on GaAs indicates an
even faster growth rate.)

![](media/image6.png){width="2.5027777777777778in"
height="2.002083333333333in"}

**Figure 6.** "Moore's law" in molecular dynamics.

In future, metacomputing on a Grid of geographically distributed
supercomputers, mass storage, and virtual environment connected via
high-speed networks will revolutionize computational research by
enabling i) very large-scale computations that are beyond the power of a
single supercomputer, and ii) collaborative, hybrid computations that
integrate distributed, multiple expertise \[33\]. A multidisciplinary
application that will soon require Grid-level computing is emerging at
the forefront of computational science and engineering. We have recently
developed such a multiscale simulation approach which seamlessly
combines continuum mechanics based on the finite-element (FE) method, MD
simulations to describe atomistic processes, and QM calculations based
on the DFT to handle breakage and formation of atomic bonds \[34\].

# Acknowledgments

This work is partially supported by AFOSR, DOE, NASA, NSF, and USC-LSU
MURI. A few million-atom simulations were performed using in-house
parallel computers at the Concurrent Computing Laboratory for Materials
Simulations at Louisiana State University. 10 million to billion atom
simulations were performed using parallel computers at Department of
Defense's Major Shared Resource Centers under a DoD Challenge project.

# References

1\. Pechenik, A., Kalia, R. K., and Vashishta, P., *Computer-Aided
Design of High-Temperature Materials*, Oxford: Oxford Univ. Press, 1999.

2\. Nakano, A., Kalia, R. K., and Vashishta, P., *Comput. Sci. Eng*.
**1** (5), 39-47 (1999).

3\. Kalia, R. K., Campbell, T. J., Chatterjee, A., Nakano, A.,
Vashishta, P., and Ogata, S., *Comput. Phys. Commun*. **128**, 245-259
(2000).

4\. Shimojo, F., Campbell, T. J., Kalia, R. K., Nakano, A., Vashishta,
P., Ogata, S., and Tsuruta, K., *Future Generation Comput. Sys.* **17**,
279-291 (2000).

5\. Hoenberg, P., and Kohn, W., *Phys. Rev*. **136**, B864-B871 (1964).

6\. Vashishta, P., Kalia, R. K., Nakano, A., Li, W., and Ebbsjö, I., in
*Amorphous Insulators and Semiconductors*, edited by M. F. Thorpe and M.
I. Mitkova, Dordrooht: Kluwer, 1996, pp. 151-213.

7\. Greengard, L., and Rokhlin, V., *J. Comput. Phys.* **73**, 325-348
(1987).

8\. Tuckeman, M. E., Berne, B. J., and Martyna, G. J., *J. Chem. Phys.*
**97**, 1990-2001 (1992).

9\. Nakano, A., *Concurrency: Practice and Experience* **11**, 343-353
(1999).

10\. Omeltchenko, A., Campbell, T. J., Kalia, R. K., Liu, X., Nakano,
A., and Vashishta, P., *Comput. Phys. Commun*. **131**, 78-85 (2000).

11\. Chelikowsky, J. R., Saad, Y., Ögüt, S., Vasiliev, I., and
Stathopoulos, A., *Phys. Stat. Sol. (b)* **217**, 173-195 (2000).

12\. Fattebert, J.-L., and Bernholc, J., *Phys. Rev. B* **62**,
1713-1722 (2000).

13\. Troullier, N., and Martins, J. L., *Phys. Rev. B* **43**, 1993-2006
(1991).

14\. Perdew, J. P., Burke, K., and Ernzerhof, M., *Phys. Rev. Lett*.
**77**, 3865-3868 (1996).

15\. Goedecker, S., *Rev. Mod. Phys*. **71**, 1085-1123 (1999).

16\. Mauri, F., and Galli, G., *Phys. Rev. B* **50**, 4316-4326 (1994).

17\. Siegel, R. W., *Scientific American*, Dec. 1996, pp. 74-79.

18\. Kalia, R. K., Nakano, A., Tsuruta, K., and Vashishta, P., *Phys.
Rev. Lett*. **78**, 689-672 (1997).

19\. Kalia, R. K., Nakano, A., Omeltchenko, A., Tsuruta, K., and
Vashishta, P., *Phys. Rev. Lett*. **78**, 2144-2147 (1997).

20\. Tsuruta, K., Nakano, A., Kalia, R. K., and Vashishta, P., *J. Am.
Ceram. S*oc. **81**, 433-436 (1998).

21\. Campbell, T. J., Kalia, R. K., Nakano, A., Shimojo, F., Ogata, S.,
and Vashishta, P., *Phys. Rev. Lett*. **82**, 4018-4021 (1999).

22\. Chatterjee, A., Kalia, R. K., Loong, C.-K., Nakano, A.,
Omeltchenko, A., Tsuruta, K., Vashishta, P., Winterer, M., and Klein,
S., *Appl. Phys. Lett.* **77**, 1132-1134 (2000).

23\. Evans, A. G., *J. Am. Ceram. Soc.* **73**, 187-206 (1990).

24\. Shimojo, F., Ebbsjö, I., Kalia, R. K., Nakano, A., Rino, J. P., and
Vashishta, P., *Phys. Rev. Lett*. **84**, 3338-3341 (2000).

25\. Kodiyalam, S., Kalia, R. K., Kikuchi, H., Nakano, A., Shimojo, F.,
and Vashishta, P., *Phys. Rev. Lett*. **86**, 55-58 (2001).

26\. Walsh, P., Kalia, R. K., Nakano, A., Vashishta, P., and Saini, S.,
*Appl. Phys. Lett*. **77**, 4332-4334 (2000).

27\. Campbell, T. J., Kalia, R. K., Nakano, A., Vashishta, P., Ogata,
S., and Rodgers, S., *Phys. Rev. Lett*. **82**, 4866-4869 (1999).

28\. Streitz, F. H., and Mintmire, J. W., *Phys. Rev. B* **50**,
11996-12003 (1994).

29\. Ogata, S., Iyetomi, H., Tsuruta, K., Shimojo, F., Kalia, R. K.,
Nakano, A., and Vashishta, P., *J. Appl. Phys*. **86**, 3036-3041
(1999).

30\. Ogata, S., Iyetomi, H., Tsuruta, K., Shimojo, F., Nakano, A.,
Kalia, R. K., and Vashishta, P., *J. Appl. Phys*. **88**, 6011-6015
(2000).

31\. Rahman, A., *Phys. Rev*. **136**, A405-A411 (1964).

32\. Car, R., and Parrinello, M., *Phys. Rev. Lett*. **55**, 2471-2474
(1985).

33\. Foster, I., and Kesselman, C., *The Grid: Blueprint for a New
Computing Infrastructure*, San Francisco: Morgan Kaufmann, 1999.

34\. Ogata, S., Lidorikis, E., Shimojo, F., Nakano, A., Vashishta, P.,
and Kalia, R. K., *Comput. Phys. Commun*., to be published.