Effective conservation of energy and momentum algorithm using switching potentials suitable for molecular dynamics simulation of thermodynamical systems

TL;DR

This paper introduces an energy and momentum conserving algorithm for molecular dynamics simulations that effectively manages energy violations during potential switching, ensuring stable thermodynamical properties.

Contribution

The novel algorithm maintains energy and momentum conservation during potential switching in MD simulations, improving accuracy over traditional methods.

Findings

The algorithm preserves system coordinates during moves.

It reduces energy violations during potential crossover.

Thermodynamical properties remain stable with the algorithm.

Abstract

During a crossover via a switching mechanism from one 2-body potential to another as might be applied in modeling (chemical) reactions in the vicinity of bond formation, energy violations would occur due to finite step size which determines the trajectory of the particles relative to the potential interactions of the unbonded state by numerical (e.g. Verlet) integration. This problem is overcome by an algorithm which preserves the coordinates of the system for each move, but corrects for energy discrepancies by ensuring both energy and momentum conservation in the dynamics. The algorithm is tested for a hysteresis loop reaction model with an without the implementation of the algorithm. The tests involve checking the rate of energy flow out of the MD simulation box; in the equilibrium state, no net rate of flows within experimental error should be observed. The temperature and pressure of the box should also be invariant within the range of fluctuation of these quantities. It is demonstrated that the algorithm satisfies these criteria.

Figures

Tables

Table 1. Table 1: Table with values for the mean heat supply per unit step and temperature. The error is one unit of standard error for the quantities.

Curve	${ϵ}_h$	${ϵ}_l$	Mean Temperature
l1	-.2274E+00 $\pm$0.19E-02	-.2295E+00 $\pm$0.21E-02	0.9063E+01 $\pm$0.62E-02
l2	-.5602E+00 $\pm$0.22E-02	-.5596E+00 $\pm$0.22E-02	0.1032E+02 $\pm$0.63E-02
l3	-.4161E-01 $\pm$0.14E-02	-.4089E-01 $\pm$0.14E-02	0.8774E+01 $\pm$0.79E-02
l4	-.5201E-01 $\pm$0.16E-02	-.5103E-01 $\pm$0.17E-02	0.8980E+01 $\pm$0.98E-02
t1	-.5312E-03 $\pm$0.92E-03	-.3334E-03 $\pm$0.76E-03	0.8082E+01 $\pm$0.49E-02
l5	0.1311E-02 $\pm$0.82E-03	0.1147E-02 $\pm$0.84E-03	0.7731E+01 $\pm$0.97E-02
t2	-.6823E-03 $\pm$0.12E-02	-.1507E-02 $\pm$0.13E-02	0.1216E+02 $\pm$0.17E-01
l6	0.7291E-02 $\pm$0.12E-02	0.6343E-02 $\pm$0.14E-02	0.1088E+02 $\pm$0.15E-01
t3	-.9348E-03 $\pm$0.18E-02	-.3379E-02 $\pm$0.17E-02	0.1622E+02 $\pm$0.18E-01
l7	0.1918E-01 $\pm$0.14E-02	0.1938E-01 $\pm$0.16E-02	0.1329E+02 $\pm$0.20E-01

Equations

$$2\text{A}\begin{array}[]{*{20}c}{k_{1}}\\ \rightleftarrows\\ {k_{-1}}\\ \end{array}\text{A}_{\text{2}}$$

$$u(r)=u_{vib}(r)s(r)+u_{LJ}\left[1-s(r)\right]$$

$$u_{vib}(r)=u_0+\frac{1}{2}k(r-r_0{)}^2$$

$$s(r)=\frac{1}{1+{\left(\frac{r}{r_{sw}}\right)}^n}$$

