Subdivision the Non Bonded Potential

In addition to the long range electrostatic contributions, $V_{\rm qr}$ and $V_{\rm qd}$, given in Eqs. (4.20,4.21), more short range forces play a significant role in the total non bonded potential energy. The latter can be written as:

$$V_{\rm nbn} = V_{\rm vdw} + V_{\rm qr} + V_{\rm qd} + V_{\rm 14}.$$ (4.36)

Where, $V_{\rm vdw}$ is the Lennard-Jones potential, namely

$$V_{\rm vdw} = {\sum_{i < j}^N}^{\prime} 4 \epsilon_{ij} \left [ \... ...j} } \right)^{12} - \left ( { \sigma_{ij} \over r_{ij} } \right )^{6} \right ].$$ (4.37)

Here, the prime on the sum indicates that interactions between atoms separated by less than three consecutive bonds must be omitted. The term $V_{\rm 14}$ is typical for force fields of complex molecular systems [3,5]. While non bonded forces between atoms involved in the same covalent bond or angle bending interaction are generally excluded, the potential between atoms separated by three covalent bonds is retained and readjusted in various ways. In all cases, the $V_{\rm 14}$ term remains in general a very stiff and, hence, a fast varying term. The computational cost of the $V_{\rm 14}$ contribution is very small compared to other non bonded interactions. Thus, it is safer to assigns this potential term to the slowest intramolecular reference system potential $V_{\rm n1}$ of Eq. (4.19).

The $V_{\rm qr}$ reciprocal lattice term, including the correction due to the excluded or partially excluded (i.e. the electrostatic part of $V_{\rm 14}$) interactions cannot be split when using SPME and must be assigned altogether to only one reference system. The time scale of the potential $V_{\rm qr}$ depends on the convergence parameter $\alpha$. Indeed, this constant controls the relative weights of the reciprocal lattice energy $V_{\rm qr}$, and of the direct lattice energy $V_{\rm qd}$. By increasing $\alpha$, one increases the weight of the reciprocal lattice contribution $V_{\rm qr}$ to the total Coulomb energy. When using SPME the cost of the reciprocal lattice sums is cut down dramatically and, therefore, the use of large $\alpha$'s becomes helpful to reduce the computational burden of the direct lattice calculation. For a value of $\alpha$ increased beyond a certain limit, there is no longer a computational gain, since the pair distances must always be evaluated in direct space until convergence of the Lennard-Jones energy (usually occurring at a 10 Å  cutoff). Furthermore, the larger is $\alpha$, the more short-ranged and fast varying becomes the potential $V_{\rm qr}$, thus requiring short time steps to integrate correctly the equations of motion. A good compromise for $\alpha$, valid for cell of any shape and size, is $\alpha =$0.4-0.5.

The direct space potential is separated [13,27] in three contributions according to the interaction distance. The overall non bonded potential breakup is therefore

$$V_{\rm n1} = V_{\rm 14}$$

$$V_{m} = V_{\rm vdw}^{(1)} + V_{\rm qd}^{(1)}$$

$$V_{l} = V_{\rm vdw}^{(2)} + V_{\rm qd}^{(2)} + V_{\rm qr}$$

$$V_{h} = V_{\rm vdw}^{(3)} + V_{\rm qd}^{(3)}$$ (4.38)

where the superscripts $m,l,h$ of the direct space term $V_{\rm vdw}$ and $V_{\rm qd}$ refer to the short, medium and long range non-bonded interactions, respectively. The $m$-th reference system includes non-bonded direct space interactions at short range, typically between 0 to $4.3$-$5.3$ Å. $V_{l}$ contains both the medium range direct space potential, with a typical range of $4.3$-$5.3$ to $7.3$-$8.5$ Å, and the reciprocal space term, $V_{\rm qr}$. Finally, the $h$-th reference system, which is the most slowly varying contains, the remaining direct space interactions from $7.3$-$8.3$ Å to cutoff distance. As the simulations proceeds the particles seen by a target particle may cross from one region to an other, while the number of two body contacts in one distance class [19] or reference system potential must be continuously updated. Instabilities caused by this flow across potential shell boundaries are generally handled by multiplying the pair potential by a group-based switching function [24]^4.3. Thus, at any distance $r$ the direct space potential $V$ can be written schematically as:

$$V = V_{1} + V_{2} + V_{3}$$ (4.39)

with

$$V_1 = V S_1$$ (4.40)

$$V_2 = V (S_2 - S_1)$$ (4.41)

$$V_3 = V (S_3 - S_2)$$ (4.42)

where $S_j$ is the switching function for the three shells, $j=m,l,h$ defined as:

$$S_{j} (R) = \left \{ \begin{array}{ll} 1 & R_{j-1} \le R < R_{j} ... ...e R < R_{j}+\lambda_{j} \\ 0 & R_{j}+\lambda_{j} < R . \\ \end{array} \right \}$$ (4.43)

Table: Potential breakup and relative time steps for complex systems with interactions modeled by the AMBER[4] force field and electrostatic computed using the SPME method

Component	Contributions	Spherical Shells	Time step
$V_{\rm n0}$	$V_{\rm stretch}+V_{\rm bend}+$	-	$\Delta t_{\rm n0}=0.33~fs$
	$+ V_{\rm i-tors}+V_{\rm p-tors}^{(h)}$
$V_{\rm n1}$	$V_{\rm p-tors} + V_{\rm 14}$	-	$\Delta t_{\rm n1}=1.0~fs$
$V_{\rm m}$	$V_{\rm LJ}^{(1)} + V_{\rm qd,\alpha = 0.43}^{(1)}$	$0 < r < 4.5~\AA$	$\Delta t_{m}=2.0~fs$
$V_{l}$	$V_{\rm LJ}^{(2)} + V_{{\rm qd},\alpha = 0.43}^{(2)} +$	$4.5 \le r < 7.5 ~\AA$	$\Delta t_{l}=4.0~fs$
	$+ V_{{\rm qr},\alpha = 0.43}$
$V_{h}$	$V_{\rm LJ}^{(3)} + V_{\rm qd,\alpha = 0.43}^{(3)}$	$7.5 \le r < 10.0 ~\AA$	$\Delta t_{h}=12.0~fs$

Here, $R$ is the intergroup distance and $\lambda_j$ is the healing interval for the $j$-th shell. While $R_{0}$ is zero, $R_{1}=R_{m}$, $R_{2}=R_{l}$, and $R_{3}=R_{h}$ are the short, medium, long range shell radius, respectively. The switching $S^{(j)}_{3p}(R)$ is 1 at $R_{j}$ and goes monotonically to 0 at $R_{j}+\lambda_{j}$. Provided that $S_{3p}^{(j)}$ and its derivatives are continuous at $R_{j}$ and $R_{j}+\lambda_{j}$, the analytical form of $S_{3p}$ in the healing interval is arbitrary [16,20,24,13]. The full breakup for an AMBER type force field along with the integration time steps, valid for any complex molecular system with strong electrostatic interactions, is summarize in table II. The corresponding five time steps integration algorithm for the NVE ensemble is given by

$$\begin{array}{lll} e^{iL\Delta t_{h}} & = \exp \left [- { \pa... ...p}_{i}} {\Delta t_{h} \over 2} \right ], \right .\\ \end{array}$$

where $n_{l} = \Delta t_{h} /\Delta t_{l}$, $n_{m} = \Delta t_{l} /\Delta t_{m}$, $n_{n1} = \Delta t_{m} /\Delta t_{n1}$, $n_{n0} = \Delta t_{n1} /\Delta t_{n0}$. The explicit integration algorithm can be easily derived applying the five-fold discrete time propagator (4.46) to the state vector $\{ {\bf p}, {\bf q} \}$ at time 0 using the rule Eq. (2.23). The efficiency and accuracy for energy conservation of this r-RESPA symplectic and reversible integrator have been discussed extensively in Refs. [13,2]. Extension of this subdivision to non NVE simulation is described in Ref. [27].