Simulation Algorithm That Conserves Energy and Momentum for Molecular Dynamics of Systems Driven by Switching Potentials

Whenever there exists a crossover from one potential to another, computational problems are introduced in Molecular Dynamics MD simulation. These problem are overcome here by an algorithm, described in detail. The algorithm is applied to a 2-body particle potential for a hysteresis loop reaction model. Extreme temperature conditions were applied to test for algorithm effectiveness by monitoring global energy, pressure and temperature discrepancies in an equilibrium system. No net rate of energy and other flows within experimental error should be observed, in addition to invariance of temperature and pressure along the MD cell for the said system. It is found that all these conditions are met only when the algorithm is applied. It is concluded that the method can easily be extended to Nonequilibrium MD NEMD simulations and to reactive systems with reversible, non-hysteresis loops.


Introduction
The packages used extensively for biophysical simulations include CHARMM, GROMACS, DL POLY, IMD, and AMBER 1-5 where routine model reactions are not included.This is a significant limitation at attempting to model processes.One of the many reasons is that in the space of potential interactions, it is not so easy to specify when a species comes into existence, and when it ceases to exist; these criteria at the molecular level seems to be arbitrary or subjective, and macroscopic level designations are not always applicable.Another reason is that the complex potentials are not single-valued and would require costly 3 body interaction calculations.Probably newer phenomena-from the simulation point of view-might be uncovered if cost-effective reactive potentials could be used.The current algorithm is primitive enough to provide first steps in this direction for very large molecular 2 Mathematical Problems in Engineering assemblies encountered in Biophysical simulation where n-body n > 2 potentials would be currently prohibitive in terms of computational time.This work presents a simple dimeric single bond reaction as an example of how reactions can be included.From these tests, it is concluded that the method can easily be extended to reversible chemical systems.The theory is general and works for n-body interactions for the coordinates involved in any crossover trajectory irrespective of the degree of interaction.Recently, a 3-body potential modeled after the method of Stillinger 6 was used to study a chemical reaction at relatively low temperatures where by assumption, no precautions were taken into account for changes due to the steep gradients in the potential.While such methods might obtain for non-synthetic MD when the temperature and other thermodynamical variables are relatively small in magnitude, implying low velocities of the particles, more care must be taken at extremely high temperatures.Synthetic MD methods-meaning those techniques where the equations of motion are solved so as to replicate the ensemble probability distribution of the specified Hamiltonian and where the temperature, pressure, and other thermodynamical variables are introduced into a pseudo-Hamiltonian directly so that the successive trajectory coordinates can be computed at any time, with fixed temperature or pressure variables, 7-10 -unlike non synthetic methods where thermostats and barostats are placed in possibly localized regions by perturbing the system, are probably more immune to energy violations due to the artificial nature of computing the particle trajectories.We note that the gerund forms such as "thermostatting" for thermostat has been used routinely to refer to the application of a method to a simulation to control temperature 11, pages 143-144, 535 .The same applies for the gerund forms of barostats 11, pages 158-160 which refer to a method for controlling the pressure of a system under simulation.However, even for synthetic simulations, conservation of energy, and momentum at crossovers could lead to less fluctuation of the quantities or variables which appear and are set in the apparent or synthetic Hamiltonian.Thus various dispersion and variance relationships, from which quantities like the Specific Heat are derived, would differ from those described by a nonsynthetic Hamiltonian.Elementary statistical mechanics and the Boltzmann H theorem all indicate that in systems with a conservative Hamiltonian, the equilibrium state would be that in which equipartition of energy and the Maxwell distribution of velocities would obtain, and this would always be the case as the system relaxes to equilibrium; hence the temperature T too, as determined by the equipartition result n 3/2 kT n i 1 1/2 m i v 2 i .For large enough systems in terms of number of free coordinates, these results imply that once a system has reached equilibrium, it would persist in that state indefinitely.If, therefore, the simulation algorithm were perfect, it would evolve about an equilibrium trajectory indefinitely.There would, therefore, be no need to "thermostat" any closed system.For canonical and other ensembles, on the other hand, where energy exchange is investigated, then thermostatting would be required even for perfect algorithms that produce for each discrete time step the exact trajectory as would be produced in principle by integrating exactly the equations of motion.For the majority of cases however, thermostats are used just to maintain a particular temperature for a system at either local or global equilibrium.Hence, one can conclude for these cases that thermostats are implemented as a corrective to imperfect trajectory algorithms that does not constantly place the system on an equilibrium trajectory, as required by the H theorem; that is, it would appear that in the majority of instances, thermostatting in equilibrium systems is related to the fact that one's algorithm is imperfect in the sense that it cannot compute a bona fide equilibrium state.This can only be the case if accumulated machine and computational errors create a trajectory not corresponding to a microcanonical or canonical ensemble of states.Hence, for these situations, thermostatting refers to the implementation of algorithms that forces the system to adopt on average a canonical or microcanonical distribution of energies among the principle components within the system.In synthetic methods, the "actual" trajectory is not traced, but one that reproduces a canonical trajectory, but even here, opinions differ as to how accurately these trajectories are traced.Indeed, recent work seems to show that external perturbations can modify the "noise" spectrum of a natural system.For instance, the presence of an external random contribution to a high-frequency periodic electric field can reduce the total noise power 12 .This suggests that some natural properties connected to time correlation functions is a function of external perturbations and so one may conclude that basic synthetic methods may not include such elements of stochasticity.Another interesting observation of 12, 13 is the use of Monte Carlo techniques to model the system.In this case, Monte Carlo is used to simulate the dynamics of electrons in the semiconductor lattice by taking into account stochastic averaging.This is to be contrasted with the method here of attempting direct and approximate integration of the equations of motion, moderated by probabilistic inputs of energy at the ends of the box to simulate a "thermostat."One guess is that such Monte Carlo methods might be suitable if the details of molecular motion are not being investigated, and that given that a particular form of behavior is accepted, then one might superimpose stochasticity upon it through a Monte Carlo algorithm to simulate scattering phenomena, which includes temperature control.One possible problem with synthetic methods is that if a phenomenon is due to the system being in a particular phase space of a particular fixed Hamiltonian, then such events may not be detected or may be underrepresented in these synthetic methods.An overview of some of the above is in order.In the Nosé-Hoover method, one defines a Lagrangian for the system coordinates {ṙ i , ṗi } as where β is the temperature parameter.This so-called Lagrangian defines the conjugate momenta to r i and s as, respectively, p i m i s 2 ṙi and p s Q ṡ.Then for this system, there results ultimately a pseudo-Hamiltonian: whose trajectory is determined by the coupled equations that must be solved: where the last equation in superfluous.Equation 1.3 is solved by special and time consuming techniques that are not typical of those used for the standard Hamiltonian, such as the well-known Verlet and Gear algorithms.An analogous set of equations can be derived for constant pressure studies 14 .Another algorithm to correct for machine errors in following a PES is temperature-coupling method of Berendsen et al. 15 which has been widely used in many systems but it is claimed 11, page 161 that the canonical distribution is not produced "exactly."In this method, the velocities are scaled every time step by factor λ given by The upshot of the above is that these algorithms can be viewed as some sort of corrective procedure used to overcome problems of trajectory calculation accuracy for the rather simplistic, single-valued potential that are used for nonreactive systems due to the nonperfect integration of the equations of motion 16 .Paradoxically perhaps, the theories that were developed never allude to the machine error basis behind equilibrium thermostatting, which is not required by the H theorem when the system relaxes to equilibrium, and thus hardly any reference is made to the error in the computations of their new equations of motion that incorporates fixed thermodynamical variables like the pressure or temperature.It may be argued that they were referring to the canonical ensemble, but a careful examination of the Nosé justification of the method refers to a microcanonical phase space trajectory with the delta function having a component form δ H − E α .This might imply that machine error was not the foremost reason for the invention of the algorithms together with the accompanying theory.To date however, there has been little-if any-development in providing corrective measures to trajectory calculations for multivalued and other potentials which require switches to transfer trajectories from one PES to another for various molecular species which involve the variables pertaining to the surfaces.This particular review refers to one such attempt, which will be described in detail in what follows.
For both synthetic and nonsynthetic methods using n-body potentials, various switches and lists would have to be created to keep track of which potential energy surface a particle can transit to, and exit from, in order to define when a molecule is formed or destroyed in a reaction.Nonsynthetic Nonequilibrium Molecular Dynamics NEMD does not presuppose a theory concerning molecular interactions and therefore if new phenomena and relationships are sought in simulation studies, making use of quantitative values for the mechanical variables, conservation algorithms would have to be employed for systems with multipotential surfaces.In such studies, especially under extreme conditions, algorithms that can control energy and momentum variations so that larger time steps could be utilized seem essential; for nonsynthetic methods, they would be essential because gradients of energy flow could be artificially induced by violation of energy and momentum conservation due to the extreme potential gradients, thereby compromising any quantitative studies in nonequilibrium energy flows in NEMD simulations where gradients of thermodynamical variables exists by imposed boundary conditions.To initiate such studies, especially at extreme conditions, an algorithm was devised to correct for such momentum and energy conservation violations at crossover points in the potential curves due to reactions.The method is applicable to any n-body interaction system; in our case, we use a 2-body interaction system with switches that can turn on the potentials at prescribed distances.The model reaction simulated may be written as where k 1 and k −1 are the forward and backward rate constant, respectively.The reaction simulation was conducted at high temperatures not used ordinarily in simulations of LJ Lennard-Jones fluids where the reduced temperatures T * all units used are reduced LJ ones 17 ranges ∼0.3-1.2, 17 whereas here, T * ∼ 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 11, 17, 18 .The global literature does not seem to cover such extreme conditions of simulation with these precautions.The simulation was at density ρ 0.70 with 4096 atomic particles.The potentials used are as given in Figure 1 where r b 1.20 for the vicinity where the bond of the dimer is broken 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.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 u r u vib r s r u LJ 1 − s r , 1.6 where u vib r is the vibrational potential given by 1.7 The switching function s r is defined as where 1.9 The switching function becomes effective when the distance between the atoms approach the value r sw see Figure 1 .Some of the other parameters used in the equations that follow include u 0 −10, r 0 1.0, k ∼ 2446 exact value is determined by the other input parameters , n 100, r f 0.85, r b 1.20, and r sw 1.11.Switches are commonly encountered in theoretical accounts of complex interactions, such as found in polymer interactions and in chemical reactions.There are many flavors of switch categories, and some are more effective than others in forcing the merging of one potential type to another for a given distance defined by a metric 19-23 .The ideal switch would resemble a Heaviside step function but such functions cannot be so easily incorporated into the dynamical equations of motion which feature continuous variables because the various orders of differentials must be defined and computable over the discrete time steps.For instance, a switching function with several known applications, including those from statistical mechanics is given by the form 19 : where a and b are defined constants and the R s represent distances.For various optimization schemes to check for global minima, such as claimed in the Hunjan-Ramaswamy global optimization method, switches such as the g t function of the following form has been used 20 : where t is a time-dependent variable.On the other hand, for cluster dynamics, a switch of the form 21, equation 7 , is used, where the parameters E and F are adjusted to minimum energy of sub-clusters according to their species partitioning scheme.Switches without explicit details have been mentioned in other complex molecular structural studies to define topologies 22 .Similarly,