$$\left\{\begin{array}[]{ll}s(r)&\rightarrow 1\text{ \ \ \ \ \ if }r<r_{sw}\\ s(r)&\rightarrow 0\text{ \ \ \ \ \ for }r>r_{sw}\end{array}\right..$$

$$⟨\sum _{i=1}^N{\mathbf{p}}_{\mathbf{i}}.{\mathbf{p}}_{\mathbf{i}}/m_i⟩=3N{k}_{B}T\text{and}P=\rho k_{B}T+W/V$$

$$\delta W_2+\delta W_1-({\varphi }_b(r_f)-{\varphi }_b(r_c))=\Delta \sum (K.E.)$$

$$\delta W_2+\delta W_1-(V^{vir}(r_f)-V^{vir}(r_c))=\Delta \sum {}^{\prime }(K.E.)$$

$$\delta V^{vir}(r)-\delta {\varphi }_b(r)=\delta \left(\Delta \sum (K.E{.}_{x,y,z})-\Delta \sum {}^{\prime }(K.E.{)}_{x,y,z}\right)$$

$${\varphi }_a(r_f)=V^{vir}(r_f);{\varphi }_a(r_c)=V^{vir}(r_c);{\varphi }_b(r_c)=V^{vir}(r_c);$$

$$\Delta ={\varphi }_b(r_f)-V^{vir}(r_f)$$

$$=$$

$$\mathbf{v_{i}\,^{\prime}}=(1+\alpha)\mathbf{v_{i}}+{\mbox{\boldmath$\beta$}}$$

$${{\mbox{\boldmath$\beta$}}}=-\alpha\sum m_{i}\mathbf{v_{i}}/\sum m_{i}$$

$$Q={\left(\sum m_i{\mathbf{v}}_{\mathbf{i}}\right)}^2/{\left(\sum m_i\right)}^2$$

$$\mathbf{v^{\prime}_{i}}^{2}=(1+\alpha)^{2}\mathbf{v_{i}}^{2}+2(1+\alpha)\mathbf{v_{i}}.{\mbox{\boldmath$\beta$}}+{\mbox{\boldmath$\beta$}}^{2}.$$

$$\sum \frac{1}{2}{m}_i{{\mathbf{v}}_{\mathbf{i}}^{\prime }}^2-\sum \frac{1}{2}{m}_i{{\mathbf{v}}_{\mathbf{i}}}^2=\Delta$$

$$\Delta =\frac{{\alpha }^2{m}_1{m}_2}{2(m_1+m_2)}\left[{{\mathbf{v}}_1}^2+{{\mathbf{v}}_2}^2-2{\mathbf{v}}_1.{\mathbf{v}}_2\right]$$

$$+\frac{2\alpha m_1{m}_2}{2(m_1+m_2)}\left[{{\mathbf{v}}_1}^2+{{\mathbf{v}}_2}^2-2{\mathbf{v}}_1.{\mathbf{v}}_2\right].$$

$${\alpha }^{2}qa+2qa\alpha -\Delta =0$$

Full text

AMS (MSC2000) Subject Classification. 00A71-2, 70H05, 80A20

1 PRELIMINARIES

The dimeric particle reaction simulated may be written

[TABLE]

where $k_{1}$ is the forward rate constant and $k_{-1}$ is the backward rate constant. The reaction simulation was conducted at extremely high temperatures which are off-scale and not used in ordinary simulations of $LJ$ (Lennard-Jones) fluids where normally [1] the reduced temperatures $T^{\ast}$ ranges $\sim 0.3-1.2$, whereas here, $T^{\ast}\sim 8.0-16.0$, well above the supercritical regime of the $LJ$ fluid At these temperatures, the normal choices for time step increments do not obtain without also taking into account energy-momentum conservation algorithms in regions where there are abrupt changes of gradient [1, 2, 3]. The global literature does not seem to cover such extreme conditions of simulation with discrete time steps using the Verlet velocity algorithm. The units used here are reduced $LJ$ ones [1]. The simulation was at density $\rho=0.70$ with $4096$ atomic particls which could react. The potentials used are as given in Fig. (1)   where $r_{b}=1.20$ for the vicinity where the bond of the dimer is broken and where 2 free particles emerge, and $r_{f}=0.85$ is the point along the hysteresis potential curve where the dimer is defined to exist for two previously free particles which collide. The reaction proceeds as follows; all particles interact with the splined $LJ$ pair potential $u_{LJ}$ except for the dimeric pair $(i,j)$ formed from particles $i$ and $j$ which interact with a harmonic-like intermolecular potential modified by a switch $u(r)$ given by

[TABLE]

where $u_{vib}(r)$ is the vibrational potential given by eq.(1.3) below

[TABLE]

The switching function $s(r)$ is defined as

[TABLE]

where

[TABLE]

The switching function becomes effective when the distance between the atoms approach the value $r_{sw}$ (see Fig. (1) ). Some of the other parameters used in the equations that follow include:

$u_{0}=-10,r_{0}=1.0,k\sim 2446$ (exact value is determined by the other input parameters), $n=100,r_{f}=0.85,r_{b}=1.20,\mbox{and}\;r_{sw}=1.11$. Particles $i$ and $j$ above also interact with all other particles not bonded to it via $u_{LJ}$ . Full simulation details are given elsewhere [2]; suffice to say the activation energy at $r_{f}$ is extremely high at approximately $17.5$. At $r_{f}$, the molecular potential is turned on where at this point there is actually a crossing of the potential curves although the gradients of the molecular and free $u_{LJ}$ potentials are ”‘very close”’. On the other hand, at $r_{b}$ , the switch forces the two curves to coalesce, but detailed examination shows that there is an energy gap of about the same magnitude as the cut-off point in a normal non-splined $LJ$ potential ($\sim 0.04$ energy units), meaning there is no crossing of the potentials. The current algorithm is applied for both these cross-over regions with their different mechanisms of cross-over. The MD cell is rectangular, with unit distance along the axis ( $x$ direction) of the cell length, whereas the breadth and height was both $1/16$, implying a thin pencil-like system where the thermostats were placed at the ends of the MD cell, and the energy supplied per unit time step $\delta t$ at both ends of the cell (orthogonal to the $x$ axis) in the vicinity of $x=0$ and $x=1$ maintained at temperatures $T_{h}$ and $T_{l}$ could be monitored, where this energy per unit step time is respectively $\epsilon_{h}$ and $\epsilon_{l}$. At equilibrium, (when $T_{h}=T_{l}$), the net energy supplied within statistical error (meaning 1-3 units of the standard error of the $\epsilon$ distributions ) is zero, i.e. $\epsilon_{l}\approx\epsilon_{h}\approx 0$. The cell is divided up uniformly into 64 rectangular regions along the $x$ axis and its thermodynamical properties of temperature and pressure are probed. The resulting values of the $\epsilon$’s and the relative invariance of the pressure and temperature profiles would be a measure of the accuracy of the algorithm from a thermodynamical point of view at the steady state. For systems with a large number of particles such thermodynamical criteria are appropriate. The synthetic thermostats now frequently used in conjunction with ”‘non-Hamiltonian”’ MD [3] cannot be employed for this type of study, where actual energy increments are sampled. The runs were for 4 million time steps, with averages taken over 100 dumps, where each variable is sampled every 20 time steps. The final averages were over the 20-100 dump values of averaged quantities.

The temperature $T$ and pressure $p$ are computed by the equipartition and Virial expression given respectively by

[TABLE]

where $W=-\frac{1}{3}\sum_{i}\sum_{j>i}w(r_{ij})$ and the intermolecular pair Virial $w(r)$ is given by $w(r)=r\frac{dv(r)}{dr}$ with $v$ being the potential.

2 ALGORITHM AND AND ANALYSIS OF NUMERICAL RESULTS

The velocity Verlet algorithm [4, p. 81]used here [1] and allied types generate a trajectory at time $n\delta t$ from that at $(n-1)\delta t$ with step increment $\delta t$ through a mapping $T_{m}$ where $(v(n\delta t),r(n\delta t))=T_{m}(v((n-1)\delta t),r((n-1)\delta t))$ which does not scale linearly with $\delta t$. For a Hamiltonian $H$ whose potential $V$ is dependent only on position $r$ having momentum components $p_{i}$, the system without external perturbation has constant energy $E$, and the normal assumption in MD (NAMD)is that for the $n^{th}$ step, $\Delta E_{n}=\left|H(n\delta t)-E\right|\leq\epsilon$ and also $\sum_{i=1}^{N}\Delta E_{i}\leq\epsilon^{s}$ for the specified $\epsilon^{\prime}s$. In the simulation under NAMD, the force fields are constant and do not change for any one time step. In these cases, the energy is a constant of the motion for any time interval $\delta t_{T}$ when no external perturbations (e.g. due to thermostat interference) are impressed. When there is a crossing of potentials at such a time interval interval from $\phi_{b}$ to $\phi_{a}$ at an inter particle distance($icd$ ) $r_{c}$ (such as points $r_{f}$ and $r_{b}$ of Fig. (1)) of general particle $1$ and $2$ (the $(1,2)$ particle pair) due to a reactive process (such as occurs in either direction of (1.1)) a bifurcation occurs where the MD program computes the next step coordinates as for the unreacted system (potential $\phi_{b}$), which needs to be corrected. Let the $icd$ at time step $i$ be $r_{i}$ (with $\phi_{b}$ potential) and at step $i+1$ after interval $\delta t$ be $r_{f}=r_{i+1}$ where $r_{f}<r_{c}<r_{i}$. Due to this crossover, a different Hamiltonian $H^{\prime}$ is operative after point $r_{c}$ is crossed, where under NAMD, the other coordinates not undergoing crossover are not affected. For what follows, subscripts refer to the particle concerned. Let the interparticle potential at $r_{f}$ be $E_{a}=E_{f}=\phi_{a}(r_{f})$ and at $r_{f}$ be $E_{b}=\phi_{b}(r_{f})$, where $\Delta=E_{b}-E_{a}$. Then if $r_{f}$ be the final coordinate due to the $\phi_{b}$ potential and force field, two questions may be asked: (i) Can the velocities of (1,2) be scaled, so that there is no energy or momentum violation during the crossover based on the $\phi_{b}$ trajectory calculation and (ii) Can a pseudo stochastic potential be imposed from coordinates $r_{c}$ (at virtual time $t_{c}$) to $r_{f}$ such that (i) above is true? For (ii) we have

Theorem 2.1

A virtual potential which scales velocities to preserve momentum and energy can be constructed about region $r_{c}.$

proof The external work done $\delta W$ on particles 1 and 2 over the time step is proportional to the distance traveled since these forces are constant and so for each of these particles $i$, $\mathbf{F_{ext,i}.\Delta r_{i}}=\delta W_{i}$ where $\mathbf{\Delta r_{i}}$ is the distance increment during at least part of the time step from $r_{c}$ to $r_{f}$. For the non-reacting trajectory over time $\lambda\delta t$ ($\lambda\leq 1)$ (virtual because it is not the correct path due to the crossover at $r_{c}$),

[TABLE]

where $\Delta\sum(K.E.)$ is the change of kinetic energy for the $(1,2)$ pair from the First Law between the end points $r_{f},\,r_{c}$. Now over time interval $t_{c}$ to $t_{f}$, for the reactive trajectory, we introduce a ”‘virtual potential”’ $V^{vir}$ that will lead to the same positional coordinates for the pair at the end of the time step with different velocities than for the non-reactive transition leading to the transition

[TABLE]

where $\Delta\sum\,^{\prime}(K.E.)$ is the change of kinetic energy for the pair with $V^{vir}$ turned on and along this trajectory, the change of potential for $V^{vir}$ is equated to the change in the K.E. of the pair as given in the results of theorem (2.2) for all three orthogonal coordinates, i.e.

[TABLE]

with momentum conservation, that is $\delta V^{vir}(r_{i})=\delta\phi_{a}(r_{i})$ for the variation along the $r_{i}$ coordinate, but $\delta\phi_{a}(r_{i})=-\delta K.E.$ along internuclear coordinate $r_{i}$ whereas $\delta V^{vir}=-K.E.$ (scaled about all three axes). Continuity of potential implies

[TABLE]

Subtracting (2.1) from (2.2) and applying b.c.’s (2.6) leads to

[TABLE]

The above shows that a conservative virtual potential could be said to be operating in the vicinity of the transition (from $t_{c}\,\text{to}\,t_{a}$) .$\bullet$

Question (i) above leads to:

Theorem 2.2

Relative to the velocities at any $r_{f}$ due to the $\phi_{b}$ potential, the rescaled velocities $\mathbf{v\,^{\prime}}$ due to the potential difference $\Delta$ leading to these final velocities due to the virtual potential can have a form given by

[TABLE]

(where $i=1,2$) for a vector $\beta$.

proof From the $\mathbf{v}$ velocities at $r_{f}$ due to $\phi_{b}$ we compute the $\mathbf{v^{\prime}}$ velocities at $r_{f}$ due to the virtual potential. Since net change of momentum is due to the external forces only, which is invariant for the $(1,2)$ pair, conservation of total momentum relating $\mathbf{v^{\prime}}$ and $\mathbf{v}$ in (2.6) yields a definition of $\beta$ ( summation from 1 to 2 for what follows, where the mass of particle $i$ is $m_{i}$)

[TABLE]

Defining for any vector \mathbf{s}^{2}=\mathbf{s.s},{\mbox{\boldmath\beta}}^{2}=\alpha^{2}Q, where

[TABLE]

then the rescaled velocities become from (2.6)

[TABLE]

With $\Delta=E_{b}-E_{a}$, Energy conservation implies

[TABLE]

The coupling of (2.9-2.10) leads, after several steps of algebra to

[TABLE]

Defining $a=(\mathbf{v_{1}}-\mathbf{v_{2}})^{2}$,   $q=m_{2}m_{1}/[{2(m_{1}+m_{2})}]$,  ($q>0,a\geq 0)$, then the above is equivalent to the quadratic equation

[TABLE]

and in simulations, only $\alpha$ is unknown and can be determined from (2.12) where real solutions exist for $\Delta/qa\geq-1$. $\bullet$ The above Inequality leads to a certain asymmetry concerning forward and backward reactions, even for reversible reactions where the region of formation and breakdown of molecules are located in the same region with the reversal of the sign of approximate $\Delta$. For this simulation, a reaction in either direction (formation or breakdown of dimer ) proceeds if (LABEL:e16) is true; if not then the trajectory follows the one for the initial trajectory without any reaction (i.e. no potential crossover).

**Interpretation of results.  **Fig. (1) shows a rapidly changing potential curve with several inflexion points used in the simulation at very high temperature (as far as I know such ranges have not been reported in the literature for non-synthetic methods) warranting smaller time steps; larger ones would introduce errors due to the rapidly changing potential and high K.E.; thus, even with the application of the algorithm between cordinates $r_{f}$ and $r_{b}$, curves l1 and l2 have too large a $\delta t$ value to achieve equilibrium - meaning flat or invariant - temperature (see Fig. (2) ) or pressure (see Fig. (3))or unit step thermostat heat supply (see Table  1)($\epsilon_{h}$ and $\epsilon_{l}$) profiles where for these curves, the $(\epsilon_{h},\epsilon_{l})$ values show net heat absorption; the curve at t1 (with $\delta t=5.0\,ep-5$ show flat profiles (within statistical fluctuations and 2 standard errors of variation) for temperature, pressure and net zero heat supply; and this choice of time step interval was found adequate for runs at much higher temperatures ($T=12$ and $T=16$) which was used to determine thermodynamical properties [2]. For this $\delta t$ value and all others, no reasonable stationary equilibrium conditions could be obtained without the application of the algorithm (curves l2,l4,l5,l6 and l7). The algorithm is seen to be effective over a wide temperature range for this complex dimer reaction simulated under extreme values of thermodynamical variables and the results here do not vary for longer runs and greater sampling statistics (e.g. 6 or 10 million time steps). The thin, pencil-like geometry of the rectangular cell with thermostats located at the ends would highlight the energy non-conservation leading to a non-flat temperature distribution, as observed and which was used to determine the regime of validity of the algorithm.