Force fields#

Bonded and non-bonded interactions#

Even with all the approximations and optimizations of the last hundred years, first principles quantum chemical calculations are still very resource demanding, even for quite small systems. The upper limit for a geometry optimization is around a thousand atoms and a dynamics simulation is only possible for systems of a few hundred atoms, and even then only for very short time scales. This means that unless new approximations are introduced, many systems that are of interest will be far beyond the capabilities of today’s computers. A small piece of protein from the human body can consist of several thousand atoms, and even very small systems can increase by hundreds of atoms when a solvent is added. A popular way of dealing with this problem when simulating the motion of atoms is to use molecular mechanics (MM).

In electronic structure theory, the [Born–Oppenheimer] approximation is utilized to separate the Schrödinger equation into electronic and nuclear parts that can be solved separately. The electrons are approximated as responding instantly to any changes in the nuclear coordinates, allowing the electronic distribution to be found around stationary nuclei, after which the problem is solved for the nuclei moving in the effective potential generated by the electrons. MM takes this one step further and combines nuclei and electrons into a single unit, unsurprisingly referred to as atoms. The potential energy of the system is then approximated as a simple analytical function of the atomic coordinates. This potential energy function is split into contributions from different kinds of interactions involving two or more atoms, such as bond stretching and angle bending. The specific form of the potential energy function as well as the parameters used in it are know as a force field. The included interactions and the functional form that they take can differ depending on the requirements of the force field, with some focused on describing specific groups of molecules well and some sacrificing accuracy for computational efficiency. A simple energy function can look like

\[ E = {E_\mathrm{s}} + {E_\theta} + {E_\omega} + {E_\mathrm{vdw}} + {E_\mathrm{el}} \]

where each term is a sum, adding up all the contribution from a certain kind of bonded or non-bonded interaction in the system.

The potential energy function of the system has as its variables the atomic coordinates, but it also depends parametrically on how the atoms are bonded and the types of the involved atoms. In MM, atom type does not just mean which element the atom belongs to, but also its chemical environment. This concept of atom types contains one of the main ideas behind MM force fields, which is that molecular structures that exist in similar environments behave in a similar manner. As an example, the carbon–carbon bond found in ethane is very similar to the two found in propane. The bond lengths are similar and the changes in potential energy when the bonds are stretched or contracted are almost the same. The parameters used in the potential energy function describing these three bonds can therefore safely be approximated as being the same. This transferability is one of the main strengths of using force fields, as parameters developed for a certain structure in one system can be reused for similar structures in other systems. For example, the parameters needed to describe butane are all that is required to describe all linear and branched alkanes in a simple force field.

While the carbon–carbon bonds in ethane and propane are similar enough for the parameters to be transferable, the same cannot be said for the carbon–carbon triple bond in acetylene, where the bond length is much shorter and the bond strength is significantly higher. Using the same parameters for all these bonds would lead to a poor description, which indicates that a finer division of atom types is needed than to just have one for each element. How fine this division is depends on what is desired from the force field. Creating accurate parameters can be a difficult and time-consuming process, with fittings made either to experimental data or higher level quantum mechanical calculations. A larger number of atom types makes it possible to give a more nuanced description of the interactions in the system, but it also strips away the transferability advantages. With more atom types the number of parameters quickly increases, and with them the likelihood that some of them will have to be created when studying a new system.

Force-field terms#

The following section describes the most common force-field interactions, using the [CHARMM] and [MM3] force fields as examples. The two main categories of interaction terms are bonded and non-bonded. Bonded interactions exist for atoms that are part of the same molecule, the structure of which is defined at the start of the simulation by the user. There are also reactive force fields, such as ReaxFF [vanduin2001], where bonds can be formed and broken dynamically. Non-bonded terms appear between atoms in different molecules as well as between atoms in the same molecule, provided that they are separated by at least a certain number of bonds—usually three. For atom pairs connected by fewer bonds, the non-bonded interactions are incorporated into the bonded interactions.

import numpy as np
import py3Dmol as p3d
import veloxchem as vlx
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit
from scipy.interpolate import interp1d

Bond stretching#

Stretching or contracting a bond from its equilibrium position causes the energy of the system to rise.

../../_images/MM_stretch.svg

