# Basics of Molecular Dynamics Simulation

## Introduction

Understanding and analyzing the atomic-level structural evolution and atomic interactions during the deformation processes are quite challenging. In this context, the computer simulation methods can be most suitable for visualizing and investigating the deformation mechanisms, thereby predicting the bulk material properties accurately [1]. In the current scenario, the available simulation methods are broadly covered under two classes: molecular dynamics (MD) simulations and Monte Carlo (MC) simulations. The MD simulation technique is preferred over the MC simulation technique for the evaluation of the dynamic property of many-body systems and its accuracy [1,2]. In this chapter, we specifically deal with the governing laws, different algorithms to set up a simulation, atomic interactions, and the post-processing of the results obtained from MD simulations.

“Molecular dynamics,” in general, implies a solution for the Newton’s law of motion. It is applied over all particles present in the system to describe their displacement and interactions numerically over time to interpret and predict the possible changes in the system [3,4]. Alder and Wainwright were the first to report a work using MD simulation in 1957 that focused on the study of phase transition behavior in a hard-sphere system [5]. In 1960, MD simulation was implemented for real material model to study the radiation damage in Cu by Gibson and another group led by Vineyard [6], whereas, in 1964, Rehman simulated a model of the atomic motion in real liquid (Argon), using MD simulations [7]. Later on, in 1967, Velocity-Verlet algorithm was developed to solve the classical equations of motion [8]; it was implemented by Rahman et al. [9] and McCammon et al. [10] to study more complex systems such as liquid water and protein. With advancement in the computing systems, multi-tasking, multi-threading, and multi-level supercomputers that controlled the massive parallel simulations and henceforth escalated the process efficiency and performance by significantly reducing the computational time were developed [11,12]. This led to the rapid development of MD simulations, which are now the most widely used computational tools to evaluate different properties that comprise the mechanical [13-17], structural [17-20], thermodynamic [17.18,20-22], physical [19], surface [23,24], vibrational [16,17], thermo-physical [25], elastic [26,27], and electronic [16,27,28] stretching to various domains of research such as material science, biology, chemistry, and physics. Specifically in the last two decades, MD simulation has gained a lot of popularity among researchers in the field of nanomaterials. Wide-range investigations have been carried out to analyze the bulk material properties and deformation behavior of nanowires [29,30], nanopillars [31,32], nanorods [33,34], nanoporous materials [35,36], nanofoams [37], nanotubes [14,16,17], nanofibers [18,38], nanoparticles [39,40], nanocomposites [26, 40], nanofilms [41,42], nanobubbles [43,44], nanobeams [45], and nanocrystalline materials [46,47]. Some of these MD-simulated studies have a great impact on modern research and scientific society. For instance, Yin et al. have studied the deformation of Cu/Al multilayers under compression [48], whereas Song et al. have analyzed the mechanical behavior of amorphous/amorphous nano-multilayers to develop and design high- strength metallic glasses [49]. Similar deformation studies using MD simulations include study of the mechanical behavior of hard amorphous Si shell around the Au core nanowire under cyclic tensile loading and compressive unloading [50], investigation of the deformation behavior of high-entropy alloy-coated metals (Ni and Cu) during nanoindentation [51,52], functioning of nanopores as a switchable gate to allow the passage of nanoparticles [53], and analysis of the mechanical behavior and morphological changes during amorphization of heat-treated nanoporous Au [54]. However, MD simulation has a number of challenges in the quantum effect studies, formulation of appropriate force field, studies in large-sized atomic systems, and long-duration processes [55-57]. Despite all the challenges, there has been progress in this research domain too, as many research groups and institutes are focused on developing accurate potentials and improving the efficiency of supercomputers in handling large systems. Recently, Robustelli et al. have successfully developed an MD force field 99SB*-disp* that is in agreement with the experimental results in the study of disordered and folded protein states [58]. On the other hand, Abraham et al. [59] have carried out large-scale atomic simulation using 12 TeraFlop Accelerated Strategic Computing Initiative (ASCI) white computer systems to investigate the underlying deformation mechanism behind the ductile failure of material using a large specimen having one billion atoms. Overall, MD simulation has no alternatives that can be best suited for the atomic- scale studies related to phase transformation, diffusional phenomena, deformation process, nucleation-growth process, solvation process, defect evolution, structural changes, crystalline-amorphous transition, creep phenomenon, ratcheting behavior, and molecular interactions.

## Molecular Interactions

Molecular interactions can be described well in terms of correct force field selection, which is a collection of equations and their respective force constants to construct the molecular geometry and their required properties. The interatomic force, *f _{h}* acting on each atom in the system can be calculated by taking the negative gradient of its total potential energy. In turn, the error-free potential energy calculation is highly capable of predicting molecular interactions in an efficient way. The total potential energy comprises non-bonded and bonded interactions and can be expressed as follows:

### Bonded Interactions

The bonded potential mainly consists of interactions of two-, three-, and four-body covalently bonded atoms. This further contributes to the potential due to all types of bonds, angles, and torsional angles, respectively. The two-body bond potential between the two covalently bonded (/,./) pairs of atoms can be expressed in terms of spring harmonic vibrational motion as:

where *щ* is the distance between two adjacent pairs of atoms and is equal to |r_{f} -i)|, r_{0} is the corresponding equilibrium distance, and *K _{h}* is the bond force constant. The bonding potential between three-body bonded atoms (

*i,j,k*) experiences oscillation about an equilibrium angle due to angular vibrational motion and can be expressed as:

where *в,* is the angle between two position vectors such as ^ and *r _{Jk}, в„* is the equilibrium angle, and

*K*is the angle force constant. The potential between four-body bonded atoms (

_{e}*i,j,k,l*) experiences a torsional rotation about the central bond and is given by

where y, is the phase shift, *n* is the dihedral multiplicity,

K_{0} is the dihedral force constant. Now, taking Eqs. 5.2 through 5.4, the net bonded potential can be expressed as:

### Non-bonded Interactions

The most commonly used non-bonded potential is Lennard-Jones (LJ) potential, which is analogous to van der Waals force of interaction. In 1924, John Lennard- Jones [60] proposed LJ interatomic potential, which most appropriately represented the interactions between a pair of inert gases or neutral atoms or molecules. In the past, LJ potential has been successfully used to study the properties of liquid argon (Ar) [7], which can be expressed as:

where cr„ is the diameter or distance at which the interparticle potential becomes zero, e,_{y} is the wall depth, and is the distance between the particles. The intermo- lecular distance (r) between the particles plays a very significant role in determining the intermolecular force. The repulsive term ^{12} in Eq. 5.6 elucidates the short-range overlapping of electron orbitals, whereas attractive term *r~ ^{6}* explains the long-range intermolecular forces. In short, it can be said that the particles repel each other at small values of and attract at higher values of

*>у.*The second most non-bonded interactions are the electrostatic interactions, which are expressed by columbic forces of interaction as:

where *q _{t}* and

*q,*are the charges,

*ц*is the distance between two charges, and e

_{0}is the permittivity of free space or vacuum permittivity. Hence, the net potential toward non-bonded interaction is given by

The determination of the correct force field is based on the validation with experimental results. If the selected force field reproduces almost the same result as the experimental results, then it can be used to characterize the bulk material properties precisely. Commonly used force fields such as Assisted Model Building with Energy Refinement (AMBER) [61,62], Groningen Molecular Simulation (GROMOS) [63], Chemistry at HARvard Macromolecular Mechanics (CHARMM) [64], and Optimized Potential for Liquid Simulation (OPLS) [65] are often used to predict the bulk material structures and properties accurately.

## Interatomic Potentials

Interatomic potentials are preferred for describing the interaction between atoms in the system. Such description is significant for developing the theoretical models of the system. These models are used for predicting the characteristics of the system and describing the system behavior in different circumstances [66]. Such interactions between the constituent atoms are described by using interatomic potentials. In order to define or construct an interatomic potential, the potential energy of the system is split into several terms based on the coordinates of individual atoms, pairs, triplets, etc, as given in Eq. 5.9 [67]:

The first term associated with the equation represents the external potential, that is, the effect of an external field. The second term representing the pair potential is based on the magnitude of pair separation *r^* given by the equation

The third term in Eq. 5.9 is the potential involving three particles and is rarely considered for computer simulations, owing to large time consumption in computing the quantities of the triplets of molecules. The variation of pair potential with the distance between the two atoms is presented in Figure 5.1.

As is evident from the figure, the potential curve shows an attractive tail at large distance values, thereby exhibiting typical features of intermolecular interactions. A negative well in the figure represents the attractive forces between the atoms due to Coulomb interactions, covalent bonding, metallic bonding, or van der Waals interactions. However, as the distance further decreases, a steep rise in the curve can be observed, signifying the repulsive forces between the atoms attributed to Columbic repulsion between the nuclei [67]. As mentioned earlier, the three-particle

potential term is excluded in the potential calculation, but the average of many- particle effect can be included by defining a term “effective” pair potential and modifying Eq. 5.9 to

Although pair potentials are extremely simple, they are widely employed in determining the behavior of many systems. Some widely used interatomic potentials are given in the following sub-sections.

### Lennard-Jones Potential

One of the first potentials describing the van der Waals interaction employed for simulation of real system, for instance, condensed noble gases (i.e., in liquid state), has been suggested by Lennard-Jones [60]. The attractive forces between the atoms are attributed to van der Waals dispersion interactions. However, the repulsive forces are attributed to the Pauli exclusion principle, which states that “two electrons belonging to a particular atom cannot have same quantum number,” thereby adding up the terms results in Eq. 5.6, which describes the functional form of the LJ potential. The depth of the potential well represented by *e* parameter is computed from the binding energy. On the other hand, the equilibrium lattice parameter determines the *a* parameter of the equation. Occasionally, the LJ potential is also used to determine the behavior of metals. This is attributed to the simplicity of the LJ potential, but failure in determining the accurate behavior occurs at very close interatomic distance, as the LJ potential is proportional to Ar LJ potential curve with varying interatomic distance is presented in Figure 5.2.

### Morse Potential

Morse suggested the potential energy curve based on the chemical bond between two atoms, given by the equation

where *D _{c}* represents the depth of the well in the potential curve,

*r*stands for the minimum energy distance when potential becomes zero, and

_{0}*a*is given by the equation

FIGURE 5.2 LJ potential of liquid Ar as a function of separation distance between Cu atoms.

The Morse potential provides analytical solution for Schrodinger’s equation and poses lesser problems with the cutoff distance, as compared with the LJ potential, on the account of the fact that the curve decays faster than the LJ potential curve (refer to Figure 5.1). As can be seen from Eq. 5.13, the Morse potential can be described based on three parameters, namely binding energy, lattice parameter, and bulk modulus. Therefore, the potential can be employed to describe the behavior of different crystal structures such as face-centered cubic, body-centered cubic, and hexagonal close packed [68,69].

