Dihedral angle bond has singularities

[jngrad]

ESPResSo uses a chain rule that introduces a term $1 / sin(\theta)$ in the expression of the dihedral force. When the multiplicity is an integer (which is required for periodicity), an alternative chain rule could be used to avoid this singularity.

Relevant literature:

• Swope, W.C. and Ferguson, D.M. (1992), Alternative expressions for energies and forces due to angle bending and torsional energy. J. Comput. Chem., 13:585-594. doi:10.1002/jcc.540130508
• Bekker, H., Berendsen, H.J.C. and van Gunsteren, W.F. (1995), Force and virial of torsional-angle-dependent potentials. J. Comput. Chem., 16:527-533. doi:10.1002/jcc.540160502
• Ailan Cheng and Kenneth M. Merz (1996), The Pressure and Pressure Tensor for Macromolecular Systems. The Journal of Physical Chemistry, 100(2):905-908. doi:10.1021/jp952646x

import espressomd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter, MultipleLocator
plt.rcParams.update({'font.size': 12})

def rotate_vector(v, k, phi):
    vrot = np.array(v) * np.cos(phi) + np.cross(k, v) * \
        np.sin(phi) + np.array(k) * np.dot(k, v) * (1.0 - np.cos(phi))
    return vrot

system = espressomd.System(box_l=[10.0, 10.0, 10.0])
system.cell_system.skin = 0.4
system.time_step = 0.1
axis = np.array([1., 0., 0.])
dh_k = 2.
N = 400  # even number to get singularities at phi=0 and phi=pi
d_phi = 2 * np.pi / N
dh_n, dh_phi0_div = (2, 3)
dh_phi0 = np.pi / dh_phi0_div
# tabulated values for the range [0, 2*pi]
phi = d_phi * np.arange(N + 1)
tab_energy = dh_k * (1. - np.cos(dh_n * phi - dh_phi0))
div = np.sin(phi)
div[0] = div[N // 2] = div[N] = 1.
tab_force = -dh_k * dh_n * np.sin(dh_n * phi - dh_phi0) / div
tab_force[0] = tab_force[N // 2] = tab_force[N] = 0.
dihedral_tabulated = espressomd.interactions.TabulatedDihedral(
    energy=tab_energy, force=tab_force)
system.bonded_inter.add(dihedral_tabulated)
p0, p1, p2, p3 = system.part.add(pos=4 * [(0., 0., 0.)])
p1.add_bond((dihedral_tabulated, p0, p2, p3))
p1.pos = [5., 5., 5.]
p2.pos = p1.pos + [1., 0., 0.]
p0.pos = p1.pos + [0., 1., 0.]

# use half the angular resolution to observe interpolation
xdata = []
ydata_energy = []
ydata_force0 = []
for i in range(2 * N - 1):
    phi = i * d_phi / 2.
    p3.pos = p2.pos + rotate_vector([0., 1., 0.], axis, phi)
    system.integrator.run(recalc_forces=True, steps=0)
    xdata.append(phi)
    ydata_energy.append(system.analysis.energy()["bonded"])
    ydata_force0.append(np.linalg.norm(p0.f))

xdata_phi = np.concatenate([np.array(xdata) - 2 * np.pi, np.array(xdata)])
plt.plot(xdata_phi, np.tile(ydata_force0, 2), label=r"$\left|\vec{f_0}\right|$")
plt.plot(xdata_phi, np.tile(ydata_energy, 2), label="energy")
ax = plt.gca()
ax.xaxis.set_major_formatter(FuncFormatter(
   lambda val,pos: r"{:.0g}$\pi$".format(val/np.pi) if val != 0 else "0"))
ax.xaxis.set_major_locator(MultipleLocator(base=np.pi))
plt.xlabel(r"Dihedral angle $\phi$")
plt.ylabel("Force or energy (simulation units)")
plt.legend()
plt.show()