Characterization of basic material properties at microscopy level has been recognized to be more and more important lately. In despite of the development and application of electron microscopy on the microstructural analysis, observation of the underlying energetic and interaction mechanism at the atomic level is not always accessible by experimental tools in most situations. In such cases, using the computer-based techniques seems to be promising [140]. Molecular dynamics (MD) and Monte Carlo (MC) simulations are two of the powerful atomic level simulate methods that are widely applied to study the atomic topology of BMGs in the recent years [30-33, 36, 92-94]. They can not only obtain strict numerical results about a complex system in microscopic scale, but can also produce macroscopic behavior under precisely controlled conditions. As compared with the MC method, time dependence is the advantage of MD to simulate the physical mechanism of rate relationship. For example, it is possible to consider the dependence of the glass properties on quench rates[141], or mechanical properties on strain rates [126].
A brief description of the MD methods is given in the following sections; the detailed presentations can be referred to textbooks elsewhere [140, 142].
3-1 Equations of motion and potential function
In the MD method, the classical Newton’s motion equations are adopted to solve for the atoms molecules as
i i
i
i F r
dt r m d
∂ Φ
−∂
=
2 =
2
, (3-1)
24
where mi, ri, Fi are mass, position vector, force vector of molecule i, respectively, and Φ is the potential of the system. The potential function Φ describes how the potential energy of a system of N atoms depends on the coordinates of the atoms, r1, r2, …, rN. It is assumed that the electrons adjust to new atomic positions much faster than the motion of the atomic nuclei (Born-Oppenheimer approximation). When the discrete time steps Δt which is usually given between 1-10 fs and a set of thermodynamic variables such as volume, pressure and temperature are chosen, the time evolution of the system is determined through a deterministic calculation of the trajectory of each particle in the phase space. This is in contrast to the MC method [140, 142], where the system evolves according to stochastic dynamics by setting up a random walk to sample configurations via a given distribution function.
It is well known that the classic equation of motion is a good approximation of the Schrödinger equation if the mass of atom is not too small and the system temperature is not too low. Once the potential of the system is taken, Eq. (3-1) should be solved by numerical methods. In principal, any of the gas, liquid, solid states, and inter-phase phenomena can be solved without the knowledge of "thermo-physical properties" such as thermal conductivity, viscosity, latent heat, saturation temperature and surface tension [143, 144].
The potential of a system Φ (r1, r2,…..rN ) can often be reasonably assumed to be the sum of the effective pair potential φ (rij) as
∑∑
>= Φ
i j i
r )ij
φ( , (3-2)
where rij is the distance between molecules i and j. The assumption of Eq. (3-2) is often employed for simplicity even though the validity is questionable. The covalent system such as carbon and silicon cannot accept the pair-potential approximation. Such categories of
25
potential is developed from the Schrödinger equation that accounts for the position of nuclei and electrons, but using computations such as first principle theory describing electronic degrees of freedom are severely computation restrictive [143, 144]. About 1,000 atoms are the utmost amount that can be solved by state-of-the-art computations. To increase a calculating length scale, various approximation strategies have been developed, in which the electronic effects are averaged out and only the nuclear coordinates are considered as atomic coordinates for instance pair potentials. Lennard-Jones potential is a typical example of pair potential which has form given below [143, 144]:
⎥⎥
⎦
⎤
⎢⎢
⎣
⎡ ⎟
⎠
⎜ ⎞
⎝
−⎛
⎟⎠
⎜ ⎞
⎝
= ⎛
6 12
4 )
(r ε σr σr
φ , (3-3)
where ε and σ are energy and length scales, respectively, and r is the intermolecular distance, as shown in Fig. 3-1. This potential has a well description on the characteristics of inert monatomic molecules (the parameters of ε and σ for each molecule are listed in Table 3-1). It also has a well investigation on the atomic behavior, which the intermolecular van der Waals interaction is the major influence in the system, of about several million-billion atoms.
However, the nature of bonding is not accounted for in pair potentials. This leads to unphysical effects in the core of lattice defects such as unusually high stacking fault energy in metals [145]. Furthermore, pair potentials can not account for environmental dependence due to the many body effects in the intermolecular interaction, e.g. atoms in the bulk is too similar to the atoms on the surface or near a defect sites, which are important in metals.
For metallic behaviors, an alternative simple but rather realistic approach to the description of bonding in metallic systems is based on the concept of local density. This allows one to account for the dependence of the strength of individual bonds on the local environment which is especially important for simulation of surfaces and defects. Many
26
methods have been proposed since early 1980s, e.g. embedded-atom method - EAM, effective medium theory - EMT, Finnis-Sinclair potential, the glue model, corrected effective medium potential - CEM, etc. and are based on different physical arguments (e.g.
tight-binding model, effective-medium theory), but lead to a similar expression for the total energy of the system of N atoms as:
∑
=
i i
total E
E , (3-4)
∑
≠+
=
i j
ij i
i
i F r
E (ρ) 21 φ( ), (3-5)
and
∑
≠=
i j
ij i
i f(r )
ρ .
(3-6)
Interpretation and the functional forms of F, f, and φ depend on a particular method.
Take the point of view of effective medium theory or the embedded-atom method for example, the sketch as shown in Fig. 3-2, the energy of the atom i is determined by the local electron density at the position of the atom and the function f describes the contribution to the electronic density at the site of the atom i from all atoms j. The sum over function f is therefore a measure of local electron density ρi. The embedding energy F is the energy associated with an atom embeded in the electron environment described by ρ. The pair-potential term φ describes their electrostatic contributions. The general form of the potential can be considered as a generalization of the basic idea of the density functional theory that the local electron density can be used to calculate the energy. Complying with different physical interpretations, the methods to determine the way function are different as well. Some functions and parameters are driven from the first principles calculations and others are guessed or fitted to the experimental data. Generally, the results are basically rather similar.
27
Good ability to describe the variation of the bond strength with coordination is the main advantage of these methods compared with pair potentials. Increase of coordination would decrease the strength of each of the individual bonds and increase the bond length, resulting in a reproduction for the bulk, surface, and defects properties. In the past decades, such category potentials have been widely used in several studies involving modeling of grain boundaries, fracture studies, study of surfaces, wear and friction.
3-2 Ensembles
Once the potential is chosen, the trajectories of atoms are determined straightforward by employing Newton’s laws to the force fields obtained from radial derivatives of potentials in molecular dynamics [140]. The classical molecular dynamics simulations are performed on a NVE (N is the number of particles, V is the volume, E is the total potential energy) ensemble, which considers the movement of a constant number of particles in a box of fixed size and shape, and the system is assumed to be free from any external force so that number of particles, volume, and the total potential energy are conserved, that is, a micro-canonical ensemble. The advantage of NVE ensemble is easy to handle but difficult to make direct comparison with experiments because laboratory experiments are usually carried out at constant pressure and temperature.
Additional canonical ensemble commonly applied to include the thermodynamic effect in MD is the isothermal-isobaric (NPT) ensemble which further considers the fluctuation of temperature (T) and pressure (P) conditions [140]. The total energy in this ensemble is fluctuated by the interaction with a piston and through the thermal contact with a heat bath.
The temperature T is related classically to the kinetic energy, and therefore the isothermal
28
condition in the ensemble is realized by adjusting the velocities vi of the particles in every time step of the simulation [140]:
∑
=
→ N
i i i B
i
i v Nk T mv
v
1
/ 2
3 . (3-7)
This velocity-scaling method is an approximation by using Gaussian constraint method to fit the correct canonical distribution in the coordinate space [146, 147]. Similarly, a simple relation of pressure in MD is expressed by the virial equation:
ex
N N
i j
ij ij B ij
P r r
r V
T
P Nk ⎟⎟+
⎠
⎞
⎜⎜
⎝
⎛
∂ Φ
− ∂
=
∑ ∑
>
1
) ( 3
1
3 , (3-8)
where NkBT is related to the average of the kinetic energy of the system as
ex N
i
N
i j
ij ij
ij i
i r p
r v r
V m
P ⎟⎟+
⎠
⎞
⎜⎜
⎝
⎛
∂ Φ
− ∂
=
∑ ∑
=1 >
2 ( )
3
1 , (3-9)
where vi is the velocity of particle i, and rij is the vector joining particle i to particle j.
ij ij
ij r r
F =−∂Φ( )/∂ is the force exerted by particle j on particle i and Pex is the external pressure acting on system. From this representation, the pressure can be adjusted by expanding or contracting the simulation box.
In the history of ensemble development, Andersen [148] first advanced a refined constant pressure ensemble by the extended Lagrangian method. In his method, the volume is viewed as a degree of freedom which fluctuates with the consistent kinetic and potential term in the total Hamiltonian of an extended system. However, the constant temperature is performed by stochastic collisions with a heat bath in the Andersen’s approximation. Nosé and Hoover [149-151] extended in the similar method based on a clever use of an extended
29
Lagrangian to treat the constant temperature condition by adding a degree of freedom which describes the coupling to the heat bath. Subsequently, Parrinello and Rahman [152] extended this technique to allow a change of both the shape and volume of the simulation box. The Parrinello and Rahman scheme made the ensemble of MD application be complete especially for studying phase transitions in solids.
3-3 Integration of the Newtonian equation
The most common MD algorithm for solving the integration of the equation of motion is the Verlet method [140], the integration scheme of Verlet can be simply derived by the Taylor expansion of the equation of motion as follows:
i i i
i
i m
t t F t
tv t r t t
r ( )
2 ) 1 ( )
( )
( +Δ = +Δ + Δ 2 , (3-10)
i i i
i
i m
t t F t
tv t r t t
r ( )
2 ) 1 ( )
( )
( −Δ = −Δ + Δ 2 . (3-11)
Through summing up and rearranging, it becomes
i i i
i
i m
t t F t t r t r t t
r ( )
) ( ) ( 2 )
( +Δ = − −Δ +Δ 2 . (3-12)
The velocities can be calculated from:
{
r t t r t t}
tt
vi( )= i( +Δ )− i( −Δ ) /2Δ . (3-13)
Another common algorithm is the leap-frog method as shown below:
i i i
i m
t tF t t
t v t
v ( )
2
2 ⎟+Δ
⎠
⎜ ⎞
⎝
⎛ −Δ
⎟=
⎠
⎜ ⎞
⎝
⎛ +Δ , (3-14)
( ) ( )
⎟⎠
⎜ ⎞
⎝
⎛ + Δ Δ
+
= Δ
+ 2
t t tv t r t t
ri i i . (3-15)
After the velocity of each particle is calculated by Eq. (3-14), the position is obtained by
30
Eq. (3-15). More elaborate integration schemes such as Gear’s predictor-corrector method and Beeman algorithms [140] are sometimes employed depending on the complexity of the potential function and the demand of the accuracy of motion in each time step.
3-4 Periodic boundary conditions
Due to the restriction of computations, it is difficult to deal with the number of atoms as practical condition in the MD system. The numerical system with a typical particle number in the range of 103 is necessarily limited in size. A system under this size should be considerably affected by the free surface effects that easily lead to a deviation on the simulated results as compared with bulk properties. Periodic boundary conditions are usually used to minimize these finite-size effects [140, 143, 144]. This approach consists in a periodic repetition of the simulation box in the three directions to fill the whole space, as illustrated in Fig. 3-3 [143].
Any interaction of molecules beyond the periodic boundary is calculated with replica molecules. This utilization removes all free surfaces. By construction, each molecule possesses infinity of periodic images.
In order to avoid the calculation of potential between a molecule and its own replica, the minimal image convention is adopted [140, 143]. Only interaction with the closest images is taken into account in this scheme. The range of the interaction is assumed to be smaller than half of the cubic box length.
3-5 List method and cut-off radius
Since most potentials for modeling metals are short-ranged force, the interactive strength of two atoms would be reduced rapidly as increasing the distance between them and would
31
decay smoothly to zero at a distance rc. That is the cut-off radius of the potential in the MD simulations. The advantage of using cut-off radius is to reduce the laborious forces calculations due to the large number of particles pairs in a system, thus the interactions are limited to the neighbors at distances smaller than rc. The calculation of the different N(N
−1)/2 distances in every integration step also wastes a significant computing time. This can be considerably reduced by making use of the some techniques such as Verlet-list and cell-list [143, 153]. To every atom, one applies a neighbor list which contains all particles within a Verlet-radius rv, chosen to be reasonably larger than rc. Only the particles belonging to its list are considered in calculating the forces acting on a given atom. The neighbor list is updated once every about 5 to 20 integration steps depending on the mobility of the particles, such that the estimated maximal displacement of the particles between two updating steps remains smaller than the difference rv − rc,as shown in Fig. 3-4 [143, 153].
Figure 3-5 [143, 153] perform another scheme, cell-list, which is often used to treat a system with larger particle numbers. In this technique, the cubic simulation box is divided into a regular lattice of ncell × ncell × ncell cells, and the side of the cell lcell = L/ncell is greater than the potential cutoff distance rcut. Thus, it is only necessary to look for atoms in the same cell as the atom of interest, and in its nearest neighbour cells. Currently, it becomes preferable to apply the Cell-list linked Verlet-list in order to enhance the efficiency for large systems with short-range forces [140, 143, 153]. A flow chart of a typical molecular dynamics code is illustrated in Fig. 3-6.
32