1.13
The above all refer to clearly defined spatial boundaries where there is a change of potential interaction type.In stochastic analysis 24, page 4 of response functions to a symmetric dichotomous switch variable ξ t having values ±1, analytical values may be derived for the flow variables.The situation here, on the other hand, is not stochastic where the particle trajectories are concerned, and the rate laws can only be determined through simulation and integration of the system Hamiltonian with the system potential given by 1.6 .Coming back to our description of the particle dynamics and the switch function, our particles i and j above also interact with all other particles not bonded to it via u LJ .Details of these potentials and their interactions are given elsewhere 18 ; here we note the high activation energy at r f of 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 nonsplined LJ potential ∼0.04 energy units , meaning there is no crossing of the potentials.It might be argued that there might be improvements of the results due to a choice of another type of switching potential, involving, for example, exponentials or the hyperbolic tanh function.The problem however, is not the smoothness of the curves and the degree of continuity with ever smaller energy gaps between states but the fact that finite time steps are used, and that the cross-over trajectory between different states of the particles from dimer to free particle and vice versais calculated according to potential before the bifurcation, so that an a posteriori calculation or algorithm must be invented to scale the velocites in such a way as to be consistent with the new potentials that are operating after the switch and transition.The current algorithm is applied for both these types of cross-over regions.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 δ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, h and 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 distributions is zero, that is, l ≈ h ≈ 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 '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 11 cannot be employed for this type of study.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 where W − 1/3 i j>i w r ij and the intermolecular pair Virial w r is given by w r r dv r /dr with v being the potential.

