The base of any molecular dynamics program is the empirical force field that is used to calculate the interactions between the atoms. Most commonly these interactions are split into two groups, bonded and non-bonded interactions. In the former group all intramolecular interactions (such as bonds, angles and dihedrals) are being taken care of, whereas in the latter interactions such as Van der Waals and electrostatics are considered. In this section all potentials available in PumMa are discussed, starting with the intramolecular potentials.

Bonded potential

To model a covalent bond in a molecular structure, many types of interaction potentials can be used, such as the Morse potential or the finitely extendible nonlinear elastic (FENE) potential. However, the most common potential to be used in any molecular dynamics program is the harmonic bond potential. The equation that describes the potential energy of the harmonic potential is given by

In this equation r0 is the reference bond length, kij the force constant, and rij the distance between atoms i and j, hence the current length of the bond. The figure to the right schematically shows some of these parameters. The harmonic potential basically is a Taylor approximation of more elaborate potentials around the reference bond length. From this potential the force acting on atom i can be deduced (assuming the we are dealing with central and conservative forces) through

Substituting Equation (1) into (2) gives us for the force acting on atom i


The force acting on atom j is, due to Newton's third law, of equal magnitude but opposite in direction to the force acting on atom i.

Angle potential

Between two bonds sharing a common atom it is possible for an angle to be defined. Such an angle is shown in the figure on the right for the set of three consecutive atoms i, j and k. The corresponding bond angle is indicated by θijk. Based on the position vectors of the atoms, the bond vectors for the two bonds can be defined as rij=ri-rj and rkj=rk-rj. It is important to notice that both bond vectors are defined towards the atom that they share (in this case atom j).

From these definitions, the angle θijk can subsequently be expressed as


In this equation rij and rkj are the lengths of the bond vectors, whereas rij and rkj are the vectors, and · denotes the scalar product. Similar to the derivation of the force with the bond potential, the force acting on atom i is given by


where the chain rule has been used to split the derivative with respect to the position vector into two parts. Using Equation (4) the rightmost derivative is easily obtained, whereas the derivative of the potential energy depends on the functional form of the potential energy function.

Two types of angle potentials are most commonly found in a molecular dynamics program. These are the harmonic and the cosine harmonic potential functions, where the former is used in the CHARMM forcefield and the latter in the GROMOS96 forcefield. PumMa can use both functional forms, however not at the same time. Prior to the simulation a choice has to be made on which functional form of the angle potential is to be used.

The cosine harmonic angle potential is given by


where θ0 is the reference angle and kCHijk is the force constant for the cosine harmonic functional form. In a similar way to the bond potential, the harmonic angle potential is expressed as


where θ0 is again the reference angle and kHijk is the force constant for the harmonic version of the angle potential. It is very important to note that the force constants for either the cosine harmonic or the harmonic potential have a different order of magnitude and different units (see the manual pages).

Combining either Equation (6) or (7) with (5) gives the expression for the force on atom i, and similar expressions can be obtained for the forces on atoms j and k.


In the CHARMM force field an additional potential function is added to the angle, in order to restrain the motions of the bonds involved in the angle. This potential function is often referred to as the Urey-Bradley potential. The functional form of this potential looks very similar to the bond potential function and is given by


Here the reference length of the distance between particles i and k is given by r0 and the corresponding force constant is given by kUBik. Basically the Urey-Bradley potential introduces a virtual bond between atoms i and k. The force associated with this potential can be derived similar to the bond potential.

Torsion potentials

Two types of torsion potentials are most commonly distinguished: dihedral angle potentials and improper torsions. Both potentials rely on a quartet of atoms, bonded in one way or the other. A dihedral angle potentials depends on four consecutive bonded atoms, whereas the improper torsion depends on three atoms centered around a fourth atom. The dihedral angle potential is mostly used to constrain the rotation around a bond. The improper torsion is used to maintain chirality on a tetrahedral extended heavy atom or to maintain planarity of certain atoms. The main difference between both torsion potentials is the definition of the torsional angle and the functional form of the potential function. Fortunately, the difference in definition of the torsional angle can be eliminated by numbering the atoms in a clever way. First the proper dihedral is discussed and subsequently the improper type.

Proper dihedrals

 In the figure on the left the definition of the torsional angle for the proper dihedral angle potential can be seen. The torsional angle φijkl is the angle between the plane going through the atoms i, j and k and the plane going through the atoms j, k and l. For each of these two planes the normal vectors can be calculated by taking the cross product of two vectors lying in this plane, where the directions of the vectors have to be chosen in a clever way to get the torsional angle in accordance with the IUPAC/IUB convention.

