1.3 Study aims
2.1.2 Force contributions and constraints
Listing 1: Velocity Verlet pseudo code. The momentum p = mv as the product of mass mand velocityv, the atoms position r, force f, time step dtand a functionforce(r)to caclulate the force acting onr.
and would require very small time-steps to faithfully represent their dynamics. Therefore, algorithms likeLINCS [126] andSETTLE [127]
have been developed which can guarantee that constraints on bond length in particular and atom to atom distances in general are ful-filled at the end of a simulation time-step. Also performed at this step is pressure coupling (to a virtualpressure bath), often by scaling the simulation box and the atom positions so that the system con-forms to a target pressure. Pressure and temperature coupling are needed to simulate an NPT ensemble where the number of atoms (N), the pressure (P) and temperature (T) are constant. Finally, the atoms position are written out to create the trajectories that are later used for analysis.
2.1.2 Force contributions and constraints
MD simulations must consider a number of non-bonded and bonded (intra-molecular) interactions to sufficiently approximate the dynamic behavior of the simulated atoms over time. However, apart from describing the forces that act on an atom, being able to compute them fast and efficiently is a prime objective when choosing the equation and models to calculate them.
2.1 molecular dynamics simulations 21
Figure6:The Lennard-Jones potential (left) is plotted for a carbon-carbon interaction with parameters from theamber99sb-ildnforce-field [129]. The Coulomb interaction (right) shows the potential for particles with equal ele-mentary charges (blue) and with the inverse sign (red).
Non-bonded interactions
The Lennard-Jones potential describes the neutral (uncharged) inter-action between a pair of atoms as a function of their distance. It consists of a weakly attractive term dominant at larger distances and a strongly repulsive term that becomes dominant when the atoms ap-proach each other (Equation 4, Figure 6). The attractive long-range term represents the van der Waals force and the short-range repul-sive term approximates the Pauli-repulsion created by overlapping electron orbitals. The Lennard-Jones potential has two parameters,ε and σ, which are specific for the type of interacting atoms. These parameters are created by fits to experimental data and provided by force-fields. Other potentials for describing neutral atom interac-tions, the Buckingham potential or the Stockmayer potential, provide a more realistic description but are more expensive to compute than the Lennard-Jones potential [128].
vLJ(ri,rj) =4ε
σ dij
12
´
σ dij
6
Equation4:The Lennard-Jones potential is used to model the neutral inter-action between a pair of atoms. The parametersεandσare specific for the types of the interacting atomsriandrj.
The Coulomb potential describes the interaction between a pair of atoms resulting from electrical charges as a function of the distance between the atoms. All atoms in charged or uncharged molecules have individual partial charges resulting in an attractive potential for charges with opposing signs and a repulsive potential for atoms with
the same signs (Equation5, Figure6). The parametersQiandQj are the partial charges of the interacting atoms and are usually calculated for force-fields by solving quantum mechanical equations for smaller parts of whole molecules.
vCoulomb(ri,rj) = QiQj
4πϵ0dij
Equation5:The Coulomb potential is a pair potential between point charges in vacuum with a permittivity of ϵ0. The parameters Qi and Qj are the charges of the interacting atoms ri and rj; dij is their distance from each other.
The non-bonded interactions must be evaluated for all pairs of atoms when applied strictly (Equation 6), which is computationally very expensive as there are 12N(N´1) such interactions in a system with N atoms. However, a closer look at the interaction potentials suggests that not all atom pairs need to be considered all the time.
At increasing distance, the forces exerted by the non-bonded inter-actions quickly get smaller. The Lennard-Jones potential approaches zero at an exponential rate in the order ofd12ij and the Coulomb po-tential falls off quadratically with respect to the inter-atom distance.
To exploit this rapid decline in interaction strength the neighbor-list data structure has been developed. Such a data structure is created for each atom (central atom) and includes all atoms that may inter-act with the central atom for a given time span. During that time, only the atoms in the neighbor-list are used when calculating the non-bonded interactions for the central atom. After the given time span, the neighbor-list is updated.
Unon-bonded=
atompairs
ÿ
i,j
vLJ(ri,rj) +vCoulomb(ri,rj)
Equation6:The non-bonded interaction potential is the sum of the Coulomb vCoulomband Lenard-JonesvLJpotentials between all pairs of atomsi,j.
Intra-molecular, bonded interactions
In addition to non-bonded interactions, there are a number of bonded interactions, that characterize the topology of the simulated molecules. These intra-molecular interactions are described by bond-stretching, angle-bending and dihedral-torsion potentials between pairs, triples and quadruples of atom positions, respectively (Equa-tion7). The actual constraining is performed during the simulation
2.1 molecular dynamics simulations 23
Figure7:Schematic display of bonded interactions for (A) bond-stretching between pairs of atoms, (B) angle-bending in atom triples and (C) dihedral-torsion in atom quadruples.
loop by algorithms like LINCS [126] for general molecules and by SETTLE [127] optimized for water molecules. They change the atom positions to guarantee that constrained atom to atom distances are fulfilled at the end of a simulation time-step. Using these distance constraints, the bond length can be constraint but also bond and di-hedral angles by expressing the angular constraints with multiple distance constraints.
Equation 7: The intra-molecular, bonded interactions are calculated be-tween pairs, triples and quadruples of atom positionsrrepresenting bond-stretching, angle-bending and dihedral-torsion potentials respectively (Fig-ure7).