Algorithm and Analysis of Numerical Results
The velocity Verlet algorithm 25, page 81 used here 17 and allied types generate a trajectory at time nδt from that at n − 1 δt with step increment δt through a mapping T m , where v nδt , r nδt T m v n − 1 δt , r n − 1 δt which does not scale linearly with δt.This follows from the form used here consisting of 3 steps in computing the trajectory at time t δt from the data at time t:

2.1
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 NEMD is that for the nth step, ΔE n |H nδt − E| ≤ and also N i 1 ΔE i ≤ s for the specified s.In the simulation under NEMD, 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 δ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 from φ b to φ a at an interparticle distance icd r c such as points r f and r b of Figure 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.5 , a bifurcation occurs where the MD program computes the next step coordinates as for the unreacted system potential φ b , which needs to be corrected.Let the icd at time step i be r i with φ b potential and at step i 1 after interval δt be r f r i 1 where r f < r c < r i .Due to this crossover, a different Hamiltonian H is operative after point r c is crossed, where under NEMD, 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 φ a r f for φ a and E b φ b r f for φ b , where Δ E b − E a .Then if r f be the final coordinate due to the φ 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 φ b trajectory calculation?and ii can a pseudostochastic 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 the following.

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 δ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, F ext,i • Δr i δW i where Δr i is the distance increment during at least part of the time step from r c to r f .For the nonreacting trajectory over time λδt λ ≤ 1 virtual because it is not the correct path due to the crossover at r c , where Δ 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 nonreactive transition leading to the transition where Δ 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, that is, with momentum conservation, that is, δV vir r i δφ a r i for the variation along the r i coordinate, but δφ a r i −δK.E.along internuclear coordinate r i whereas δV vir −K.E.scaled about all three axes .Continuity of potential implies Subtracting 2.2 from 2.3 and applying b.c.'s 2.5 leads to