For the plane going through the atoms i, j and k this normal vector m is expressed as


where rij and rjk are the bond vectors lying in the plane ijk. Similar for the plane going through the atoms j, k and l, the normal vector n is given by


For the torsional angle a similar definition as in (4) can be used. The cosine of the torsional angle φijkl is expressed as


and the sine of the same angle φijkl as


where m, n, and rjk are the lengths of the respective vectors. Using the two definitions for the cosine and sine and taking the IUPAC/IUB convention into account, the torsional angle is given by


To model the rotation barriers around bonds in accordance with thermodynamic data, the cosine form of the dihedral potential was introduced, which is expressed as


where kCijkl is the force constant belonging to the cosine type of potential, φ0 the angle where the potential passes through its minimum value, and nijkl is the multiplicity, which indicates the number of minima as the bond is rotated through 360o. The multiplicity is a nonzero, positive integer number. However, it is not uncommon that the rotation around a bond has local and global minimums. In order to accomodate this type of behavior the potential is split in its underlying harmonic functions, each having its own force constant, multiplicity and reference angle. In PumMa each harmonic contribution to the dihedral potential of a specific quartet of atoms has to be specified separately in the force field parameter set.

Besides the cosine functional form also the harmonic form is available in PumMa. It has to be mentioned that this dihedral potential is only valid when the corresponding dihedral has only one minimum. Whether this is the case in fully atomistic simulations is arguable, but it is believed that it models dihedrals in coarse grained simulations well. The harmonic dihedral potential is expressed as


where kHijkl is the corresponding force constant, which obviously differs from the force constant in the cosine functional form. The reference dihedral angle φ0 should be between -π and π.

The force that is excerted due to any of the two types of functional forms of the proper dihedral potential can be expressed by



For the improper torsion the torsional angle definition is shown in the figure on the right. The angle φijkl still depends on the same two planes ijk and jkl, since the atoms have been chosen in the clever way as can be seen in the figure with the atom i in the center instead on one of the ends of the dihedral chain. By using this definition the same equations (9)-(13) can be used to compute the torsional angle, and, consequently, not many changes have to be applied to molecular dynamics code.

Since the improper torsion potential is mainly used to maintain planarity in a molecular structure. Hence, it only has one minimum and a harmonic potential may be used. Therefore the improper torsion functional form is given by


The derivation of the force based on this potential is again given by (16).

Van der Waals potential

Van der Waals interactions are often referred to as the combination of attractive and repulsive forces between two atoms, which are not bonded to each other. The energy arising from this interaction varies with the separation distance between the two atoms. This energy is zero at infinite distance, but as the separation is reduced the energy decreases, passing through a minimum and from there on increasing rapidly.

The most common potential to model the Van der Waals interactions is known as the Lennard-Jones potential, which depends only on two parameters and is expressed as


where εij is the minimum (well depth) of the potential for the interaction between atom i and j, σij the collision diameter (the separation for which the energy is zero) and rij=|rirj| the separation distance. The attractive part of the potential (the part which contains the power 6) has been experimentally validated. For the repulsive part different powers are suitable, however, the power 12 is usedmost often, since as a result the potential can be more easily computed, being the square of the attractive part. Commonly the parameters (such as σij and εij) are not given for a specific pair, but as a parameter of the atom itself. To be able to calculate the Van der Waals interactions the Lorentz—Berthelot mixing rules are applied to determine the values for the parameters, which are then given by


where σii, σjj, εii and εjj are the collision diameters for the atoms i and j and the well depths for the atoms i and j, respectively. In PumMa, however, not the parameters for an atom type are given (as is common in CHARMM), but the parameters for a pair of atoms, see the parameter file for more information.

The force excerted by the Lennard-Jones potential on, for instance, atom i, can be derived through a similar expression as equation (2), and is expressed as


Electrostatic potential

 Electronegative atoms attract electrons more than less electronegative atoms, giving rise to an uneven distribution of charge in a molecule. This distribution can be represented in many ways, but one of the most common ways is to represent the charge distribution as point charges localised throughout the molecule, mostly coinciding with the nuclei of the atoms. For ions their charge is also considered to be a point charge.

The electrostatic interaction (also known as Coulomb interaction) between two atoms, in the same or a different molecule, is expressed by using Coulombs law


where qi and qj are the charges of atoms i and j respectively, ε0 is the permittivity of vacuum (ε0=8.8542·10-12 C2N−1m−2), εR the dielectric constant (εR=1 for a vacuum by definition) and rij is the distance between atoms i and j. Analogous to the derivation of the force with the Van der Waals potential we can write for the force excerted on atom i as an effect of the Coulomb potential


<< | Theory | Ensembles >>