MATERIALS MODEL-
ING VIA MASSIVELY 
PARALLEL MOLECU-
LAR DYNAMICS

At Los Alamos, we have embarked 
upon a program to study the behavior 
of materials under the influence of 
external stresses and heat flow, using 
the method of molecular dynamics 
(MD), where the motion of millions of 
strongly interacting atoms or molecules 
is followed on a computer. Because the 
complex flow phenomena occur on 
much larger distance scales than do 
molecular sizes and spacings, and 
because the time scales of interest are 
much longer than are vibrational peri-
ods or mean collision times, we are 
using the massively parallel CM-5 at 
the ACL for these simulations. On the 
earlier version of this machine (CM-2), 
we were able to carry out preliminary 
simulations of 106 atoms for over 1000 
vibrational periods, using some 50 
hours of cpu time. For these 2D calcula-
tions, the linear dimension is 0.3 
micrometers, and the physical time 
duration is a nanosecond. In fact, we 
have shown that 8 x 106 atoms are feasi-
ble on the CM-2, and 3 x 107 atoms can 
be simulated on the CM-5, which is 
estimated to be capable of 10 times the 
clock speed of the CM-2 by the end of 
1992. This order of magnitude increase 
in both size and speed of the CM-5 rela-
tive to the CM-2 is comparable to the 
relative advantages of the CM-2 over 
the single-processor Cray Y-MP. On the 
CM-5, we may be able to extend the 
physical time simulated to 0.1 ms, for 
cpu times of less than 200 hours.

 The Fortran 90 MD code on the Con-
nection Machines has been written in 
such a way that a variety of geometries 
(initial conditions) and external driving 
forces (boundary conditions) are easily 
implemented. For bulk materials, 
where periodic boundary conditions 
are appropriate, three distinct kinds of 
dynamics are available: (1) the usual 
Newtonian equations of motion (micro-
canonical ensemble - NVE - constant 
number of particles N, constant volume 
V, constant total energy E, as well as 
constant total linear momentum); (2) 
isothermal Nose'-Hoover MD (canoni-
cal ensemble - NVT - where T is the 
temperature of the heat bath); and (3) 
isothermal-isobaric Nose'-Hoover MD 
(NPT - where P is the diagonal of the 
pressure tensor). In addition to these 
special ensemble dynamics, the peri-
odic boundaries under Newtonian 
mechanics can be moved at constant 
independent velocities. If correspond-
ing velocity gradients are imposed 
upon the peculiar momenta of the 
atoms, homogeneous adiabatic expan-
sion can be simulated; otherwise, inho-
mogeneous shock waves can be 
generated. Moreover, boundary-driven 
flows, such as either Couette shear flow 
or heat flow between hot and cold 
walls, can be simulated. Simple geo-
metric objects (e.g., free surfaces, 
spheres, and plates) can be generated 
and given body velocities. Input to the 
code at setup time and for restarts is 
uncomplicated, and modifications of 
the MD code to do more complicated 
initial and boundary conditions is not 
difficult.

 The particle positions and momenta 
are updated by integrating the equa-
tions of motion via the Stoermer (or 
Verlet) method of finite central differ-
ences. In the case of Nose'-Hoover 
dynamics, the handful of global control 
variables for coupling of the thermostat 
and barostat are also easily handled by 
the Stoermer method. In addition to 
being time-reversible and very robust 
(e.g., energy and linear momentum are 
well conserved over long trajectory 
times), Stoermer integration minimizes 
memory requirements to arrays for the 
phase (particle coordinates and 
momenta) and forces, plus the few sca-
lars for the Nose'-Hoover coupling 
variables.

 We have applied this parallel MD code 
to shock-wave propagation, annealing 
of polycrystalline materials, and spalla-
tion phenomena (fracture). Shock-wave 
simulations revealed similar induced 
plasticity for ductile materials, whose 
atoms interact via a many-body embed-
ded-atom method (EAM) potential 
appropriate for metals, compared with 
brittle pair-potential materials - at least 
in 2D, where dislocations are more eas-
ily formed than in 3D. The annealing 
calculations demonstrate that coarsen-
ing of grains proceeds by the motion of 
grain boundaries in the direction of 
concavity (big grains gobble up little 
ones), so that the average linear dimen-
sion grows with time t like t0.3. Spalla-
tion occurs when a thin plate of 
material hits a thicker target at suffi-
cient velocity. The shock waves gener-
ated on impact are relieved by 
rarefaction waves at the free surfaces 
opposite the impact plane. When the 
rarefactions meet at the spall plane, the 
material is put into tension; if the 
impact is hard enough, the material 
expands until it breaks. Our MD simu-
lations showed that homogeneous 
uniaxial adiabatic expansion (linear 
velocity profile about the spall plane) is 
an excellent approximation. Moreover, 
a simple dislocation-motion model 
accounts well for the stress at failure, as 
a function of the imposed strain rate.

 We are planning to apply this versatile 
MD code (in both 2D and 3D) to a wide 
variety of interesting phenomena, 
including fracture under various exter-
nal driving forces, dislocation genera-
tion at crack tips, shock-wave-induced 
plasticity in granular (polycrystalline) 
solids, ablation of surfaces by radiation 
(in collaboration with 3M), chemically 
reactive flow (as in detonation waves, 
in collaboration with the Naval 
Research Laboratory), and polymer 
dynamics problems, including simula-
tion of the flow of colloidal suspensions 
(in collaboration with du Pont). Many 
of these applications directly impact the 
development and manufacturing of 
novel materials by American industry.

The first sequence (TXY) shows a MD simula-
tion of shock wave in a 2D Lennard-
Jones solid. Particle velocities are col-
ored by a rainbow ranging from -up 
(deep blue) to 0 (blue-green) to +up 
(bright red), where up is the piston 
velocity. The shock front moves to the 
right at the shock velocity (us) which 
is greater than the speed of sound; the 
piston is at the far left. An elastic pre-
cursor has almost reached the right-
hand edge of the window, which 
includes about 1/3 of the total number 
of atoms (65K). The plastic wave is 
near the midpoint of the window.

The second sequence (V) shows a MD shock-
wave simulation. Local crystal orienta-
tion for each particle (hexagonal sym-
metry is assumed) is colored by a 
rainbow ranging from -30 deg. relative to 
the horizontal (deep blue) to 0 deg. (blue-
green) to + 30 deg. (bright red). Shock-
induced dislocation generation is 
apparent near the midpoint of the 
window (plastic wave front); disloca-
tion cores appear as localized spots. 

Acknowledgement: Brad Holian, LANL, T-12