2.6
The above shows that a conservative virtual potential could be said to be operating in the vicinity of the transition from t c to t a .
Question i above leads to the following.

Mathematical Problems in Engineering
Theorem 2.2.Relative to the velocities at any r f due to the φ b potential, the rescaled velocities v due to the potential difference Δ leading to these final velocities due to the virtual potential can have a form given by (where i 1, 2) for a vector β.
Proof.From the v velocities at r f due to φ b , we compute the v 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 v and v in 2.7 yields a definition of β summation from 1 to 2 for what follows, where the mass of particle i is m i Defining for any vector s, s 2 s.s, β 2 α 2 Q, where then the rescaled velocities become from 2.7 With Δ E b − E a , energy conservation implies The coupling of 2.10 and 2.11 leads, after several steps of algebra to

2.12
Defining a v 1 − v 2 2 , q m 2 m 1 / 2 m 1 m 2 , q > 0, a ≥ 0 , then the above is equivalent to the quadratic equation: and in simulations, only α is unknown and can be determined from 2.13 where real solutions exist for Δ/qa ≥ −1.The above Inequality leads to a certain asymmetry concerning forward and backward reactions, even for reversible reactions where the regions of formation and breakdown of molecules are located in the same region with the reversal of the sign of approximate Δ.For this simulation, a reaction in either direction formation or breakdown of dimer proceeds if 2.12 is true for real α; if not, then the trajectory follows the one for the initial trajectory without any reaction i.e., no potential crossover .We would like to suggest that the real reasons for shifted potentials showing "instability in the numerical solution of the differential equations" 25, page 146, line 5 has nothing to do with the forces being discontinuous.It will be recalled that this potential v S has the form v S r ij v r ij − v c for 0 < r ij ≤ r c , and v S 0 otherwise with v c v r c .This is because the potentials are continuous and by Newton's Third Law there can be no net change in the momentum due to intermolecular forces implying momentum conservation.Further, the energies both kinetic and that for the continuous potential cannot change over an instantaneous change of the forces over zero distance.Thus there is also energy conservation.The reason for the instabilities is due to the fact that the change of position is calculated from the forces of the previous step before the sudden change in force, that puts the particle position away from the PES where there is no mechanical algorithm to correct for the violation in energy conservation with respect to the PES and the kinetic energy.Hence, the problem has nothing to do with the mechanical or dynamical setup of the potential and the forces, but with the MD move algorithm that cannot handle effectively discontinuities of the forces.

Interpretation of Results
Figure 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 nonsynthetic 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 about coordinates r f and r b , curves l1 and l2 have too large δt value to achieve equilibrium-meaning flat or invariant-temperature see Figure 2 or pressure see Figure 3 or unit step thermostat heat supply see Table 1 .h and l profiles where for these curves, the h , l values show net heat absorption.
The curve at t1 with δt 5.0 ep − 5 shows 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 18 .For this δ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 nonconservation leading to a nonflat temperature distribution, as observed and which was used to determine the regime of validity of the algorithm.

Conclusions
Without difficulty, one can easily construct a reversible system where r f and r b coincide, and this will be investigated next.Such systems would typically have most of the particles in the molecular or dimer state, and accumulated machine computational errors would be one factor to consider which this algorithm should effectively address.The two body potentials considered here saves time but the methodology is general and applies to all n-body interactions, because the essential kinetics and dynamics of all physical phenomena are governed by the principle of conservation of energy and momentum without exception.This element has often been bypassed or has received little emphasis in non-Hamiltonian and other synthetic methodologies used currently.

Figure 1 :
Figure 1: Potentials used for this work.

,
switching functions SW to demarcate potential boundaries 23, equations 5, 6 about bonding angles θ and bond distances r having the forms below have been used: SW θ abc 1 − cos16 θ abc 4 SW r ab 1 − tanh a r ab − R e r ab b 8 .

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