All potentials discussed earlier are pair potentials, where the force calculations are simplified based on the fact that Т,_{(} = *V _{/r}* since the potential is a function of the distance between the two atoms (i.e., /;,■ =|/;

_{;}| = |r

_{/(}-|). But in many-body potentials,

*Vjj *Vjj,*thereby making the calculations of potential more complicated.

### Embedded-Atom Method

Embedded-atom method (EAM), based on the density functional theory (DFT), was proposed first by Daw and Baskes [70]. It has been used for simulating interatomic interaction of materials. The DFT helps to calculate the interatomic forces based on the structure of the material. It is determined by an electric field formed by the distribution of nuclei and electron charges.

Based on this principle, Hohenberg and Kohn proposed that the state of an electron in a multielectronic system can be determined by means of electron density, which is a function of coordinates. In EAM, the energy of an atom in the system is the resultant sum of the energy of the pair interaction of the atom under consideration with other atoms in the system, as well the embedded function.

The embedded function considers the effect of the surroundings. In order to determine the EAM potential, the following three factors are required to be evaluated:

- • The energy of the pair interaction (
- • Embedded function (F)
- • Electron density (p)

The three factors of the EAM potential are given by the following equations:

The terms *cp*_{e}, *p, y,* and *r _{c}* are determined by fitting to the characteristics of the material;

*E*is the energy associated with the formation of a vacancy

_{v}

where *В* denotes the bulk modulus, *V _{a}* indicates atomic volume, and £

_{0}is the binding energy.

Although EAM successfully predicts the properties of metals and alloys, it fails miserably to determine the properties of structures with covalent bonds. The EAM potentials are widely used to determine the atomic behavior of metallic specimens [71,72]. The EAM potentials developed for various systems (atoms and compounds) can be easily found in the databases [73].

### Modified Embedded-Atom Method

Similar to EAM, modified EAM (MEAM) also comprises three functions, namely pair interaction, embedded function, and electron density but with modified definition. For instance, the electron density in EAM is spherical and uniform for all constituent atoms. However, each component defining the electron density in case of MEAM consists of various directions of the bond or takes into consideration the directionality of the bond. Therefore, unlike EAM, MEAM can be employed to compute the atomic interaction of a specimen with different types of bonds. The electron density in case of MEAM can be mathematically represented as

where *pf ^{h)}* denotes the atom electron densities from

*j*atom at

*R*distance from

_{4}*i*site and

*R°i*is the vector component of distance vector between two atoms

*i*and

*j*(where

*a*can be

*x, y,*or

*z*component of the vector). Modified EAM is employed to determine the physical and mechanical behaviors of metals, alloys, oxides, and semiconductors, owing to its ability to discern between the different types of bonds [74-77].

The modified embedded function employed in MEAM is given as follows:

where *A* is a constant, *E _{c}* is the binding energy, p is the electron density, and p

_{0}is the electron density corresponding to a standard reference atomic structure.

Finally, the pair potential term is modified by adding a screening function to the term, taking into account the influence of the first-nearest neighboring atom or second-nearest neighboring atom on the pair-interaction function, which are termed 1NN-MEAM and 2NN-MEAM, respectively (Figure 5.3).

Atom *X* is inside the ellipse, the outer bound of which is defined by C_{min}; atom *X *will screen the interaction between atom *A* and atom *В* completely. If atom *X* is outside the ellipse, the outer bound of which is defined by C_{max}, there will be no screening of interaction between atom *A* and atom *В* via atom *X.* The screening effect will gradually vary when atom *X's* position is in between C_{mi}„ and C_{max}.

### Charge-Optimized Many-Body Potentials

The ability of the potential to model dissimilar materials and appreciable transfer- ability of the potentials are some of the advantageous features of charge-optimized many-body (COMB) potentials. Transferability of a potential may be defined as its

FIGURE 5.3 Screening function mechanism.

ability to compute accurate or near-accurate results despite the alteration in the simulating conditions from those used in fitting procedure. The COMB potential is bond- order and dynamic-charge-based potential. It can be considered as a sum of energies, given in Eq. 5.19 [78].

The term £^{self} stands for self-energy, constituting the ionization energy of an isolated atom and the variation in the ionization energy of the ion when it is embedded in a crystal unit or a molecule. The field effect term (indicated by £^{tl,ul}) includes electron-electron interaction and nuclear-electron charge interaction between two atoms. Furthermore, the energy due to self-dipole, charge-dipole, and dipole-dipole interactions is contributed by the molecular polarization term (£^{polar}). The nonbonding interactions are included via van der Waals energy term (£^{vdW}).

The bond-order term denoted by £^{bond} constitutes the overall effect of pair-wise repulsive and attractive energy, which decays exponentially with increasing interatomic distance [79]. In addition, the energy associated with bond angle, torsion, conjugation, and coordination is also included in the bond-order term. Finally, the £^{0,hers} expressed in Eq. 5.19 can be considered a correction term used in the equation associated with bond angle.