Molecular Dynamics Simulations
Molecular Dynamics simulation (MD) is a computational method which calculates the trajectories of objects such as atoms and molcules
in phase space. The calculations are performed according to Newton's equations of motion for N interacting particles. The interaction energy is
calculated according to atom position only, using an effective interatomic potential which, in turn, is calculated using various methods and approximations.
In most molecular dynamics simulations the potential is given as an input file and isn't calculated as part of the simulation.
The MD simulation method uses theoretical elements from the field of statistical mechanics in order to simulate different aspects of material behaviour.
It is focused, usually, on the smallest of elements - atomic or molecular movement.
In order to create a simulation, some assumptions must be made. First, we treat all of the atoms as though they are an ensumble of interacting point particles. Each point particle has its own position and velocity vector in space. The integration is performed by solving the system of N equations of motion for each time step until a certain criteria are met, such as number of runs, minimal change and so forth.
In this project, molecular dynamics is used for building and simulating the motion of atoms from the crystal lattice of Si3N4 or from the amorphous phase.
Interatomic potential: The Tersoff potential
The interatomic potential describes the interactions between atoms. The potential is usually a function of the position of the nuclei which represent the potential energy of the system. Mostly, and also in our case, only two-body and three-body interactions are taken into consideration. A cutoff radius is usually used so the computer has to calculat only the interaction of each atom with the atoms within its cutoff radius. This enable us to neglect the interactions between atoms which are far from each other, and therefore have a very weak interaction.
The Tersoff potential was developed for modelling of predominantly covalent bonded materials. While this potential is less accurate then ab initio methods, it is far less complicated and thus can be used to model larger and more complex systems. In this potential, the potential energy is modeled as a sum of pair interactions, where the attractive term is dependant on the local enrivonment, giving an effective many-body potential. Thus, we take several parameters for triplets of atoms. The first atom is the central atom, the second is the bonded atom and the third is the atom which influence the two interacting atoms, creating an effective environment. The energy as a function of the atomic coordinates is taken to be
where the potential energy is decomposed unto the site energy Ei and a bonding potential Vij. rij is the distance between the atom i and the atom j. fA and fB are the attractive and repulsive pair potentials respectively, while fC is a smooth cutoff function.
The parameters R and D are chosen to include the first neighbor shell only for a selected structure. the function fC decreases from 1 to 0 in the range of R-D < r < R+D
As stated earlier, this potential takes into account the presence of an environment atom as well, which translates as a variation of the strength of each bond, depanding on the local environment. Thus, we modify the force according to:
The specific parameters for silicon nitride can be found in F.de Brito Mota et al., Phys. Rev. B 58(13), pp. 8323 (1998).
further information on the tersoff potential can be found on Lammps Manual: pair_style tersoff command
The Virial Stress
The virial stress is a measure of the macroscopic stress derived from interatomic potential in molecular dynamics simulations. Taking a control volume Ω we can write an average of the stress (using the virial average, and hence the name virial stress) to get:
where m(i) is the mass of the i-th atom in the control volume, x(i) and u(i) are the position and velocity in the cartesin system.
f(i,j) is the force which the i-th atom exerts on the j-th atom.
The topic of virial stress and the translation of the calculated virial stress to the actual stress appearing in the simulation is still under research today.
Further reading can start at The physical interpertation of the virial stress
In the next page we shall explain how to calculate the virial stress in LAMMPS.