Fig. 20 Bond stretching.#

This change in energy is approximated by the bond stretch terms, usually represented by a simple harmonic oscillator of the form

\[ {E_\mathrm{s}} = \frac{k_\text{s}}{2} (l - l_0)^2 \]

where \(k_\text{s}\) is a force constant, defining the strength of the bond, \(l\) is the distance between the two bonded atoms and \(l_0\) is the equilibrium distance for the interaction. Both \(k_\text{s}\) and \(l_0\) are constants specific to the combination of atom types involved in the interaction. It should be noted that while \(l_0\) is the bond length that minimizes the equation, it is not necessarily the equilibrium distance found in a molecule. As the two bonded atoms are most likely also involved in several other interactions, the molecular equilibrium distance will be the one that minimizes the full potential.

We will now calculate the B3LYP potenteial energy surface of CO, and se how different fitted force fields can replicate the reference potential.

Note

We will here consider a larger bond distance, in order to more clearly see how different choices of force field fittings affect the general shape of the resulting potential. For practical studies, these large distortions of the CO bond will not be sampled, as 10 kJ/mol corresponds to a temperatuer of about 1200 K. Any regions of the PES above this will thus hardly be sampled for most simulations.

# CO with variable C position
co_base = """
O       0.0       0.0       0.0
C         x       0.0       0.0
"""

bond_e = []
bond_l = np.arange(1.05, 1.30, 0.01)

for d in bond_l:
    # change C position
    co_tmp = co_base.replace("x", "{}".format(d))
    mol_vlx = vlx.Molecule.read_str(co_tmp)
    bas_vlx = vlx.MolecularBasis.read(mol_vlx, "6-31G")

    # calculate B3LYP energy
    scf_settings = {}
    method_settings = {"xcfun": "b3lyp"}
    scf_drv = vlx.ScfRestrictedDriver()
    scf_drv.update_settings(scf_settings, method_settings)
    scf_drv.compute(mol_vlx, bas_vlx)
    bond_e.append(scf_drv.get_scf_energy())

# index for minimal energy
bond_indx = bond_e.index(min(bond_e))
# convert to kJ/mol and shift to zero
bond_e = 2600.0 * (np.array(bond_e) - min(bond_e))

Structures and B3LYP potential energy surface:

plt.figure(figsize=(6,3))
plt.plot(bond_l, bond_e, 'o--')
plt.xlabel('Bond distance [Å]')
plt.ylabel('Energy [kJ/mol]')
plt.tight_layout()
plt.show()

viewer = p3d.view(width=400, height=200)
viewer.setViewStyle({"style": "outline", "color": "black", "width": 0.1})
viewer.addModelsAsFrames(bond_xyz)
viewer.animate({"loop": "forward"})
viewer.setStyle({"stick": {}})
viewer.show()
../../_images/forcefields_6_0.png

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

Harmonic fit#

This will now be fitted using the using three different approaches:

  1. Using the full energy curve calculated above

  2. Using a symmetric distance interval around the minima

  3. Using an interval from the same energies

# function to be fitted to
def harm_func(l, l0, ks):
    return (ks / 2.0) * (l - l0) ** 2

# feed with equilibrium distance and a stiff force constant
# fitting for full curve
opt_harm1 = curve_fit(harm_func, bond_l, bond_e, [bond_l[bond_indx], 1000])
print(
    "l_0 / k_s from fit #1: {:04.3f} / {:06.1f}".format(
        opt_harm1[0][0], opt_harm1[0][1]
    )
)

# fitting for minima +/-0.025 Å
for i in np.arange(len(bond_l)):
    if bond_l[i] < bond_l[bond_indx] - 0.025:
        i_min_d = i
    if bond_l[i] < bond_l[bond_indx] + 0.025:
        i_max_d = i + 1
opt_harm2 = curve_fit(
    harm_func, bond_l[i_min_d:i_max_d], bond_e[i_min_d:i_max_d], [bond_l[bond_indx], 1000]
)
print(
    "l_0 / k_s from fit #2: {:04.3f} / {:06.1f}".format(
        opt_harm2[0][0], opt_harm2[0][1]
    )
)

