Note: Smooth torsional potentials for degenerate dihedral angles

Torsional potentials impose steric or geometric constraints between four bonded particles, e.g., the rotational barriers between conformations of the carbon backbone in atomistic models for hydrocarbons.^1,2 Typically, torsional potentials are functions of the dihedral angle $ϕ$ between a plane A containing three particles at {r₁, r₂, r₃} and a plane B containing particle positions {r₂, r₃, r₄} (Fig. 1). The planes A and B have normal vectors $nA=r12×r23$ and $nB=r23×r34$⁠, respectively, where r_ij = r_j − r_i. The dihedral angle between the planes is³ 

In most cases, this definition can be straightforwardly applied to compute the torsional energy and forces.⁴ However, a degenerate state exists when {r₁, r₂, r₃} or {r₂, r₃, r₄} are collinear, and it is no longer possible to uniquely construct the required planes. The dihedral angle and torsional potential are then undefined, and the torsional forces are singular. Near this degenerate state, particles can be subjected to large forces, resulting in numerical integration errors or a complete crash for a molecular dynamics simulation.³ The singular nature of the torsional forces may also lead to unphysical motion or undesired restriction of conformations.

Fortunately, such degenerate conformations are rarely, if ever, observed in atomistic models for particles participating in torsional potentials because the bond angle $𝜃ijk$ between them (Fig. 1),³ 

is restricted by additional interactions to values far from the collinear conformations. However, some coarse-grained polymer models impose significantly softer interactions,^7–9 leading to a finite probability of observing a degenerate conformation, $𝜃ijk=180°$⁠. It is necessary to somehow modify the torsional potential in such a case in order to obtain physically, well-defined behavior in molecular dynamics simulations.

The steric interactions approximated by a torsional potential are expected to weaken when the end particles are farther apart and should be negligible when $??ijk=180°$⁠. Singular torsional forces can then be removed by modifying the potential so that forces go smoothly to zero when $𝜃ijk$ approaches $180°$⁠. Other authors have demonstrated the effectiveness of this approach for certain polymer models.^10–13 However, the functional forms of the potentials chosen in their works offer little control over the range of bond angles for which the potential is modified, which may be problematic for some geometries.

In this note, we propose a simple, general approach for smoothing a torsional potential as a function of the bond angle. The proposed method can be applied to existing torsional potentials over a controllable smoothing range and gives well-defined behavior in the degenerate case. Motivated by a common method for smoothing truncated pair potentials in molecular simulations, we define a smoothed torsional potential,

where $V0(ϕ)$ is the unsmoothed potential, V_s is the value to which the potential is smoothed, and s is a smoothing function. The force on particle n is

where $F0,n(ϕ)$ is the force on particle n in the absence of smoothing and $∇n$ is the gradient with respect to r_n. The smoothing function must satisfy the following conditions:

1. $s(cos𝜃s)=1$ at the onset of smoothing (angle $𝜃s$⁠);

2. $s(cos𝜃c)=0$ at the end of smoothing (angle $𝜃c$⁠);

3. the first derivatives of s are zero at $𝜃s$ and $𝜃c$ so that F_n = F_0,n at $𝜃s$ and $Fn=0$ at $𝜃c$⁠; and

4. s must decay sufficiently quickly so that V and F_n have well-defined limits for all bond angles.

A suitable polynomial for s is

where $xs=cos𝜃s$⁠, $xc=cos𝜃c$⁠, and $m≥2$⁠. (In writing s piecewise, we assume without loss of generality that $𝜃c>𝜃s$⁠.) The values of $𝜃s$⁠, $𝜃c$⁠, and m set the range and steepness of smoothing. The gradients of s can be computed by chain rule, $∇ns(cos𝜃ijk)=s′(cos𝜃ijk)∇ncos𝜃ijk$⁠, where

and

Here $∇n⊗rij$ denotes the outer product and is equal to the identity matrix I when n = j, −I when n = i, and the zero matrix 0 otherwise. In the degenerate limit, the unsmoothed torsional forces⁴ diverge as $1∕sin𝜃ijk$⁠. If $𝜃c=0°$ or $180°$⁠, $s′∼(1−cos2𝜃ijk)m−1∼sin2(m−1)𝜃ijk$⁠, and F_n has a well-defined limit for $m≥2$⁠.

The proposed smoothing function was tested for a simple torsional potential, $V0(ϕ)=−cosϕ∕2+cos2ϕ$⁠. The dihedral angle was initially fixed at $ϕ=120°$ with r₁ = (−1, 0, 1), r₂ = (−1, 0, 0), r₃ = (0, 0, 0), and $r4=(0,−sinϕ,cosϕ)$⁠. The angle $𝜃234$ was then varied from its initial value of $90°$ by rotating r₄ around r₃ in plane B, causing the dihedral angle to jump to $−60°$ as $𝜃234$ passed through $180°$⁠. Figure 2 shows the z-component of the torsional force on particle 4, F_4,z, as a function of $𝜃234$ for different smoothing ranges. As expected, the unsmoothed force diverges at $𝜃234=180°$⁠. The smoothed force is equal to the unsmoothed force until $𝜃s$⁠, at which point F_4,z goes continuously to zero at $𝜃234=180°$ before returning to the unsmoothed value at $𝜃s$ for $ϕ=−60°$⁠.

The chosen form of s is particularly computationally convenient because (1) s and $∇ns$ can be calculated from quantities that are already required for the unsmoothed potential and (2) no additional square-root evaluations are necessary if smoothing is applied symmetrically around $90°$ to eliminate both possible degeneracies (i.e., $𝜃c=180°$ also implies smoothing toward $𝜃c=0°$⁠). We encourage the use of the smoothed torsional potential especially during the development of coarse-grained models, using s as an additional adjustable parameter, in order to obtain correct simulation behavior and physics.

See supplementary material for an example implementation of the smoothed torsional potential.

This work was supported by the Princeton Center for Complex Materials (PCCM), a U.S. National Science Foundation Materials Research Science and Engineering Center (Award No. DMR-1420541). This research is part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (Award Nos. OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications.