Subdivision of the ``Bonded'' Potential

As we have seen in Sec. 2.3 the idea behind the multiple time step scheme is that of the reference system which propagates for a certain amount of time under the influence of some unperturbed reference Hamiltonian, and then undergoes an impulsive correction brought by the remainder of the potential. The exact trajectory spanned by the complete Hamiltonian is recovered by applying this impulsive correction onto the ``reference'' trajectory. We have also seen in the same section that, by subdividing the interaction potential, we can determine as many ``nested'' reference systems as we wish. The first step in defining a general protocol for the subdivision of the bonded potential for complex molecular systems consists in identifying the various time scales and their connection to the potential. The interaction bonded potential in almost all popular force fields is given as a function of the stretching, bending and torsion internal coordinates and has the general form

$$V_{\rm bnd} = V_{\rm stretch} + V_{\rm bend} + V_{\rm tors},$$ (4.2)

where

$$V_{\rm stretch} = \sum_{\rm Bonds} K_r \left ( r - r_0 \right )^2$$

$$V_{\rm bend} = \sum_{\rm Angles} K_{\theta} \left ( \theta - \theta_0 \right )^2.$$

$$V_{\rm tors} = \sum_{\rm Dihedrals} V_{\phi} \left [1 + \cos \left ( n \phi - \gamma \right)\right ].$$ (4.3)

Here, $K_r$ and $K_{\theta}$ are the bonded force constants associated with bond stretching and angles bending respectively, while $r_0$ and $\theta_0$ are their respective equilibrium values. In the torsional potential, $V_{tors}$, $\phi$ is the dihedral angle, while $K_{\phi}$, $n$ and $\gamma$ are constants.

The characteristic time scale of a particular internal degrees of freedom can be estimated assuming that this coordinate behaves like a harmonic oscillator, uncoupled form the rest the other internal degrees of freedom. Thus, the criterion for guiding the subdivision of the potential in Eq. (4.2) is given by the characteristic frequency of this uncoupled oscillator. We now give, for each type of degree of freedom, practical formula to evaluate the harmonic frequency from the force field constants given in Eq. (4.3).

Stretching: The stretching frequencies are given by the well known expression

$$\nu_{s} = {1 \over 2 \pi} \left ( {K_r \over \mu} \right )^{(1/2)},$$ (4.4)

where $\mu$ is reduced mass.

Bending: We shall assume for the sake of simplicity that the uncoupled bending frequencies depends on the masses of the atom 1 and 3 (see Fig. 4.1), that is mass 2 is assumed to be infinity. This turns out to be in general an excellent approximation for bending involving hydrogen and a good approximation for external bendings in large molecules involving masses of comparable magnitude. The frequency is obtained by writing the Lagrangian in polar coordinates for the mechanical system depicted in Fig. 4.1. The Cartesian coordinates are expressed in terms of the polar coordinates as

$$x_{1} = - r_{12} \sin(\alpha/2) ~~~~ y_{1} = r_{12} \cos(\alpha/2)$$ (4.5)

$$x_{3} = r_{32} \sin(\alpha/2) ~~~~ y_{3} = r_{32} \cos(\alpha/2)$$ (4.6)

where the distance $r_{32}$ and $r_{12}$ are constrained to the equilibrium values. The velocities are then

$$\dot x_{1} = -r_{12} \cos(\alpha/2) {\dot \alpha \over 2} ~~~~ \dot y_{1} = - r_{12} {\rm sin}(\alpha/2) {\dot \alpha \over 2}$$ (4.7)

$$\dot x_{3} = r_{32} \cos(\alpha/2) {\dot \alpha \over 2} ~~~~ \dot y_{3} = - r_{32} {\rm sin}(\alpha/2) {\dot \alpha \over 2}$$ (4.8)

The Lagrangian for the uncoupled bending is then

$${\cal L} = {1 \over 2} \left ( m_{1} \dot x_{1}^{2} + m_{1} \dot y_{1}^{2} + m_{3} \dot x_{3}^{2} + m_{3} \dot y_{3}^{2} \right )- V_{bend}$$ (4.9)

$$= {1 \over 8} \left ( m_{1}r_{12}^{2} + m_{2} r_{32}^{2}) \right )\dot \alpha^{2} - {1 \over 2} k_{\theta} (\alpha - \alpha_{0})^2.$$ (4.10)

The equation of motion ${d \over dt} {\partial {\cal L} \over \dot \alpha} - {\partial {\cal L} \over \alpha} =0$ for the $\alpha$ coordinate is given by

$$\ddot \alpha + {4 K_{\theta} \over I_{b}} (\alpha -\alpha_{0})^{2} = 0.$$ (4.11)

Where, $I_{b} = m_{1}r_{12}^{2} + m_{3} r_{32}^{2}$ is the moment of inertia about an axis passing by atom 3 and perpendicular to the bending plane. Finally, the uncoupled bending frequency is given by

$$\nu_{b} = {1 \over 2 \pi} \left ( {4 K_{\theta} \over m_{1}r_{12}^{2} + m_{3} r_{32}^{2}} \right )^{(1/2)}$$ (4.12)

Figure: Bending and dihedral angles

Torsion: We limit our analysis to a purely torsional system (see Fig. 4.1b ) where atoms 2 and 3 are held fixed, and all bond distances and the angle $\theta$ are constrained to their equilibrium values. The system has only one degree of freedom, the dihedral angle $\phi$ driven by the torsional potential $V_{\phi}$. Again we rewrite the kinetic energy in terms of the bond distances, the dihedral angle and the constant bend angle $\theta$. For the kinetic energy, the only relevant coordinates are now those of atoms 1 and 4:

$$x_{1} = - d_{12}\cos\theta + {d_{23} \over 2} ~~~~~~ x_{4} = d_{34}\cos\theta + {d_{23} \over 2}$$

$$y_{1} = d_{12}\sin\theta \cos(\phi/2) ~~~~~~ y_{4} = d_{34}\sin\theta \cos(\phi/2)$$

$$z_{1} = - d_{12}\sin\theta \sin(\phi/2) ~~~~~~ z_{4} = d_{34}\sin\theta \sin(\phi/2).$$ (4.13)

The Lagrangian in terms of the dihedral angle coordinate is then

$${\cal L} = {1 \over 8} I_{t} \dot \phi^{2} - V_{\phi} \left [1 + \cos \left ( n \phi - \gamma \right)\right ],$$ (4.14)

where

$$I_{t} = sin^{2}\theta \left ( m_{1} d_{12}^{2}+ m_{4} d_{34}^{2} \right ).$$ (4.15)

Assuming small oscillations, the potential may be approximated by a second order expansion around the corresponding equilibrium dihedral angle $\phi_{0}$

$$V_{tors} = {1 \over 2} \left ( {\partial ^{2} V_{tors} \over \partial \phi^{... ...i_{0}} (\phi - \phi_{0})^{2} = {1 \over 2} V_{\phi} n^{2} (\phi - \phi_{0})^{2}$$ (4.16)

Substituting (4.16) into Eq . (4.14) and then writing the Lagrange equation of motion for the coordinate $\phi$, one obtains again a differential equation of a harmonic oscillator, namely

$$\ddot \phi + {4 V_{\phi} n^{2} \over I_{t} } (\phi -\phi_{0}) = 0.$$ (4.17)

Thus, the uncoupled torsional frequency is given by

$$\nu_{t} = {n \over 2 \pi} \left ( 4 {V_{\phi} \over sin^{2}\theta \left ( m_{1} d_{12}^{2}+ m_{4} d_{34}^{2} \right )} \right )^{(1/2)}.$$ (4.18)

For many all-atom force fields, improper torsions [3,6] are modeled using a potential identical to that of the proper torsion in Eq. (4.3) and hence in these cases Eq. (4.18) applies also to the improper torsion uncoupled frequency, provided that indices 1 and 4 refer to the lighter atoms. In figure (4.2) we report the distribution of frequencies for the hydrated protein Bovine Pancreatin Trypsin Inhibitor (BPTI) using the AMBER [4] force field. The distributions might be thought as a density of the uncoupled intramolecular states of the system. As we can see in the figure there is a relevant degree of overlap for the various internal degrees of freedom. For example, ``slow'' degrees of freedom such as torsions may be found up to 600 wavenumber, well inside the ``bending'' region; these are usually improper or proper torsions involving hydrogen. It is then inappropriate to assign such ``fast'' torsions involving hydrogen to a slow reference system. We recall that in a multiple time simulation the integration of a supposedly slow degree of freedom with a excessively large time step is enough to undermine the entire simulation.

Figure: density of the uncoupled (see text) states for stretching, bending, proper and improper torsion obtained with the AMBER force field on the protein bovine pancreatic trypsin inhibitor (BPTI). Frequencies were calculated according to Eqs. (4.4,4.12,4.18)

In Fig. 4.2 we also notice that almost all the proper torsions fall below 350 cm$^{-1}$. An efficient and simple separation of the intramolecular AMBER potential [27] assigns all bendings stretching and the improper or proper torsions involving hydrogen to a ``fast'' reference system labeled $n0$ and all proper torsions to a slower reference system labeled $n1$. The subdivision is then

$$V_{\rm n0} = V_{\rm stretch} + V_{\rm bend} + V_{\rm i-tors} + V_{\rm p-tors}^{(h)}$$

$$V_{\rm n1} = V_{\rm\rm p-tors}$$ (4.19)

Where with $V_{\rm p-tors}^{(h)}$ we indicate proper torsions involving hydrogen. For the reference system $V_{\rm n0}$, the hydrogen stretching frequencies are the fastest motions and the $\Delta t_{\rm n0}$ time step must be set to 0.2-0.3 fs. The computational burden of this part of the potential is very limited, since it involves mostly two or three body forces. For the reference system $V_{\rm n1}$, the fastest motion is around 300 cm$^{-1}$ and the time step $\Delta t_{\rm n1}$ should be set to 1-1.5 fs. The computational effort for the reference system potential $V_{n1}$ is more important because of the numerous proper torsions of complex molecular systems which involve more expensive four body forces calculations. One may also notice that some of the bendings which were assigned to the n0 reference system fall in the torsion frequency region and could be therefore integrated with a time step much larger than $\Delta t_{\rm n0} \simeq$ 0.2-0.3. However, in a multiple time step integration, this overlap is just inefficient, but certainly not dangerous. Indeed, no instability may derive for integrating slow degrees of freedom with exceedingly small time steps.