# fitting for minima +10 kJ/mol
for i in np.arange(len(bond_e)):
    if (bond_e[i] > 10.0) and (i < bond_indx):
        i_min_e = i
    if (bond_e[i] < 10.0) and (i > bond_indx):
        i_max_e = i
opt_harm3 = curve_fit(
    harm_func, bond_l[i_min_e:i_max_e], bond_e[i_min_e:i_max_e], [bond_l[bond_indx], 1000]
)
print(
    "l_0 / k_s from fit #3: {:04.3f} / {:06.1f}".format(
        opt_harm3[0][0], opt_harm3[0][1]
    )
)
print("B3LYP bond distance: {:04.3f} Å".format(bond_l[bond_indx]))

plt.figure(figsize=(6, 3))
x = np.arange(1.06, 1.29, 0.001)
y = interp1d(bond_l, bond_e, kind="cubic")
plt.plot(x, y(x), "k-", lw=2.0) # plot spline B3LYP results
plt.plot(x, harm_func(x, opt_harm1[0][0], opt_harm1[0][1]))
plt.plot(x, harm_func(x, opt_harm2[0][0], opt_harm2[0][1]))
plt.plot(x, harm_func(x, opt_harm3[0][0], opt_harm3[0][1]))
plt.ylim(-1, 40)
plt.xlim((1.05, 1.275))
plt.legend(("B3LYP", "Full", "Distance symm.", "Energy symm."))
plt.xlabel("Bond distance [Å]")
plt.ylabel("Energy [kJ/mol]")
plt.show()
l_0 / k_s from fit #1: 1.167 / 9888.5
l_0 / k_s from fit #2: 1.157 / 10094.0
l_0 / k_s from fit #3: 1.158 / 10536.4
B3LYP bond distance: 1.160 Å
../../_images/forcefields_8_1.png

It is clear that using the full SCF curve we calculated yield unrealistic fits, as it puts weight on parts of the PES that will no be sampled by most dynamics. Using fits with intervals symmetric in either energy or bond distance yield similar results, especially within the region of most interest (\(<10\) kJ/mol).

Anharmonic fit#

Generally, contracting a bond produces a steeper rise in energy than stretching it. This is accounted for in some force fields, such as MM3, which introduce anharmonicity into the bond energy expression by including higher order terms based on the [Morse potential]. The MM3 bond stretching term has the form

\[ {E_\mathrm{s}} = \frac{k_\text{s}}{2} (l - l_0)^2 \left( 1 - \alpha (l - l_0) + \alpha \frac{7}{12} (l - l_0)^2 \right) \]

where \(\alpha\) is derived as 2.55 from the Morse potential, but could be taken as an additional parameter. Performing curve fits for \(l_0\), \(k_s\), and \(\alpha\) yield:

def anharm_func(l, l0, ks, a):
    return (
        (ks / 2.0)
        * (l - l0) ** 2
        * (1 - a * (l - l0) + a * ((7.0 / 12.0) * (l - l0) ** 2))
    )

# feed with equilibrium distance, a stiff force constant, and a=2.55
# fitting for full curve
opt_anharm1 = curve_fit(anharm_func, bond_l, bond_e, [bond_l[bond_indx], 1000, 2.55])
print(
    "l_0 / k_s / alpha from fit #1: {:04.3f} / {:06.1f} / {:04.3f}".format(
        opt_anharm1[0][0], opt_anharm1[0][1], opt_anharm1[0][2]
    )
)

opt_anharm2 = curve_fit(
    anharm_func, bond_l[i_min_d:i_max_d], bond_e[i_min_d:i_max_d], [bond_l[bond_indx], 1000, 2.55]
)
print(
    "l_0 / k_s / alpha from fit #2: {:04.3f} / {:06.1f} / {:04.3f}".format(
        opt_anharm2[0][0], opt_anharm2[0][1], opt_anharm2[0][2]
    )
)

opt_anharm3 = curve_fit(
    anharm_func, bond_l[i_min_e:i_max_e], bond_e[i_min_e:i_max_e], [bond_l[bond_indx], 1000, 2.55]
)
print(
    "l_0 / k_s / alpha from fit #3: {:04.3f} / {:06.1f} / {:04.3f}".format(
        opt_anharm3[0][0], opt_anharm3[0][1], opt_anharm3[0][2]
    )
)

plt.figure(figsize=(6, 3))
plt.plot(x, y(x), "k-", lw=2.0) # plot spline B3LYP results
plt.plot(x, anharm_func(x, opt_anharm1[0][0], opt_anharm1[0][1], opt_anharm1[0][2]))
plt.plot(x, anharm_func(x, opt_anharm2[0][0], opt_anharm2[0][1], opt_anharm2[0][2]))
plt.plot(x, anharm_func(x, opt_anharm3[0][0], opt_anharm3[0][1], opt_anharm3[0][2]))
plt.ylim(-1, 40)
plt.xlim((1.05, 1.275))
plt.legend(("B3LYP", "Full", "Distance symm.", "Energy symm."))
plt.xlabel("Bond distance [Å]")
plt.ylabel("Energy [kJ/mol]")
plt.show()
l_0 / k_s / alpha from fit #1: 1.157 / 10480.0 / 2.234
l_0 / k_s / alpha from fit #2: 1.156 / 10026.5 / 2.747
l_0 / k_s / alpha from fit #3: 1.156 / 10181.5 / 2.671
../../_images/forcefields_10_1.png

Comparing anharmonic and harmonic fits (fitted with energy-symmetric), as well as using force constant:

Note

Calculation of force constant is to be added.

plt.figure(figsize=(6, 3))
plt.plot(x, y(x), "k-", lw=2.0) # plot spline B3LYP results
plt.plot(x, harm_func(x, opt_harm3[0][0], opt_harm3[0][1]))
plt.plot(x, anharm_func(x, opt_anharm3[0][0], opt_anharm3[0][1], opt_anharm3[0][2]))
plt.ylim(-0.5, 10.0)
plt.xlim((1.11, 1.21))
plt.legend(("B3LYP", "Harmonic", "Anharmonic"))
plt.xlabel("Bond distance [Å]")
plt.ylabel("Energy [kJ/mol]")
plt.show()
../../_images/forcefields_12_0.png

Angle bending#

Same as for bond stretching, distorting an angle formed between two atoms bound to a common third atom from its equilibrium causes the energy of the system to rise. This too is commonly approximated by a harmonic oscillator of the form

\[ {E_\theta} = \frac{k_\theta}{2} (\theta - \theta_0)^2 \]

where \(k_\theta\) is a force constant, \(\theta\) is the angle formed by the atoms and \(\theta_0\) is the equilibrium angle.

../../_images/MM_bend.svg

Fig. 21 Angle bending.#

We do the fitting for the bond angle of water. As a scan creates output files which are then to be read, we instead “scan” one angle at a time and save the results:

# Set up the molecule and basis set
water_xyz = """
O       0.0000000000     0.0000000000     0.1178336003
H      -0.7595754146    -0.0000000000    -0.4713344012
H       0.7595754146     0.0000000000    -0.4713344012
"""

# Create veloxchem mol and basis objects
mol_vlx = vlx.Molecule.read_str(water_xyz)
bas_vlx = vlx.MolecularBasis.read(mol_vlx, "6-31G")

scf_settings = {}
method_settings = {}
scf_drv = vlx.ScfRestrictedDriver()
scf_drv.update_settings(scf_settings, method_settings)
scf_drv.compute(mol_vlx, bas_vlx)

# Set up the gradient and optimization dirvers:
grad_drv = vlx.ScfGradientDriver(scf_drv)
opt_drv = vlx.OptimizationDriver(grad_drv)

angle_a = np.arange(60,180,5)
angle_e = []
angle_xyz = []
for angle in angle_a:
    angle_scan = "angle 2 1 3 {} {} 1".format(angle,angle) 
    opt_drv.constraints = ["$scan", angle_scan]

    # Run the scan
    opt_molecule = opt_drv.compute(mol_vlx, bas_vlx)
    angle_xyz.append(opt_molecule.get_coordinates())
    angle_e.append(scf_drv.get_scf_energy())

# index for minimal energy
angle_indx = angle_e.index(min(angle_e))
# convert to kJ/mol and shift to zero
angle_e = 2600.0 * (np.array(angle_e) - min(angle_e))
plt.figure(figsize=(6,3))
plt.plot(angle_a, angle_e, 'o--')
plt.xlabel('Bond angle')
plt.ylabel('Energy [kJ/mol]')
plt.tight_layout()
plt.show()


