MOLECULAR DYNA.MICS SIMULATION
USING THE MODIFIED EMBEDDED
ATOM METHOD
By
KALYAN K. MAVULETI
Bachelor of Technology
lawaharlal Nehru Technological University
Kakinada, India
1998
Submitted to the Faculty of the
Graduate College of the
Oklahoma State University
in partial fulfillment
the requirements for
the Degree of
Master of Science
December, 2000
MOLECULAR DYNAMICS SIMULATION
USING THE MODIFIED EMBEDDED
ATOM METHOD
Thesis Approved:
TheSiSf#A~
0'rrud' _=_~ _
.(] < 1;._, ~.~ k
11
ACKNOWLEDGEMENTS
First and foremost, I would like to thank my parents for their love, inspiration and
emotional support at every stage of my life to help me grow and attain higher education.
Special thanks are due to my fiancee whose love has provided me with the strength,
confidence and motivation.
I would like to express my sincere appreciation to my advisor. Dr. Ranga
Komanduri, for his constructive guidance, motivation, inspiration and support throughout
my coursework. I would also like to express my sincere appreciation lO Dr. Raff, for his
intelligent advice and help with Molecular Dynamics (M D) and also for the numerous
discussions we had during the present study and \vitlwut wlwsl.' help this project wouldn'l
have succeeded. I would also Ii ke w thank Dr. Price for agreeing to serve on my
committee and for his invaluahlc guidance and advice.
The Molecular Dynamics (MD) code USIng the Modified Embedded Atom
,I\, kthod (ME/\M) was developed and tested with Mr. David Stokes, a graduate student
\\ llLling under the guidance of Dr. Komanduri. Uniaxial tension simulations were
conducted on various metals using this code to verify the validity of the developed
MEAM model for MD simulation, which is the subject of the current report. I extend my
III
gratitude for his valuable assistance and cooperation during this project. ] would also
like to thank Mr. Naga Chandrasekaran, for his advice, encouragement, and valuable
suggestions. Finally, I would like to thank all my friends and colleagues for their support
and encouragement.
IV
TABLE OF CONTENTS
Chapter
I. INTRODUCTION
1.1 Introduction
1.2 Molecular Dynamics Simulation
II. LITERATURE REVIEW
2.1 Introduction to Tensile Testing at Microlevel
2.2 Molecular Dynamics Simulation
2.3 Modified Embedded Atom Method
III. PROBLEM STATEMENT
IV. MOLECULAR DYNAMICS SIMULATON
4.1 Molecular Dynamics Modeling
4.2 Nwnerical Integration
4.3 Choice of the Model
V. POTENTIAL MODEL
5.1 Modified Embedded Atom Method
5.2 Determination of Parameters
5.3 Determination of Forces
VI. VALIDATION OF THE POTENTI AL MODEL
6.1 Numerical vs. Analytical Forces
6.2 Conservation of Lnergy
6.3 Back Integration
VII. RESULTS /\NI) DISCUSSION
7.1 Introduction
7.2 Tension Tests
7.3 Shear Test of Nickel
VI]] CONCLUSIONS
8. I Testing for the Accuracy of the Software
8.2 Uniaxial Tension
8.3 Shear Test of Nickel
8.4 Future Work
REFERENCES
v
Page
1
1
2
5
5
7
11
14
16
16
]8
20
22
22
23
26
27
27
28
33
35
35
35
39
48
48
49
49
50
52
Table
LIST OF TABLES
Page
6.1 Comparison of analytical and numerical forces on atom 1 in xdirection 29
6.2 Comparison of analytical and numerical forces on atom I in ydirection 30
6.3 Comparison of analytical and numerical forces on atom I in zdirection 3 I
6.4 Conservation of energy test 32
7.1 Experimental conditions 35
7.2 Comparison of ultimate tensile strength values obtained by using the 38
MEAM potential with those obtained by Komanduri et al (2000) and
those quoted by Hertzberg (1996)
vi
LIST OF FIGURES
Page
21
32
34
34
Illustration of RungeKutta method in action
Conservation of energy test
Path traced by atom 1 during forward and backward integration
Potential, kinetic and total energies during forward and backward
integration
7.1(a)(d) Snapshots ofthe animation during tension test of nickel using 41
the MEAM potential
Figure
4.1
6.1
6.2
6.3
7.2(a)(d) Snapshots of the animation during tension test of nickel using 42
the Morse potential
7.3
7.4
Tension test ofnicke1 using the Morse pOlential and MLAM
potential
Stress vs. strain curves or niek(:I. copper, and iron using
43
43
the MLt\\tl potential
7 5( a)( d) Snapshots 0 f th~ animation during tension test of copper using 44
the MLAM potential
7J)(aHd) Snapshots of the animation during tension test of iron using
the MEAM potential
7.7 (a)(d) Snapshots of nickel in shear
45
46
VlI
7.8 Comparison of the shear stress/elastic shear modulus VS. shear
Strain for nickel using the MEAM potential with that obtained
by Horstemeyer (1999)
VIII
47
NOMENCLATURE
Ai scaling factor for embedding energy
£;0 sublimation energy (eV)
£101 total energy of the system (eV)
FXI) xcomponent of the force on atom i due to atom j
FYI) yeomponent of the force on atom i due to atom j
FZij .l.component ofthe force on atom i due to atomj
N total number of atoms
R distance between the atoms (A0)
R/ equilibrium nearest neighbor distance (A0)
RI) distance between atoms i and j (A0)
R/ a component of the distance vector between atoms i and j.
V'o, total potential energy of the system (eV)
ZJ number of nearest neighbors
m, mass of atnm i (atomic mass unit)
Pi mOIl1t:ntlll1l or atom i
I), force on atom i (eY/ AD)
(I: Cartesian coordinates of atom i
(i, velocity of atom i (Ao/time unit)
ix
s, (1),0 geometry factors
t/') weighting factors for the atomic densities
Xo initial position of the atom
Xi xcoordinate of atom i
xtJ xcomponent of the distance between the atoms i andj
Xij a ratio of a component of the distance vector between atoms i and j and the distance
between atoms i and j.
y, ycoordinate of atom i
YtJ' ycomponent of the distance between the atoms i and j
Zi zcoordinate of atom i
Zij zcomponent of the distance between the atoms i andj
a; exponential decay factor for the universal energy function
Pj
(I) exponential decay factors for the atomic densities
p, background electron density of atom i
p, total background electron density
P, (I) partial background electron dCIlsit)' nr atom i
x
Chapter 1
INTRODUCTION
1.1 INTRODUCTION
The industry of micromachining is in its infancy today, just as the Very Large
Scale Integrated (VLSI) industry has been developing quickly since the late 70's. As
design tools made the development of the Integrated Circuit (lC) industry possihle,
design tools will make the development ofnew components possible, \'\'hich will cnmhine
the physical world needs of sensing and actuators with tht: rapidly growing Glpabilities of
infonnation technology, Microelectromechanical syslems (MEMS) are miniature
electromechanical sensor and actuator ~ystem~ c.kvelopcd from the mature batchfabricated
processes of VI.Sl technologies. MLMS have wide applications such as
miniature inertial meaSLlTL'mcnl units, biochemical analysis on a chip, arrayed
micromanirulation or parts, optical displays and microprobes for neural recording. The
current 3m] increasing success of MEMS sterns from their promise of better perfonnance.
10\\ manufacturing costs, miniaturization and their capacity for integration with
electronic circuits. The MEMS market is conservatively projected to reach between $12
and $14 billion by the end of this year. Microoptics and MEMS are paving the way for
MicroOptoElectroMechanical Systems (MOEMS). Using MOEMS technology, microoptical
elements are batchfabricated on chips concurrently with microsensors and
microactuators to fonn integrated microsystems. MOEMS teclmology is highly attractive
for commercial applications, since it leverages the integrated circuit infrastructure, which
enables high volume production of microsystem components at a low manufacturing cost.
Uniaxial tension is the most direct way of evaluating mechanical properties of
materials such as the elastic properties, the character and extent of plastic deformation,
yield and tensile strengths, and toughness. Tensile tests are most common in determining
the mechanical properties at macro level. However, when it is applied tu thin film
materials used in MEMS devices, many problems arise:
I. Alignment of the specimen in the testing machine is not easy to perform.
2. Gluing the specimen to the machine is not reliable.
3. Manipulation of the thin film specimen may cause irreparable mcdwllical damage,
4. It is very costly to perfonn nanoregime tensile testing due to colnplcxity of the
equipment.
S, Production of defect free materials in the forlll of lensile specimens is very difficult.
Thus it is very difficult 10 pcrfmm lensile testing of these devices. An alternative
approach would be Molecul<:H Dynamics Simulation, which is very easy to per[onn and is
inexpensl \'C.
2
1.2 MOLECULAR DYNAMICS SIMULATION
The essence of Molecular Dynamics (MO) is simply stated: nwnerically solve the
Nbody problem of classical mechanics. Molecular Dynamics started way back in the
1950's, but widespread attention was given only in the late 1970's. Even then it was
possible only in some of the big national labs which had supercomputers. But today with
the advent of low cost powerful workstations with fast processors (e.g. the Digital Alpha
workstation with 500 MHz clock speed used in this project) and parallel computing, it is
possible to construct large scale MD simulations.
Molecular Dynamics Simulation brings together ideas from several disciplines.
Knowledge of classical mechanics, vector analysis, numerical analysis, thermodynamic:,,
and programming is essential. Also, a good understanding of the manufacturing processes
is required to analyse the results. MO simulation is basically calculating the trajccluries
of the atoms by solving the differential equations of motion. MD predicts the motion of a
given number of atoms governed by their mutual inlcratol11it.: interactions described by a
continuous potential function and requires the numerical integration of Hamilton's
Classical equation of motion. /\ potential model is required to determine the forces on
each atom due to its neighbours. Until recently, the Morse Potential was used to represent
the potential hetw<:<:n tv,o atoms. It is a pair potential and represents FCC materials fairly
well. hut when it comes to BCC materials it does not represent the deformation behaviour
\\cll( Komanduri ~t al (2000)). So, a better potential model is required to represent all the
materials. The Modified Embedded Atom Method (MEAM) developed by Baskes (1992)
:I
represents the material properties better than the Morse Potential and is applicable to
almost the whole range of metals. Pair potentials like the Morse Potential yield the total
energy directly, but need the volume dependant energy to describe the elastic properties
of a metal. If the volume is not represented properly, it may invalidate the results of a
pairpotential calculation because the elastic properties of the solid are not represented
accurately. In the Modified Embedded Atom Method, every atom is considered as an
impurity, embedded in a host lattice consisting of all other atoms. This allows
calculations using electron densities and allows realistic treatment of impurities in
structures that include cracks, surfaces, and alloying additions.
4
Chapter 2
LITERATURE REVIEW
2.1 INTRODUCTION TO TENSILE TESTING AT MICROLEVEL
With the advent of thin film materials like silicon wafers used in the semiconductor
industry, it is becoming increasingly necessary to perform tensile testing at micro level to
determine how the materials would react under load. Sato et al. (1998) performed
uniaxial tensile testing of a singlecrystal silicon film on a silicon chip. They proposed a
tensile testing procedure, in which the external load is applied perpendicular to the
loading lever, by which the file specimen is uniaxially stretched in the horizontal
direction which is given in figures 2.1 and 2.2. This method allowed the tcnsile Icsting of
single crystal silicon film having any arbitrary orientation. Thcy shU\'vcd that the load
linearly increased until the specimen fractured When the fracture occurred, the load
dropped to equal that of the rotatioll3l stiffness or the torsion bars. They performed
uniaxial tensile testing Oil three dilTcrcntly oriented specimens and measured Young's
modulus and fracture strain for each orientation. The measured values were in close
agreemcnt \\ llh Ihe calculated values of bulk materials.
1111: dil\.'cl tension tests, like the one described above, are effective only when
properly performed. The set up requirements for testing, such as alignment and deflection
5
Specimen
Supportmg frame
.Loaamg .lever
TcrSlOn bar
Test material
Figure 2.1 The onchip tensile testing method: chip structure [Sato et al. (J 998)]
.. ~ ... _~..~.
_ ... ..
Spe...."'imen
( Tensile force)
( Elongation)
Torsion bar
Loading point
( Load J
( Displacement )
Lvading lever
( Spring reaction)
Figure 2.2 Th~ onchir tcnsik testing method: cross section
of the chip [Sato et al. (1998)]
measuremcnt. Jrc Jit'Jicult to meet for microscale test samples. Yi et aJ. (2000) proposed
:l nc\\ mcthud in which a load cell measured the induced force and the strain was
m~::J.sureJ by :l laser interferometry system. The advantage of the optical method was that
6
the strain was obtained without physical contact to the sample. The specimen was etched
by four different common silic.on etchants  KOH, EDP, TMAH and XeF2 The Young's
modulus measured, 169.2±3.5 GPa, was very close to the widely accepted value of
168.9 GPa for silicon.
2.2 MOLECULAR DYNAMICS SIMULATION
Eyring and his colleagues (1944) performed the first trajectory calculations of H and
H2 molecules. All the calculations were performed manually since there were no
computers available during that time. They found that the system was caught in a local
minimum and oscillated in that minimum without escaping out. Later. Alder and
Wainwright (1959) devised molecular dynamics in the late 1950's, which is one of the
forms of equilibrium molecular dynamics. It is typically applied to an isolated system
containing a fixed number of molecules N in a fixed volume V. Becausl: th~ system is
isolated, the total energy E (sum of Potential and Kinetic energies) is also C(lnst~nt. The
first applications of MD techniques for molecular simulation were l1l;Jdc ror sinlple
fluids. Another form of molecular dynamics. nonequilihrium molecular dynamics, first
appeared in the early 1970's. In these methods. an eXkrnal force is applied to the system
to establish the nonequilibrium situation of interest, and the system's response to the
force is determined from tht simulation. MD simulation has been applied to various
fields like crystal growth. reactive scattering and simulation of complex liquids in
chemistry. simulations for energetic and structural features of biological systems, and in
7
machining for tension, indentation, cutting and friction at the atomic scale. The available
literature is vast, only the literature concerning the tension and shear is reviewed.
LyndenBell (1994, ]995) investigated the behavior of FCC crystals of the metals
platinum, gold, rhodium and silver under uniaxial tension using the FinnisSinclair
potential. The study was conducted for the variation of potential energy and longitudinal
stress with strain for the above materials at four different temperatures 0.04T~, 0.35 Tm.
0.55 Tm. 0.7 Tm, where Tm was the buJk melting temperature of the material. Void
fonnation and growth of nanocracks were reported which were the causes of failure. The
stress was reported to increase to a maximum at low temperatures and then decreased due
to a series of structural rearrangements. Platinum and gold, which were highly ductile,
were reported to develop local regions of disorder first when compared to the not50ductile
rhodium and silver. However, at temperatures abDve half the melting point, all the
metals were reported tD be disordered before failure, by void formation. For (h4.'
investigation at different temperatures, it appeared that both the shorl and long range
tenns fDr interactions were needed. So, a potential model with manyhouy terms would
have better represented the behavior of the material in buIk.
Rentsch and Inaski (1995) cDnducted MD simulation of silicon under uniaxial tension
using the TersoffpDtential. They rerorted a linear stressstrain relationship followed by a
sudden break down to zero. Anisotropic defonnation of silicon was reported. The
Young's muuulus was found to be 171 GPa and the specific surface energy to be 0.393
.rIII .: .
8
Kitamura et al (1997) perfonned MD simulation of nanDsingle crystal of nickel in
tension using the embedded atom potential. Tension tests were carried out under two
conditions: (1) tension without constraint of transverse deformation and (2) tension with
constraint of transverse deformation. In the first case, yielding was brought about by the
crystallographic slip on the {Ill} planes at a strain of 0.1. The yield stress in tension was
about] 5  20 OPa and very little differences were noticed among the wire, film and bulk
samples. The multiple slip on the {] II} planes continued to take place after the yield.
The plastic deformation caused ductile shear fracture. With constraints, the yield stress
reached 40 OPa. No plastic strain was generated. A cleavage crack initiated and brought
about brittle fracture. It was reported that the constraint changed the fracture mode.
Heino et al (1998) studied the mechanical properties of copper by MD using the
effective medium theory as the potential model. Simulations of point defects, grain
boundary, and a larger void, which served as the seed for crack propagation, were studied
at room temperature. A decrease in fracture stress and strain, and tensile modulus was
reported with an increasing number of defects. Systems with larger nUlllhcr of dckets
were reported to appear more isotro ic that ordered systems, in terms of tensile modulus.
They reported that the systems with grain boundaries were weaker than the ordered
systems in terms of modulus, fracture stress and fracture strain. With thick systems, with
free boundaries and an initial large void, the {Ill} slip plane was reponed.1.o opropagate
in a (110) direction with a speed of about 60% of the longitudinal speed of the sound for
the speci fic crystal orientatio~..~ith thin systems, including a crack seed and having free
boundaries, crack propagation was reported m the (II~) direction by microvoid
coalescence.
Komanduri et al (2000) studied the uniaxial tension of singlecrystal materials, both
FCC [AI, Cu, and Ni] and BCC [Fe, Cr, and W] using the Morse Potential. They reported
a rapid increase in stress up to a maximum followed by a ~radual drop to zero when the
specimen failed by ductile fracture. They also reported that the radius of the neck
increased with an increase in the .defo:m:ation of thespedmen and decreased as the
ductility of the material decrease~:..g~pid fluctuations in the force values were reported.
The strain to fracture was reported to be lower for BCC materials than FCC materials.
Tungsten had the highest strength and aluminum had the lowest strength. The ultimate
~
tensile strength of Cu, Ni, AI, Fe, Cr, and W were reported to be 28, 36, 13,29, 31, and
51 GPa, respectively. The strain to fracture was reported to be 2.17, 1.67, 3.2, 1.52.
1.52, and 1.4 for Cu, Ni, AI, Fe, Cr, and W, respectively. They also reported a good
correlation between the D~ and a parameters of the Morse potential with the ultimate
tensile strength and the strain to fracture for the FCC materials, and no such correlation
for the Bee materials.t:he~ggested that an alternate potential model should he used
for BCC materials since the deformation behavior wasn't represented well by the Morse
potential. ,! '
Horstemeyer and Baskes (J 999) performed atomistic finite deformation using the
Embedded Atom Method. They observed a spatial size scale effect on the yield stress.
They ohservcd that the mechanical yield point occurred from dislocation initiation at the
10
edge of the numerical specimens. They also observed that~ as the spatial length scale
increased, the continuum rotational effect coupled with the increase in the dislocation /
I
population reduced the oscillatory behavior. They proposed a length scale bridging idea ,
i.
by relating a continuum single degree of freedom loss coefficient, which related the _J
plastic energy to the total strain energy, to varying sizes of blocks of atoms.
By the above review, it is clear that a potential model which is good for all the metals,
FCC, BCC, diamond structure and HCP, must be used for MD simulation. Modified
Embedded Atom Method (MEAM) is one such potential model and is considered in this
project.
2.3 MODIFIED EMBEDDED ATOM METHOD
The first step towards the present day MEAM was the quasiatom theory proposed
by Stott and Zaremba (1982), that was used successfully to calculate the characteristics of
hydrogen in metals. They proposed that "the energy of an impurity in a host is a function
of the electron density of the unperturbed (i.e. without irppurity) ho~t".
Daw and Baskes (1984) generalized the quasiatom theory to treat all atoms in a
unified way, and called it the embedded atom method. hey proposed that. every atom is
considered as an impurity, embedded in a host lattice consisting of all other atoms. This
allowed calculations using electron densities and allowed realistic treatment of impurities
in structures that include cracks, surfaces. and alloying additions. The electron densities
II

were approximated by the linear superposition of spherically averaged atomic electron
densities. They also proposed that the embedding energy of the atom depends only on the
environment immediately around the impurity, i.e. impurity experiences a locally
uniform electron density. The energy was given by
(1)
Where Pi was the electron density of the host atom without atom i, and ~ was the shortrange
electrostatic pair potential. From this expression for the total energy of a given
metal, several ground state properties like lattice constant, elastic constants, sublimation
energy, and vacancyformation energy were calculated. The validity of the above
functions was tested by computing a wide range of properties, like the formation volume
and migration energies of vacancies, the formation energy, the migration energy of
divacancies and selfinterstitials, the surface energy and geometries of the lowindex
surfaces of the pure metals, and the segregation energy of substantial impurities to {I OO}
surfaces. The embedded atom method developed by Daw and Baskes was based on the
density functional theory and the electron density was approximated as a linear
superposition of electron densities. These assumptions are better approximations for FCC ._.
metals but not for BCC metals.
Adams and Foiles (1990) extended the embedded atom method to BCe material
Vanadium. Since the electron density in BCC metals was not well approximated by linear
superposition of electron densities, the authors used the adj uslahle electron density
proposed by Voter and Chen (1987). This theory was adopted since spherically
12
symmetric electron densities were easier to incorporate in the model. The pair tenn in Eq
.~...  ~ ~••  . _ ......... _'0"_ '.' I
(1) was asswned to have the fonn of the Morse Potential.
Later on, Baskes (1987) modified the embedded atom method to include
directional bonding and applied it to silicon. Baskes, Nelson and Wright (1989) extended
the silicon embedded atom method to the silicongermanium system. This extended
method was called modified embedded atom method. This method had a few deficiencies
when it was initially developed. There was inward relaxation at a vacancy, an extremely
large stacking fault energy and only qualitatively accurate small cluster predictions.
These deficiencies were partially resolved in the later paper by Baskes (1992). The
common attribute of all the papers is that the interaction between two atoms depends on
......  •• &'
the local environment.
The MEAM has been applied to metals and semiconductors and also for diatomic
gaseous elements. In this method, simplification to the first nearest neighbors is possible,
which reduces the computational time. The difference between the EAM and MEAM is
..._.. j ....... _il ••.. ~ ......
that the .Pi. which is. give~ as the line~rly superposition of spherically averaged atomic
electron densities in EAM is augmented by an angularly dependent term in MEAM. . . "." . •  """Ia ~~~..."._ ... to... ,, ••• •
Baskes and Nelson (1994) have extended MEAM to HC? materials. So, MEAM is much
more versatile and can be used for FCC, BCC, diamond structures, and Her metals.
I ''
Chapter 3
PROBLEM STATEMENT
MD simulation has been used to conduct uniaxiaJ tension of different materials using
different potential models like the FinnisSinclair potential (LyndenBell [1994,1995]),
the Tersoff potential (Rentsch and Inaski [1995]), the Morse potential (Komanduri et al.
[2000]), and the effective medium theory (Heino et al. [1998]). Using the FinnisSinc1air
potential, it appeared that a potential model with many body terms would have better
represented the behavior of the material in bulk. When the Morse potential was used for
MD simulations, good correlation was found between the D and a parameters of the
Morse potential with the ultimate tensile strength and strain to fracture: for the FCC
materials. No such correlation was found for BCe materials. So, a potential model which
can represent the bulk properties of metals fairly accurately, and which can be used for
the whole range of metals should be used for MD simulations. The embedded atom
method and the MEAM are such potential models. The difference between these models
is that the linear superposition of atomic electron densities in the embedded atom method
is augmented by an angularly dependent term in the MEAM. Thus MEAM polelltial j"
chosen as the potential model for MD simulation. The objcctuvc (d' thi" slUdy i~ La:
11
1. Develop the software for MD usmg the MEAM as the potential model for the
trajectory calculations
:2. To find the forces on the atoms by calculating the derivatives of the total potential
with respect to the three coordinate axes, x, y and z..
3. To validate the software by different testing procedures like the numerical vs. the
analytical force test, the conservation of energy test and the back integration test, to
validate the accuracy of the model.
4. To perform a shear test of Nickel to evaluate the shear stress and strain and to observe
the deformation during the simulation.
5. To perfonn uniaxial tension calculations using MD for various FCC and Bee metals
and to find the ultimate tensile strength and the strain to fracture of these metals and,
also, to observe the deformation and necking during the simulation.
Chapter 4
MOLECULAR DYNAMICS SIMULATION
Molecular dynamics simulation basically involves the calculation of trajectories.
This calculation involves the numerical integration of classical equations of motion for a
system of interacting atoms over a period of time. The time step used in this integration is
of the order of 1015 sec, which is less than the period of vibration of the atoms.
4.1 MOLECULAR DYNAMICS MODELING
Consider an isolated system comprising N bodies with the coordinates (Xi, Yi, lj)
where i = 1, 2, 3, .,.",N. Given a set of N independent generalized coordinates and
velocities {q"q, }that describe the state of a conservative system (one in which all the
forces derive from some potential energy function U ), so that L =L({q, }, {qJ }, t) ,then L
can be shown to satisfy the Lagrange equations
~(8L]_ ~ = O,i =1,,,.,,N
dt 8q, 8q,
(4.1 )
These equations are the starting point for many of the suoscquent dcvelopments,
Newton's second law is a simple consequence or this result. \\here. if qi denotes a
II)

component of the Cartesian coordinates for one of the atoms (and assuming identical
masses m):
So that equation (4.1) becomes
mq, =av/aq, = F,
Where Fi is the corresponding force component.
(4.2)
(4.3)
By replacing the generalized velocities {q in the Lagrange formulation by the
generalized momenta Pi = ai/aq, (if the coordinates are Cartesian, thenp, = mqi) and
consider the Hamiltonian H = H({qj}' {P, }, t) defined by
The two first order equations of motion associated with each coordinate are
. aH
p=
I aq,
(4.4)
(4.5)
(4.6)
If H has no explicit time dependence, then H =0, and H  the total energy  is a
conserved quantity.
From the above equations (4.5) and (4.6), we get the following dillcrcntial
equations for each atom in the three coordinate systems.
17
BHjBPx, =dx;/dt = pXj 1m,
BHjBPy; =dy;/dt =py, 1m,
8HjBPz, =dz;/dt = pz, Imj
8HjOx, =8V/Ox; = dpx;/dt
BHjBy, =8V/Byj =dpy)dt
8Hj8z; = 8V/8z, = dpz)dt
(4.7)
(4.8)
(4.9)
(4.10)
(4.11 )
(4.12)
So the total number of differential equations to be solved are 6N. The following are the
system of units used in the above equations
1 mass unit = 1 atomic mass unit = 1.007/6.023 x 1023
1 distance unit = 1 A0 = 10·gcm
1 energy unit = 1 eV = 23.06 kcaVmol = 4.184 x 23.06 kllmol
1 time unit = 1 t.u. = 1.018 X 10 14 sec.
4.2 NUMERICAL INTEGRATION
The calculation of trajectories requires the numerical integration of the equations
(4.7) through (4.12) from an initial state in the configuration space identified as reactants
to some final state associated with products. The various numerical techniques used are:
1. Fourth order RungeKutta method (which is self starting)
2. Fifth, sixth etc order predictorcorrector methods (non self starting)
3. Variable step size methods.
RungeKutta method has several advantages when compared to the other methods.
IR
1. It is selfstarting, so it is unnecessary to know the values of elements prior to t = 1{).
2. The integration error is very small, of the fourth order. 0 (h4). So the error can be
neglected without affecting the results significantly.
3. The method is stable and fairly easy to program.
Its main disadvantage is:
1. The need to compute a large number of derivatives (24N) for each integration step,
which demands a lot of computer time
Though this is a disadvantage, it is worth it because of the greater accuracy and the
stability of the method. Other methods, like PredictorCorrector methods, have the
advantage of providing an automatic error estimate at each integration step, thus allowing
the program to use variable step size to achieve the specified accuracy. However, these
methods are not selfstarting and require the use of RungeKutta method to start the
integration. So, the RungeKutta method was chosen for the calculations.
If there are N particles whose initial derivatives 01 = f (xo. Yo) are known and we
integrate them over a period of time t, the Runge·Kutta method involves the calculation
by the following steps.
1. Move them to a new position P2 USIng Oland a time step of 1/2, so that
D2 = f (xo+~x12. yo+D1" ~x/2)
2. Get them back to their original positions and move them to their new positions P3
using the derivatives 02 and calculate 03 = f (xo+~x12. Yo+02* ~x/2)
3. Get them back to their original positions and move them to their nc\\ positions ha~eJ
l1n derivatives 03 and a time step oft and calculate 04 = f (xrf'llx.i:2 )q+f)V"l\\/2)
19

4. Compute the average derivative <D> = [01 + D2"'2 +03"'2 + 04]/6.0 and then
compute the !1y = <D>*!1x and y 1 = Yo +!1y , to get them to their final positions
Figure 4.1 illustrates this point. This is a simultaneous approach in the sense that all the
particles are moved during the above steps. Thus, for a set of N particles, there will be 4N
intermediately calculated positions and forces for a RungeKutta procedure time step.
The interatomic potential used in the simulation to model the materials plays an
important role in determining the accuracy of simulation results. The modified embedded
atom method is used in the present calculations, which is explained in the next chapter.
4.3 CHOICE OF THE MODEL
The following are some of the important considerations in the choice of the model
for MD simulations:
1. The number of differential equations to be evaluated for a system of N atoms, which
is 6N. So, if there are 1000 atoms, we need to integrate 6000 first order coupled
differential equations. This provides a restriction on the computer memory and
computational time. So, we need to choose a model within these restrictions
2. Number of terms in the potential energy hypersurface that are to be considered. For
example, in a pair potential it is N(Nl )/2 and for a many body potential like the
Modified embedded atom method, it is N(Nl)
3. Size convergence should be considered so that the results are independent or our
choice ofN.
• 20
PI
P3
Final Position
O..'
...,."
.............
,;.......
./
.......
.................
.....
02
P2
D4
P4
Figure 4. J Illustration of RungeKutta Method in action
21
Chapter 5
POTENTIAL MODEL
The choice of the potential model to be used in the trajectory calculation is very
important since it has to represent the properties of materials in bulk. The modified
embedded atom method is a good choice in the sense that in represents the properties of
bulk materials better than other pair potentials, like Morse or Lennard Jones potentials.
5.1 MODIFIED EMBEDDED ATOM METHOD
The total energy £ of a system of atoms given by the embedded atom method is:
I
£/u/ = L F; (P, (R;)) +  L ¢J(R,1 )
I 2'J·
(5.1 )
The first term in the right side of the equation is the embedding function, i.e., the energy
to embed an atom of type i into the background electron density at site i, P" and the
second term is the pair interaction between the atoms i and j, whose separation is given
by Ri}. According to the embedded atom method, P; is the linear supposition of
spherically averaged atomic electron densities. In the modified embedded mdhod. it is
augmented by an angularly dependent term. The energy for an individual a10111 is given
22
(5.2)
The background electron density is renonnalized by dividing it by the number of nearest
neighbors Z;.
The pair interaction is given by the equation
¢;; (R) =2 {E/ (R)  F, (Pi O(R) I Zj)}
Zj
From the above three equations the energy of an atom Ej is given by
(5.3)
(5.4)
IJ· ,j \ I ,,
The first part of the above equation is the average of the energy per atom, of the reference
lattice at each of the nearestneighbor distances. The second part is the difference
between the embedding energy at the background electron density actually seen by the
atom i and the average embedding energy of this atom in the reference lattice, at each of
the nearestneighbor distances.
5.2 DETERMINATION OF PARAMETERS
The equation (5,3) has three tenns and each of the tenns can be detennined in the
following way. In the first term Zj is the number of nearest neighbors and Et is given by
o a*=a,(RIRj I)
Where
23
(5.5)
(5.6)
ai =exponential decay factor for the universal energy function 7
Et = sublimation energy (eV)
R = distance between the atoms (A~
RiD = equilibrium nearest neighbor distance (A0)
In the second term, the function F is given by
So the second part becomes
Where
Pi = total background electron density
Ai = scaling factor for embedding energy
Ei
D
= sublimation energy (eV)
The total background electron density is given by
J
(Pi)2 ="Lt/'\p//)2
1=0
Where
f/I) = weighting factors for the atomic densities
The first partial background electron density at site i is given by
P,(0) = 'L"J'P 0(0) (R ) J IJ
j(.,i)
(5.7)
(5.8)
(5.9)
(5.10)
Where the atomic electron density oftypej atom at a distance Rij from site i is given by
0(1) (R ) = h' Pj ij e
eh' = f3 (/)(R / RO 1)
) IJ I
24
(5.11)
(5. 12)
f3j (I) =exponential decay factors for the atomic densities
RIO =equilibriwn nearest distance (A0)
RIj =distance between atoms i and} (Ao)
Similarly the second, third and fourth partial electron densities are given by
2 2
(Pi (2) ) 2 = "L.[L".XijaXij fJ Pu(2) (Rij )]  31"[L".PjU(2)(Rij )]
a.fJ K"i) )(,<ij
(5.13)
(5.14)
(5.15)
Where Xij a =Rij a / Rij , and Rij ais the a component of the distance vector between atoms
i and j. The above equations are chosen so that the partial background electron densities
are invariant to lattice translation and rotation, scale simply with atomic electron density
for homogeneous deformation, and equal zero for a cubic lattice,
In the third term, the function F is given by
Where
1
(PO, (RIJ))2 =L~t.(/I ) S,(/),o(P,.U(l)(R'I ))2
I~O
Sf (1),0 =geometry factors
The rest of the terms are the same as explained above.
25
(5.16)
(5.17)

All the prior tenns are put together to obtain the final equation, i.e., the total potential of
the system. The next step is to detennine the forces between the atoms from the potential
function.
5.3 DETERMINATION OF FORCES
The detennination of forces is crucial and time consummg because of the
complexity of the potential function. The force on an atom i due to atom} is obtained by
differentiating the total potential with respect to x, y, and zcomponents of the distance
between the atoms i and}, to get the force in x, y and z directions respectively.
Force in xdirection,
Fx. = _ av,o/
IJ Ox ..
I)
Similarly, in y and zdirections
F = _ av/o/
Yij rh'
VYlj
(5.18)
(5.19)
(5.20)
Once the potential and forces are obtained, the model is to be tested to validate the
accuracy of the model.
26
Chapter 6'
VALIDATION OF THE POTENTIAL MODEL
Once the software for the MD simulation is developed, i.e., the potential and the
forces are determined; it has to be tested before it is used. There are many ways of testing
the model for different parameters like numerical vs. analytical forces test for the
validation of the forces, the conservation of energy test and the back integration test.
6.1 NUMERICAL vs. ANALYTICAL FORCES
The first and foremost thing is to validate the force function. This is done by
comparing the analytical forces, got from the derivatives of the potential, to the numerical
forces derived from the formula:
dVI de f~'n =0.75 x 5 1 0.15 x 5 2 + 0.01666666667 x 53
51 = [V(xo+ Lli)  (xo Lli)]/ Lli
52 = [V(Xo + 2Lli)  (Xo  2Lli )]/ ~x
53 = [V(Xo+ 3Lli) (xo  3Lli)]/ Lli
27
(6.6)
(6.7)
(6.8)
(6.9)
The testing sample contains five atoms, which are placed in the following configuration
(2.0,2.0,0.0), (2.0, 2.0, 0.0), (2.0, 2.0, 0.0), (2.0, 2.0,0.0), (0.0, 0.0, 2.0). The atom 1
is moved to 1 AO in steps of 0.1 AO and the forces are calculated at each step. The
numerical and analytical forces agree well, demonstrating that the analytical forces are
right. The tables 6.1 through 6.3 show the comparison between analytical vs. numerical
forces in x, y and z directions respectively
6.2 CONSERVATION OF ENERGY
After the validation of the force function, it is necessary to check if the system
conserves energy. If there are no external forces acting on the system, and if the model is
allowed to integrate for the given period of time, the sum of potential and kinetic energies
which is the total energy must remain constant. This test is necessary, as it not only
validates the accuracy of the potential and forces functions but also the integration
procedure, i.e., the RungeKutta procedure.
For this test, the sample is one lattice of Nickel, which has 14 atoms. This sample
is allowed to integrate for 10 time units and the values are given in the table 6.4. The total
potential is constant for up to 12 significant digits. The same is shown in the figure 6.1
28
:l
~
cr' ~
0\ n0a
'0
~
:=!.
Vl
0
;:l
0.......,
~
;:l
~
'<.....
(i"
~
b
0
N ;:l
'D t:
3
~
:J.
()
~
0...,'
()
~
til
0
::l
~.....
0a
::l
;..<
I
0 ...,
<1l
~
0
::J
Distance moved by Analytical force on I Numerical force on Analytical force on Numerical force on
I
atom 1 in the atom 1 I atom 1 atom 2 atom 2
positive x direction in the xdirection I in the xdirection in the xdirection in the xdirection
0 0.317994486191 j 0.317994486191 ,i 0.1462990104611 0.146299010461
o1 0392243426220 0.392243426220 0.1502041368581 0.150204136858
02 0443425449382 0443425449382 0.156172952299 1 0.156172952299
o3 0478482723558 0.478482723558 0.163469672011 ; 0.163469672011
0.4 0501392714854 0.501392714854 0.171538306544 0.171538306544
05 0514487845642 0.514487845642 0.179955816315 0.179955816315
06 0519309615179 0519300615179 0188397761492 0.188397761492
0.7 0517073341564 0.517073341564 0.196614201357 0.196614201357
0.8 0.508875262056 0.508875262056 0.204413074192 0.204413074192
0.9 0495757231715 0.495757231715 O.211648656965! 0.211648656965
1 0478706598832 0.478706598832 0.218213345221 ' 0.218213345221
Distance moved by Analytical force on Numerical force on Analytical force on
,
Numerical force on
atom 1 in the atom 3 atom 3 atom 4 atom 4
positive x direction in the xdirection in the xdirection in the xdirection in the xdirection
0 0257069230915 0257069230915 0236635445138' 0,236635445138
01 0.338090086813 0338090086813 0,232026780657 0.232026780657
0.2 0 396303411500 0.396303411500 0.228303822004 0,228303822004
0.3 0437962686834 0.437962686834 0.225355906459 0225355906459
04 0.466695160543 0.466695160543 0.223087644108 0.223087644108
05 0.484697915969 0484697915969 0.221418205942 0.221418205942
0.6 0493505232742 0493505232742 0220279968193 0.220279968193
07 0.494398598502 0494398598501 0.219616899031 0.219616899031
08 0488581739721 0.488581739721 0.219382936439 ' 0.219382936440
0.9 0.477227183899 0.477227183900 0.219540480476' 0.219540480476
1 0.461465660331 0,461465660331 0.220059039661 1 0.220059039661
,
....,
~
0" n
O'.
........
na3
"0
~
:J.
v a
;:J
a
'"""T'.
t>l
;:J
~ '..<..
f)'
2..
[3
0.
w ;:J 0 c3
n
:J.
("l e:
.O.,'
("l
n
V>
0
;:J
.t>..l.
03
>
;:J
'<:
I
.0.,.
n
("l .. o'
;:J
Distance moved by I Analytical force on II Numerical force on Analytical force on Numerical force on
I , I
atom 1 in the I atom 1 atom 1 atom 2 atom 2
positive x direction! in the ydirection in the ydirection I in the ydirection in the ydirection
Oi 0054107368965 0.054107368965 0.080636084097 0.080636084097
o 1 0041623846943 0041623846942 0075181102209 1 0.075181102209
02 0 025029638329 0025029638329 0.066073809105 0.066073809105
03 0005801407510 0005801407510 o054151400966 0.054151400966
04 0.014951936750 0.014951936750 0.040057342538 0.040057342538
05 0036413427791 0.036413427791 0.024293594021 0.024293594021
0.6 0057983435612 0.057983435612 0.007260434499 0007260434499
0.7 0079217767652 0079217767652 0010714927052 0.010714927052
08 0099781729090 0.099781729090 0 029358504450 0.029358504450
09 o 119418518149 0.119418518149 0.048435958660 0.048435958660
1 o 137928567431 0.137928567431 0.067743288466 i 0.067743288465
,
Distance moved by Analytical force on Numerical force on Analytical force on Numerical force on
atom 1 in the atom 3 atom 3 atom 4 atom 4
positive x direction in the ydirection in the ydirection in the ydirection in the yclirectjon
0 0051063618610 0051063618610 0019561987246 0.019561987246
0.1 0040987459133 0 040987459133 0021888614212 0.021888614212
02 0028741038951 0028741038951 0.026941023364 0.026941023364
0.3 0015132293624 0015132293624 0.033992787708 0.033992787709
04 0000657430361 0000657430360 0.042541731626 0.042541731626
0.5 0.014379828210 0.014379828210 0.052235196542 0.052235196542
06 0029786079394 0.029786079394 0.062817203350 0.062817203350
0.7 0.045428224009 0.045428224009 0074094373119 0.074094373119
0.8 0.061206002845 0061206002845 0.085915200241 0.085915200241
0.9 0077038261115 0077038261115 0098157757934 0.098157757934
1 0.092856669310 0092856669310 0 110722307002 0.110722307002
vl
.....,
I>l cr
~
0'
(J
o3
'0
p)
~.
(fJ o
;:l
o......,
I>l
;:l e:..
'<......
n
p)
p)
:::l
0..
;:l3~
::1.
ne:.
O'
""1 n
~
V>
o
;:l
p)o3
s·
Ne:
~n
o
:::l
Distance moved by i
atom 1 in the I
positive x direction i
o
0.1
02
03
04
05
06
07
0.8
09
1
Distance moved by
atom 1 in the
positive x direction
o
01
0.2
03
04
0.5
0.6
0.7
08
09
1
Analytical force on I
atom 1
in the zdirection
0.119698158477·
0.127825866682
0 136446272533
0144953537850
0152895817351
0 159959002693
0 165938142358
0.170708506925
0 174202060831
0 176390634542
0.177274897688
Analytical furce on
atom 3
in the zdirection
0.046221264055
0.050412089960
0054611433443
0058718867904
0062715834715
0066624396184
0.070480257711
0.074318051068
0078164992095
0.082039377610
0 085951445681
Numerical force on
atom 1
in the zdirection
0.119698158477
0.127825866682
0 136416272533
0 144953537850
0.152895817350
0 159959002693
0 165938142358
0 170708506925
0.174202060831
0 176390634542
0.177274897688
Numerical force on
atom 3
in the zdirection
0046221264055
0050412089960 ;
0 054611433442
0.058718867904
0062715834715
0.066624396184
0070480257711
0074318051068
0078164992095
0082039377610
0.085951445681
Analytical force on
atom 2
in the zdirection
0.113413198955:
0.111954133477 1
0.110826785940
0 110028566566
0 109543211377
0109346965477
0.109413322326
0.109715992768 .
0.110230426767
0.110934382926
0.111807971497
Analytical force on
atom 4
in the zdirection
0.096603114301
0.095602641243 !
0.094868811573
0.094418927064 ,
0 094258323386
0.094384779920
0094792331587 1
0.095473873734 '
0096422598917
0.0976325665561
0099098730598
Numerical force on
atom 2
in the zdirection
0.113413198955
0.111954133477
0.110826785940
0 110028566566
0 109543211377
0.109346965477
0 109413322326
0.109715992768
0.110230426767
0.110934382927
0.111807971497
Numerical force on
atom 4
in the zdirection
0.096603114301
0.095602641243
0094868811573
0.094418927064
0.094258323386
0 094384779920
0.094792331587
0095473873734
0096422598916
..() 097632566556
0.099098730598
Table 6.4 Conservation of energy test
Conaarvatlon 01 Enargy Taat
10.00 . ,
500
000
5.00
'>
j ·10.00
~
.15.00
Total Energy
20.00
2500
Potential Energy
30.00 ' '
Figure 6.1 Conservation of energy test
32
6.3 BACK INTEGRATION
Back integration is the most sensitive of all the tests. In this test, the model, which
is an isolated system, is allowed to integrate over a period of time, and then the time step
is made negative, and the model is allowed to back integrate for the same period of time.
The model should trace the same potential, kinetic and total energy curves in both the
forward and backward integration. Also the atoms must trace the same path in both the
forward and backward integration.
For this test, one Lattice of Nickel is taken and allowed to integrate for a period of
100 time steps. Then the time step is made negative and the structure is allowed to back
integrate for the same period of time. The potential, kinetic and total energy curves for
both fonvard and backward integration are given in the figure 6.3. The curves overlap
each other almost exactly and look as if they are the same curve. Also the position of one
of the atoms in the structure is plotted against time and it traces the same path in both the
directions and can be seen in the figure 6.2. These three tests validate the potential and
force functions and also the integration procedure.
33
!
I
10 60 /0 10
Forward Integration
Backward Integrati
30
o.
o.
o.
~
g 02
II
'll
.;,]
0
'0
02
0" ; \\, J '~
00
11rrw shipja.I.u.J
Figure 6.2 Path traced by atom 1 during forward and backward
integration
'0 10 110
Forward Integration
Backward Integration
'0
Total Energy
'cr
I
·'0II
."
·20
."
30
Tn. .lap (&t.ul
Figure 6.3 Potential, kinetic and total energies during forward and backward
integration
34
Chapter 7
RESULTS AND DISCUSSION
7.1 INTRODUCTION
The potential model developed is good for any material with FCC, BCC, diamond
structure, and HCP metals. The values are given for some of these materials in the paper
by Baskes (1992). Of these, nickel, copper, and iron are selected and subjected to uniaxial
tension.
7.2 TENSION TESTS
Table 7.1 shows the conditions in which the experiments are conducted.
Configuration 3 Dimensional
Potentials Used Morse, MEAM
Wark Material Dimension 4a x 4a x 6a (a = lattice constant)
Tensile Loading Condition Uniaxial
Tension Rate 500 m/s
Bulk Temperature 293 K
3S
The specimen has boundary (which are fIxed) atoms at the top and bottom surfaces.
Immediately after the boundary atoms, is a layer of peripheral atoms. The rest are moving
atoms. The tension simulations are conducted at the rate of 500m/s to achieve reasonable
computational time. Consequently the system temperature will increase significantly
which is dissipated by means of the peripheral atoms. The motion of the atoms in the
moving zone is determined solely by the forces produced by the interaction potential and
the direct solution of the classical Hamiltonian equations of motion. The motion of the
peripheral atoms is also calculated from the solution of Hamiltonian equations, but
modifIed by the presence of velocity reset functions associated with each atom in the
peripheral zone. In this method, the Cartesian velocity components of each peripheral
lattice atoms is reset at periodic time intervals, !:!.t, using the following algorithm:
new =(1 )112 old 112V(T J=) Val W Val +W .,=, (8.1)
Where Vaoildl.S the acomponent (a=x, y, or z) of velocity of lattice atom
resulting from the solution of the Hamiltonian equations of motion, and va/
cw is the reset
a velocity component. 'w' is a parameter that controls the strength of the reset with w=o
corresponding to no reset and w= I being a complete reset. V (T, ~) is a randomly chosen
velocity from a MaxwellBoltzmann distribution at temperature T. ~ is a random number
whose distribution is uniform on the interval [0,1] that controls the random selection.
This procedure simulates the thermostatic effect of the bulk and guarantees that the
equilibrium temperature will approach the desired value, which is 293 K in these
calculations.
~ 36
The animations for nickel using the MEAM and the Morse Potentials are given in
figures 7.l(a)(d) and 7.2(a)(d). After relaxation, a light bulge in the specimen is
observed, using the Morse potential as well as the MEAM potential. This is because,
when the crystal relaxes, it tries to take the minimum energy position and tries to attain a
spherical shape. But, since it is constrained at the top and bottom layers by boundary
atoms, it bulges only on the +ve and ve x and ydirections. The bulge is more in the case
of the Morse potential than the MEAM potential. MEAM potential represents crystal
surfaces better than the Morse potential. This is because MEAM potential involves
calculations using electron densities and allows realistic treatment of impurities in
structures that include cracks, surfaces, and alloying additions. When the tension test
begins, similar behavior is observed in both the cases. The bulge decreases gradually and
the specimen starts to neck. Because of the high tension rates and also because the
specimen is very small, the crystal becomes amorphous almost immediately after the
experiment begins. The necking continues and the diameter of the neck decreases with
the increasing strain. This process continues until the specimen fails due to fracture. The
stressstrain curves for the tension test of nickel using the Morse and the MEAM
potentials, is given in the figure 7.3. The curves behave almost similarly until reaching
the peak of the tension curves. After that, the tension curve of nickel, using the MEAM
potential drops rapidly when compared to the one using the Morse potential. This can be
attributed to the fact that the atoms in the tension test using the Morse potential are
bonded by a pairwise potential and these bonds exist until the all the bonded pairs in the
center of the specimen are out of the cutoff radius. But, in the case of the MEAM
potential, the volume of the material is represented well, because the calculations are
37
done using the electron densities. This results in the rapid drop down in the tension curve
simulating the behavior of the bulk material. Using the MEAM potential, the ultjmate
tensile strength is 44.46 OPa approximately at a strain of 0.31 and in the case of Morse
Potential the ultimate tensile strength is 54.25 OPa approximately at a strain of 0.232.
The stressstrain curves using the MEAM potential for nickel, copper, and iron
are shown in figure 7.4. The animations for the tension tests of copper, and iron using the
MEAM potential, are given in figures 7.5(a)(d) and 7.6(a)(d). For FCC materials, nickel
and copper, the simulations show reasonable behavior. Copper has a lesser ultimate
tensile strength than nickel as is shoVv11 in figure 7.4. It is 25.280Pa approximately at a
Ultimate Tensile Ultimate Tensile Ultimate Tensile
Strength values Strength values Strength values
quoted by Hertzberg obtained by obtained by MEAM
(1996) Komanduri et al Potential
(2000)
(OPa) (OPa) (GPa)
Nickel 33.4 36.0 44.46
Copper 19.1 28.0 26.20
Iron 31.8 29.0 32.16
Table 7.2 Comparison of ultimate tensile strength values obtained by using the MEAM
potential with those obtained by Komanduri et al (2000) and those quoted by Hertzberg
(1996)
strain of 0.159. The strain to fracture is 0.49.The fact that the copper has lesser strength
than nickel is clearly illustrated by these observations. The ultimate tensile strength of
iron is 32.16 OPa approximately at a strain of 0.133 and the strain to fracture is 0.455.
38
Iron is found to have lesser ultimate tensile strength than Nickel but more than Copper.
The ultimate tensile strength values of nickel, copper, and iron are compared with those
obtained from Komanduri et al (2000) and those quoted by Hertzberg (1996) in the table
7.2. The values obtained from MEAM have the same ranking as the above
7.3 SHEAR TEST OF NICKEL
Horstemeyer and Baskes (1999) conducted a shear test of a nickel specimen of
dimensions 4a x 2a x 2a at a speed of 1.0mJs and observed the shear stress vs. shear strain
curves for these samples (figure 7.8). The shear stress was normalized by the elastic shear
modulus of nickel (124.8 GPa). They observed a spatial size scale effect on yield stress.
They used the embedded atom method as the potential model. Since embedded atom
method is supported by strong physical arguments, shear of Nickel using the MEAM
potential was done to compare and validate the results with those obtained by
Horstemeyer et al.
A nickel sample with dimensions of 4a x 4a x 6a is taken and is subjected to shear
at the speed of 500 mJs (this speed is chosen to keep the computational times reasonable),
until a shear strain of 0.3. The animations are shown in figures 7.7(a)(d). The shear
stress/elastic modulus vs. the shear strain is shown in figure 7.8. It is observed that both
the curves reach a peak ofO.l(approximately) for a shear strain of 0.125 (approximately).
After that, the curve obtained by Horstemeyer (1999), drops down rapidly when
compared to the one obtained by using the MEAM potential. This is because the shear
39
test of nickel using the MEAM potential is conducted at the speed of 500 mJs and this
generates a lot of heat. This causes the material to become more ductile and hence
increases the value of strain to fracture. The maximwn shear stress/elastic modulus
values are in reasonable agreement for both the curves.
The software developed is good for single crystal materials. The
development and validation of the MD using MEAM was done along with Mr. David
Stokes. The project consisted of five important steps, developing the code for potential
function (done by me), developing the code for force function (50% by me and 50% by
David), validating the potential model (50% by me, 50% by David). conducting the
tension tests (done by me), and conducting the shear test (done by me). This software can
be used to model other manufacturing applications like cutting, milling, and indentation
by using the Morse potential as the interface potential between the tool and workpiece.
40
Figure 7.la
Figure 7.lc
Figure 7.1 b
Figure 7.ld
Figures 7.1 (a)(d) Snapshots of animation during tension test of nickel using
the MEAM potential
41
Figure 7.2 a
Figure 7.2 c
Figure 7.2 b
Figure 7.2 d
Figure 7.2(a)(d) Snapshots of animation during tension test of nickel using the
Morse potential
42
60, ,
0.8 0.8 1.2
Morse potential
0.4
t ••rrt+TMEAM potential
02
40
10
·'0' '
Strain
Figure 7.3 Tension test of nickel using the Morse potential and the
MEAM potential
40
30
10
Nickel
Copper
o.
02 04 0.8 0.8 , 2 ,. 18
·10 L '
Strain
Figure 7.4 Stress vs. strain curves of nickel, copper, and iron using
the MEAM potential
43
Figure 7.5a
Figure 7.5e
Figure 7.5b
Figure 7.5d
Figures 7.5(a)(d) Snapshots of animation during tension test of copper
44

Figure 7.6a
Figure 7.6c
Figure 7.6b
Figure 7.6d
Figures 7.6(a)(d)Snapshots of animation during tension test of iron
45
i • • •
ill. II i
II II II.
SCII
II II
II II • It •
~ ~ II • • II II W 5
S.Il •••• IlE
Ii II II II II II III tj
II "II II II II
II II
II
!II II ..
~ .1111 Ill~
Figure 7.7a Figure 7.7b
II • II
Figure 7.7c Figure 7.7d
Figures 7.7(a)(d) Snapshots of nickel in shear
46

0.150
CIl
::3 :; 0. 'CO
"'Cl
0~
.... ro
Q,) ...c 0.050
(/'J
.~ til ro
r;§ 0.000
(J)
(J)
.Q..,.)
ci5
.... ro ~.050
Q,) ...c
(/'J
Horsteme er (199~MEAM
\~r, ~
\
f  . 
0.' 0.'5
.Q.,CO L J
Shear Strain
Figure 7.8 Comparison of the shear stress/elastic shear modulus vs. shear
strain for nickel using the MEAM potential with that obtained by
Horstemeyer (1999)
47
Chapter 8
CONCLUSIONS
8.1 TESTING FOR ACCURACY OF THE SOFTWARE
The "Modified Embedded Atom Method" developed by Baskes (1992) was used
In the Molecular Dynamics Simulation. The following tests were conducted for
determining the accuracy of the software:
1. Numerical vs. analytical forces test for the validation of the forces: A system of fiveNickel
atoms was taken and one of the atoms was moved in increments of 0.1 AU and
both the analytical and numerical forces were calculated. They were in excellent
agreement up to 12 significant digits
2. Conservation of energy test: Based on the fact that an isolated system should
conserve energy, one lattice of Nickel was taken and allowed to integrate for a time
period of 10 time steps. The total energy (sum of potential and kinetic energies)
remained constant up to 12 significant digits.
3. Back Integration test: One lattice of Nickel was taken and allowed to integrate up to
100 time units and then the time step was changed to a ve value and the system was
hack integrated. It traced back the same curve for kinetic, potential and total energies
and the same path for the atom positions.
48
Based on above tests it was proved that the software developed was accurate and ready to
be used.
8.2 UNIAXIAL TENSION
Uniaxial tension experiments were done on four materials: two FCC (nickel, copper) and
a BeC (iron) material at SOOmJs using the MEAM potential. Also, uniaxial tension of
nickel, using the Morse potential, was performed and compared with the resuJts obtained
from the MEAM potential. The foHowing observations were made:
1. The measured ultimate tensile strengths were 26.20 OPa (copper), 44.46 OPa
(nickel), and 32.16 GPa (iron). These values have the same ranking as the theoretical
\'alues of defect free whiskers given in Hertzberg (] 998) and also as those obtained
by Komanduri et al (2000).
2. Copper has the least strength (26.20 GPa) while nickel has the highest strength (44.46
OPa) among the three materials. In bulk also, copper has the least strength of the
three materials and nickel has the highest.
3. The strain to fracture at room temperature for Cu, Ni, and Fe, 0.493, 0.43. and 0.461
respectively. Copper exhibited the maximum ductility undergoing a maximum strain
0[0.493 before failure.
By the above observations it is clear that the MEAM potentials represent both the FCC
and BCC materials fairly well and can be used to perform tensile testing of these
materials,
49
8.3 SHEAR TEST OF NICKEL
Shear test of Nickel was conducted at a speed of SOOm/s and the shear stress, which
was nonnalized by dividing it with the elastic shear modulus of Nickel (124.8 GPa) vs.
shear strain was plotted. The curve reached a maximum of 0.1 at a shear strain of
~O.125. The value of maximum shear stress/ elastic modulus was in reasonable
agreement with that obtained by Baskes and Horstemeyer (1999). After reaching the
peak, the curve obtained by MEAM dropped less rapidly than that obtained by
Horstemeyer. The difference is attributed to the fact that the shear test conducted by
using the MEAM potential was conducted at a high speed, and this increased the
temperature of the specimen. This lead to higher ductility and increased value of strain to
fracture.
8.4 FUTURE WORK
I. The existence of stray atoms might be because of the small specimen sizes. Due to
time constraint on this project, the experiments were conducted with smaller sizes.
Experiments must be conducted using large sizes to find out the sizescale effects on
the ultimate tensile strength and strain to fracture.
2. The model developed above is good for single crystal materials. It can be used to
define the potential between like atoms. An alternate potential model like the Morse
potential can be used to define the potential between unlike atoms. Then this model
can be extended to other manufacturing applications like cutting, milling, and
indentation, where the Morse potential can be used to define the interactions between
tool and workpiece (unlike atoms).
50
3. Crystals with defects and cracks can be modeled using the MEAM potential. This can
be done by introducing a defect by removing few atoms randomly from the
workpiece to simulate a defect, or by removing few atoms from the workpiece in the
shape of the crack to be studied to simulate a crack. Uniaxial tension and shear
experiments should then be conducted on materials to study their effects on the
tensile properties.
51
REFERENCES
Adams, lB., Foiles, S.M, "Development of an Embedded Atom Potential for a BCC
Metal: Vanadium," Physical Review B, 41/6, pp.33163327, 1990.
Allen, M., and Tildesley, D., "Computer Simulation of Liquids," Oxford University
Press, Oxford, U.K, 1991 .
• Alder, B.l., Wainwright, T.E., "Studies in Molecular Dynamics. I. General Method,"
The Journal of Chemical Physics, 31/2, pp.459466, 1959.
Baskes, M.L, Daw, M.S., "Semiempirical, Quantum Mechanical Calculation of
Hydrogen Embrittlementin Metals:' Physical Review Letters, 50117, pp.1285, 1983.
Baskes, M.l., Daw. M.S., "Embedded Atom Method: Derivation and Application to
Impurities, Surfaces, and other Defects in Metals," Physical Review B, 29/12, pp.64436453,
1984.
Baskes, M.L, Foiles, S.M., Daw, M.S., uEmhedded Atom Methodfunctionsfor Cu, Ag,
Au, Ni, Pd, Pt, and their alloys," Physical Review B, 33/12, pp.79837991, 1986.
Baskes, M.I., "Application of the EmbeddedAtom Method to Covalent Materials: A
Semiempirical Potentialfor Silicon," Physical Review Letters, 59, pp.2666, 1987.
Baskes, M.I., Foiles, M.F., Daw, M.S., "Atomistic Studies of Interfacial Structure and
Properties," Materials Research Soc. Symp. Proc., 122, pp.343353, 1988.
Baskes, M.I., Nelson, 1.S., Wright, A.F., "Semiempirical Modified Embedded Atom
Potentials for Silicon and Germanium," Physical Review B, 40/9, pp.60856100, 1989.
Baskes, M.I., "Modified Embedded Atom Potentials for Cubic Materials and
Impurities," Physical Review B, 46/5, pp.2727, 1992.
52
Baskes, M.l.,"MEAM for HCP Metals," Modeling Simulation Mater. Sci. Eng., 2,
pp.l47, 1994.
Baskes, M.I., "Determination of Modified Embedded Atom Method Parameters for
Nickel," Materials Chemistry and Physics, 50, pp.152158, 1997.
Belak, J., Boercker, D.E., Stowers, I.F., "Simulation of NanometerScale Deformation
ofMetallic and Ceramic Surfaces," MRS Bulletin, pp.5560, 1993.
~ Chandrasekaran, N., "Length Restricted Molecular Dynamics Simulation of
Nanometric Cutting," Thesis, Oklahoma State University, 1997.
Daw, M.l., "Model of Metallic Cohesion: The Embedded Atom Method," Physical
Review B, 39/11, pp. 7441745 L 1989.
p Eyring, H., Walter, 1., Kimball, G.E., ., Quantum Chemistry," John Wiley & Sons, Inc.,
1944.
o Haile, J., "Molecular Dynamics Simulation Elementary Methods," John Wiley & Sons,
Inc, 1992.
Heino, P., Hakkinen, H., Kaski, K., "Molecular Dynamics Study of Mechanical
Properties of Copper," Europhysics Letters, 41/3, pp.273278, 1998.
Hertzberg, R.W., "Deformation and Fracture Mechanics ofEngineering Materials," 4th
Edition, John Wiley & Sons, NY, 1996.
Horstemeyer, M.F., Baskes, M.I, "Atomistic finite deformation simulations: A
discussion on length scale effects in relation to mechanical stresses," Journal of
Engineering Materials and Technology, 121/2, pp.11411 9, 1999.
53
Kitamura, T., Yashiro, K., Ohtani, R., "Atomistic Simulation on Deformation and
Fracture ofNanoSingle Crystal ofNickel in Tension," JSME Int. A, 40/4, pp.430435,
1997.
Komanduri, R., Chandrasekaran, N., Raff, L.M, {'Molecular Dynamics Simulation of
Uniaxial Tension," submitted for publication.
LyndenBell, R.M., "Computer Simulation of Fracture at Atomistic Level," 1. Phy.
Condens. Matter, 4, pp.17041705, 1994.
LyndenBell, R.M., " A Simulation Study ofInduced Disorder, Failure and Fracture of
Perfect Metal Crystals under Uniaxial Tension," 1. Phy. Condens. Matter, 7, pp.46034624,1995.
Raff, L.M., and Thompson, D.A., "The Classical Trajectory Approach to Reactive
Scattering," Chapter 1, Volume 3 in Theory of Chemical Reaction Dynamics, Ed.
Michael Baer, CRC Press, 2121, 1986.
. Rappaport, D.C., "The art of Molecular Dynamics Simulation," Cambridge University
Press, 1995.
Rentsch, R., Inaski, I., "Indentation Simulation on Brittle Materials by Molecular
Dynamics," Annals ofCIRP, 44/1, pp.295, 1995.
Sato, K., Yoshioka, T., Ando, T., Shikida, M., Kawabata, T., "Tensile Testing ofSilicon
Film Having Different CrystaLlographic Orientations Carried out on a Silicon Chip,"
Sensors and Actuators, A 70, pp.148152, 1998.
Stewart, R., {{Investigation on Molecular Dynamics Simulation of Nanometric
Cutting," Thesis, Oklahoma State University, 1998.
54
Stott, M.1., Zaremba, E.. "Quasiatoms: An Approach 10 Atoms in Nonuniform
electronic systems," Physical Review B, 22/4, pp.15641583, 1982.
Voter, A.F., Chen, S.P., Mater. Res. Soc. Symp. Proc., 82, pp.17S, 1987.
Vi, T., Li, 1., Kim, C.J., "Microscale Material Testing of Single Crystal Silicon:
Process Effects on Surface Morphology and Tensile Strength," Sensors and Actuators,
83, pp.172178, 2000.
55
, J
\ VITA
Kalyan K. Mavuleti
Candidate for the Degree of
Master of Science
Thesis: MOLECULAR DYNAMICS SIMULAnON USING THE MODIFIED
EMBEDDED ATOM METHOD
Major Field: Mechanical Engineering
Biographical:
Personal Data: Born in Gudivada, India, on January 3, 1976, the son of Ganga
Raju Mavuleti and V.V. Subramanyeswari Mavuleti
Education: Received Bachelor of Technology degree in Mechanical Engineering
from Jawaharlal Nehru Technological University, Kakinada in June 1998.
Completed the requirements for the Master of Science degree with a major
in Mechanical Engineering at Oklahoma State University in Oct, 2000.
Experience: Research Assistant in Oklahoma State University from August'98December
2000.