angle_struct = ''
b2a =  0.5291772
for i in np.arange(len(angle_a)):
    angle_struct += '3\n*\n'
    angle_struct += 'O ' + str(np.around(b2a*angle_xyz[i][0][0],5)) + ' ' + str(np.around(b2a*angle_xyz[i][0][1],5)) + ' ' + str(np.around(b2a*angle_xyz[i][0][2],5)) + '\n'
    angle_struct += 'H ' + str(np.around(b2a*angle_xyz[i][1][0],5)) + ' ' + str(np.around(b2a*angle_xyz[i][1][1],5)) + ' ' + str(np.around(b2a*angle_xyz[i][1][2],5)) + '\n'
    angle_struct += 'H ' + str(np.around(b2a*angle_xyz[i][2][0],5)) + ' ' + str(np.around(b2a*angle_xyz[i][2][1],5)) + ' ' + str(np.around(b2a*angle_xyz[i][2][2],5)) + '\n'

viewer = p3d.view(width=400, height=200)
viewer.setViewStyle({"style": "outline", "color": "black", "width": 0.1})
viewer.addModelsAsFrames(angle_struct)
viewer.animate({"loop": "forward"})
viewer.setStyle({"stick": {}})
viewer.rotate(-90, "x")
viewer.show()
../../_images/forcefields_15_0.png

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol

Harmonic fit#

Fitting with a harmonic function, focusing on an interval with equal energy end-points:

# fitting for minima +20 kJ/mol
for i in np.arange(len(angle_e)):
    if (angle_e[i] > 20.0) and (i < angle_indx):
        i_min_a = i
    if (angle_e[i] < 20.0) and (i > angle_indx):
        i_max_a = i
opt_harm_angle = curve_fit(
    harm_func, angle_a[i_min_a:i_max_a], angle_e[i_min_a:i_max_a]
)
print(
    "l_0 / k_s from fit: {:04.3f} / {:7f}".format(
        opt_harm_angle[0][0], opt_harm_angle[0][1]
    )
)
print("B3LYP angle: {:04.1f}".format(angle_a[angle_indx]))

plt.figure(figsize=(6, 3))
x = np.arange(75,150, 0.1)
y = interp1d(angle_a, angle_e, kind="cubic")
plt.plot(x, y(x), "k-", lw=2.0) # plot spline B3LYP results
plt.plot(x, harm_func(x, opt_harm_angle[0][0], opt_harm_angle[0][1]))
plt.ylim(-1,20)
plt.xlim((90,135))
plt.legend(("B3LYP", "Harmonic fit"))
plt.xlabel("Angle")
plt.ylabel("Energy [kJ/mol]")
plt.show()
l_0 / k_s from fit: 112.149 / 0.134355
B3LYP angle: 110.0
../../_images/forcefields_17_1.png

The real angle bending interaction is clearly not a harmonic oscillator, as exemplified by the potential energy curve in below Figure, showing the change in energy when distorting the H–O–H angle found in water. Compressing the angle causes the two hydrogen atoms to come close together, leading to a steep curve near \(\theta = 0^\circ\), while increasing the angle leads to a cusp at \(\theta = 180^\circ\) due to symmetry. Again, some force fields account for this anharmonicity with additional terms, such as the MM3 angle bending term:

\[\begin{equation*} \begin{split} {E_\theta} = \frac{k_\theta}{2} (\theta - \theta_0)^2 \bigl( 1 &- \alpha (\theta - \theta_0) + \beta (\theta - \theta_0)^2 \\ & - \gamma (\theta - \theta_0)^3 + \delta (\theta - \theta_0)^4 \bigr) \end{split} \end{equation*}\]

where the constants have been fitted to experimental data but are the same for all bending interactions. This gives a better agreement with the actual angle bending behaviour, as seen in Figure below, but again, for lower temperatures the harmonic oscillator is still a valid and computationally efficient approximation.

This gives a better agreement with the actual angle bending behaviour, as seen in Figure below, but again, for lower temperatures the harmonic oscillator is still a valid and computationally efficient approximation.

../../_images/MM_waterbend_even.svg

Fig. 22 Potential energy curves for distortion of the water angle. The reference ab initio energy was calculated at the MP2/cc-pVTZ level of theory while the other curves were derived from the harmonic oscillator used in CHARMM and the anharmonic expression used in MM3.#

Torsion#

When two atoms are bonded to atoms on opposite sides of a central bond, as for the opposing hydrogen atoms in ethane, then twisting these around the central bond can cause a change in the potential energy. The torsional energy is thus a function of the dihedral angle, \(\omega\), formed between the first and fourth atom, while looking through the bond between the second and third atoms. For a visualization of this, see Figure. The torsion energy term is usually expressed in terms of a Fourier expansion of the form

\[ {E_\omega} = k_\omega \left(1 + \cos(n\omega -\delta)\right) \]

where \(k_\omega\) is a force constant, \(n\) is the periodicity and \(\delta\) is a phase shift. Typically, several such terms are included for each set of four atoms in order to combine different multiplicities. The MM3 force field always includes terms of multiplicity 1, 2 and 3, with \(\delta\) set to zero for each, while the CHARMM force field allows a more varied combination of terms.

Figure shows the potential energy curve when twisting butane around the central carbon–carbon bond. The two identical lower peaks represent the conformations when the one hydrogen atom on each side is eclipsed by the carbon atom on the other side, while the main peak comes from the conformation when all atoms are eclipsed by an equivalent atom on the other side. It should be noted that the torsional energy of the force field does not come from a single interaction in this case, but nine. One for the two carbon atoms, four for the carbon atoms meeting the hydrogen atoms on the other side and four for the combinations of hydrogen atoms on opposite sides. It is easy to see that torsional interactions can quickly become very complicated when there are many atom types involved. This is further complicated by the fact that most force field also include non-bonded interactions between atoms separated by more than two bonds, meaning that the true rotational barrier will be a sum of all torsional and non-bonded interactions.

The conformation of a molecule is defined by its current dihedral angles, with two identifiable conformations for butane visible as minima in below Figure: the global minimum at 180 degrees and two identical minima at \(\pm70\) degrees. For calculations where both the conformations themselves as well as the transitions between them are important, it is necessary to recreate the full potential energy curve and not just focus on the regions around the minima, as was done for bond stretching and angle bending terms.

../../_images/MM_torsion.svg

Fig. 23 Torsion.#

../../_images/MM_butanetorsion_even.svg

Fig. 24 Potential energy curves for the torsional barrier found when twisting butane around the central bond. The reference ab initio energy was calculated at the MP2/cc-pVTZ level of theory while the MM curve was obtained as a least squares fit with multiplicities 1, 2, 3, 4 and 5.#

# Set up the molecule and basis set
ethane_xyz = """
C       0.0000000000     0.0000000000     0.7616375820
C       0.0000000000     0.0000000000    -0.7616375820
H       1.0145603421     0.0000000000     1.1557300969
H      -0.5072801711    -0.8786350299     1.1557300969
H      -0.5072801711     0.8786350299     1.1557300969
H       0.5072801711    -0.8786350299    -1.1557300969
H       0.5072801711     0.8786350299    -1.1557300969
H      -1.0145603421    -0.0000000000    -1.1557300969
"""

# Create veloxchem mol and basis objects
mol_vlx = vlx.Molecule.read_str(ethane_xyz)
bas_vlx = vlx.MolecularBasis.read(mol_vlx, "STO-3G")

scf_settings = {}
method_settings = {}
scf_drv = vlx.ScfRestrictedDriver()
scf_drv.update_settings(scf_settings, method_settings)
scf_drv.compute(mol_vlx, bas_vlx)

# Set up the gradient and optimization dirvers:
grad_drv = vlx.ScfGradientDriver(scf_drv)
opt_drv = vlx.OptimizationDriver(grad_drv)

dihed_a = np.arange(0,65,2.5)
dihed_e = []
dihed_xyz = []
for dihedral in dihed_a:
    dihedral_scan = "dihedral 3 1 2 6 {} {} 1".format(dihedral,dihedral) 
    opt_drv.constraints = ["$scan", dihedral_scan]

    # Run the scan
    opt_molecule = opt_drv.compute(mol_vlx, bas_vlx)
    dihed_xyz.append(opt_molecule.get_coordinates())
    dihed_e.append(scf_drv.get_scf_energy())

# index for minimal energy
dihed_indx = dihed_e.index(min(dihed_e))
# convert to kJ/mol and shift to zero
dihed_e = 2600.0 * (np.array(dihed_e) - min(dihed_e))

plt.figure(figsize=(6,3))
plt.plot(dihed_a, dihed_e,'o--',color='#1f77b4')
plt.plot(-dihed_a, dihed_e,'o--',color='#1f77b4')
plt.xlabel('Dihedral angle')
plt.ylabel('Energy [kJ/mol]')
plt.xlim((-60,60))
plt.tight_layout()
plt.show()
../../_images/forcefields_22_0.png

Electrostatic interaction#

The non-bonded force field interactions are divided into electrostatic interactions, caused by static multipoles, and van der Waals interactions, containing everything else. For most force fields, only point charges are used in the electrostatic interactions, resulting in an energy given simply by the Coulomb interaction

\[ {E_\mathrm{el}} = \frac{q_i q_j}{4\pi \epsilon_0 r_{ij}}, \]

where \(q_i\) and \(q_j\) are the partial charges of the interacting atoms, \(r_{ij}\) is the distance between them and \(\epsilon_0\) is the vacuum permittivity. It is also possible to use higher order multipoles, an example of which can be found in the MM3 force field, where dipole moments are assigned to each bonded pair. This gives the advantage of a more flexible description, as the dipoles depend on the types of both involved atoms, but requires additional parametrization and has difficulty dealing with ionic systems.

van der Waals interaction#

The van der Waals interaction is typically represented by two terms. One due to exchange interaction, causing a strong repulsion at close distances, and one that is due to London dispersion, resulting in a weak attraction at longer distances. The most commonly used form is the [Lennard-Jones potential]:

\[ {E_\mathrm{vdw}} = \epsilon \left[ \left( \frac{r_\text{min}}{r} \right)^{12} - 2 \left( \frac{r_\text{min}}{r} \right)^6 \right] \]

\noindent where \(\epsilon\) is the depth of the potential well, usually created as the mean of two parameters specific to the two interacting atom types, \(r_\text{min}\) is the position of the bottom of the well and \(r\) is the interatomic distance. While it is possible to motivate the use of \(r^{-6}\) for the attractive part, the use of \(r^{-12}\) for the repulsive part is more pragmatic in origin, having that form simply because it is the square of \(r^{-6}\), making its computation highly efficient. It can be argued[cite jensen2006] that an exponential form is more realistic, such as the [Buckingham potential]used in MM3:

\[ {E_\mathrm{vdw}} = A e^{- B r} - \frac{C}{r^{-6}} \]

where \(A\), \(B\), and \(C\) are constants. Calculating this expression is more demanding than the Lennard-Jones (LJ) potential, but produces a repulsive force that is closer to what is found using \emph{ab initio} methods. For very short interatomic distances, however, the Buckingham potential becomes strongly attractive, which can lead to highly unphysical behaviour for systems in unfavourable starting positions.

Additional terms#

Depending on what the force field is intended for, many additional terms can be added to the potential energy function. Cross-terms are a common addition, accounting for more complex couplings between structural properties. Examples of these can be found in the MM3 force field, which contains stretch–bend, torsion–stretch and bend–bend interactions. As an example, the stretch–bend interaction accounts for the fact that compressing an angle formed by three bonded atoms will cause the repulsion between the two end atoms to increase, resulting in an elongation of their bonds to the middle atom. Another example of such structural corrections comes in the form of the out-of-plane term. If three or four atoms are bonded to a central atom with their equilibrium position in a plane, as for e.g.\ formaledehyde, the change in energy caused by a distortion of the central atom out of the plane is hard to model by just angle bending terms, as the individual angle distortions are small. The out-of-plane term instead depends on the angles formed between the plane and the bonds connecting the central and outer atoms. This is shown in below Figure, where the \(\chi\) is an example of an out-of-plane angle.

../../_images/MM_outofplane.svg

Fig. 25 Out-of-plane bending.#