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

Abstract

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.

1. 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 assemblies encountered in Biophysical simulation where 𝑛-body (𝑛>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 𝑛-body interactions for the coordinates involved in any cross-over 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 𝐻 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 𝑇 too, as determined by the equipartition result 𝑛(3/2)𝑘𝑇=∑^𝑛_𝑖=1⟨(1/2)𝑚_𝑖𝑣²_𝑖⟩. 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 𝐻 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 {𝐫̇_𝑖,𝐩̇_𝑖} as ℒ_Nose=^𝑁_𝑖=1𝑚_𝑖2𝑠²𝐫̇²_𝑖−𝒰(𝐫^𝑁𝑄)+2𝑠²̇−𝐿𝛽ln𝑠, (1.1) where 𝛽 is the temperature parameter. This so-called Lagrangian defines the conjugate momenta to 𝐫_𝑖 and 𝑠 as, respectively, 𝐩_𝑖=𝑚_𝑖𝑠²𝐫̇_𝑖 and 𝑝_𝑠=𝑄𝑠̇. Then for this system, there results ultimately a pseudo-Hamiltonian: ℋ_Nose=^𝑁_𝑖=1𝑚_𝑖2𝑠²𝐩²_𝑖−𝒰(𝐫^𝑁𝑄)+2𝜉²+𝐿𝛽ln𝑠, (1.2) whose trajectory is determined by the coupled equations that must be solved: 𝐫̇_𝑖=𝐩_𝑖𝑚_𝑖,𝐩̇_𝑖=𝜕𝒰(𝐫_𝑖)𝜕𝐫̇_𝑖−𝜉𝐩_𝑖,𝜉̇=(∑_𝑖𝑝²_𝑖/𝑚−𝐿/𝛽)𝑄,𝑠̇𝑠=𝑑ln𝑠𝑑𝑡=𝜉, (1.3) 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 𝜆=[1+𝛿𝑡𝜏_𝑇(𝑇₀𝑇−1)]^1/2. (1.4) 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 𝐻 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 𝛿(𝐻−𝐸+𝛼). 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 𝑛-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 𝑛-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 2𝐴^𝑘1_⇄𝑘−1𝐴₂, (1.5) where 𝑘₁ and 𝑘₋₁ 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 𝑇^∗ (all units used are reduced LJ ones [17]) ranges ∼0.3–1.2, [17] whereas here, 𝑇^∗∼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 𝑟_𝑏=1.20 for the vicinity where the bond of the dimer is broken (2 free particles emerge) and 𝑟_𝑓=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 𝑢_LJ except for the dimeric pair (𝑖,𝑗) formed from particles 𝑖 and 𝑗 which interact with a harmonic-like intermolecular potential modified by a switch 𝑢(𝑟) given by 𝑢(𝑟)=𝑢_vib(𝑟)𝑠(𝑟)+𝑢_LJ[1−𝑠(𝑟)], (1.6) where 𝑢_vib(𝑟) is the vibrational potential given by (1.7) 𝑢_vib(𝑟)=𝑢₀+12𝑘(𝑟−𝑟₀)². (1.7) The switching function 𝑠(𝑟) is defined as 1𝑠(𝑟)=1+(𝑟/𝑟_sw)^𝑛, (1.8) where 𝑠(𝑟)→1if𝑟<𝑟_sw,𝑠(𝑟)→0for𝑟>𝑟_sw. (1.9) The switching function becomes effective when the distance between the atoms approach the value 𝑟_sw (see Figure 1). Some of the other parameters used in the equations that follow include 𝑢₀=−10,𝑟₀=1.0,𝑘∼2446 (exact value is determined by the other input parameters), 𝑛=100,𝑟_𝑓=0.85,𝑟_𝑏=1.20,and𝑟_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]: 𝑆(𝑅)=1−tanh[𝑎(𝑅−𝑅_𝑒)^𝑏], (1.10) where 𝑎 and 𝑏 are defined constants and the 𝑅^′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 𝑔(𝑡) function of the following form has been used [20]: 𝑔(𝑡)=exp(−𝜁𝑡)cos²(3𝜋𝜁𝑡)(1−𝜆𝜁𝑡), (1.11) where 𝑡 is a time-dependent variable. On the other hand, for cluster dynamics, a switch of the form [21, equation 7], Φ(𝑟)=tanh(𝑟−𝐸𝐹), (1.12) is used, where the parameters 𝐸 and 𝐹 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, switching functions SW to demarcate potential boundaries [23, equations 5, 6] about bonding angles 𝜃 and bond distances 𝑟 having the forms below have been used: SW(𝜃_𝑎𝑏𝑐)=(1−cos¹⁶𝜃_𝑎𝑏𝑐)⁴,SW(𝑟_𝑎𝑏)=1−tanh[𝑎(𝑟_𝑎𝑏−𝑅_𝑒)(𝑟_𝑎𝑏+𝑏)⁸]. (1.13)

Figure 1 

Potentials used for this work.

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 𝜉(𝑡) 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 𝑖 and 𝑗 above also interact with all other particles not bonded to it via 𝑢_LJ. Details of these potentials and their interactions are given elsewhere [18]; here we note the high activation energy at 𝑟_𝑓 of approximately 17.5. At 𝑟_𝑓, 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 𝑢_LJ potentials are "very close." On the other hand, at 𝑟_𝑏, 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 versa—is 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 (𝑥 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 𝛿𝑡 at both ends of the cell (orthogonal to the 𝑥 axis) in the vicinity of 𝑥=0 and 𝑥=1 maintained at temperatures 𝑇_ℎ and 𝑇_𝑙 could be monitored, where this energy per unit step time is, respectively, 𝜖_ℎ and 𝜖_𝑙. At equilibrium (when 𝑇_ℎ=𝑇_𝑙) the net energy supplied within statistical error (meaning 1–3 units of the standard error of the 𝜖 distributions) is zero, that is, 𝜖_𝑙≈𝜖_ℎ≈0. The cell is divided up uniformly into 64 rectangular regions along the 𝑥 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 𝑇 and pressure 𝑝 are computed by the equipartition and Virial expression given, respectively, by ⟨∑^𝑁_𝑖=1𝐩_𝐢⋅𝐩_𝐢𝑚_𝑖⟩=3𝑁𝑘_𝐵𝑇,𝑃=𝜌𝑘_𝐵𝑊𝑇+𝑉, (1.14) where 𝑊=−(1/3)∑_𝑖∑_𝑗>𝑖𝑤(𝑟_𝑖𝑗) and the intermolecular pair Virial 𝑤(𝑟) is given by 𝑤(𝑟)=𝑟(𝑑𝑣(𝑟)/𝑑𝑟) with 𝑣 being the potential.

2. Algorithm and Analysis of Numerical Results

The velocity Verlet algorithm [25, page 81] used here [17] and allied types generate a trajectory at time 𝑛𝛿𝑡 from that at (𝑛−1)𝛿𝑡 with step increment 𝛿𝑡 through a mapping 𝑇_𝑚, where (𝑣(𝑛𝛿𝑡),𝑟(𝑛𝛿𝑡))=𝑇_𝑚(𝑣((𝑛−1)𝛿𝑡),𝑟((𝑛−1)𝛿𝑡)) which does not scale linearly with 𝛿𝑡. This follows from the form used here consisting of 3 steps in computing the trajectory at time 𝑡+𝛿𝑡 from the data at time 𝑡: 1𝐯(𝑡+21𝛿𝑡)=𝐯(𝑡)+21𝛿𝑡𝐚(𝑡),𝐫(𝑡+𝛿𝑡)=𝐫(𝑡)+𝛿𝑡𝐯(𝑡+21𝛿𝑡),𝐯(𝑡+𝛿𝑡)=𝐯(𝑡+21𝛿𝑡)+2𝛿𝑡𝐚(𝑡+𝛿𝑡). (2.1)

For a Hamiltonian 𝐻 whose potential 𝑉 is dependent only on position 𝑟 having momentum components 𝑝_𝑖, the system without external perturbation has constant energy 𝐸, and the normal assumption in MD (NEMD) is that for the 𝑛th step, Δ𝐸_𝑛=|𝐻(𝑛𝛿𝑡)−𝐸|≤𝜖 and also ∑^𝑁_𝑖=1Δ𝐸_𝑖≤𝜖^𝑠 for the specified 𝜖^′𝑠. 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 𝛿𝑡_𝑇 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 𝜙_𝑏 to 𝜙_𝑎 at an interparticle distance (icd) 𝑟_𝑐 (such as points 𝑟_𝑓 and 𝑟_𝑏 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 𝜙_𝑏), which needs to be corrected. Let the icd at time step 𝑖 be 𝑟_𝑖 (with 𝜙_𝑏 potential) and at step 𝑖+1 after interval 𝛿𝑡 be 𝑟_𝑓=𝑟_𝑖+1 where 𝑟_𝑓<𝑟_𝑐<𝑟_𝑖. Due to this crossover, a different Hamiltonian 𝐻^′ is operative after point 𝑟_𝑐 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 𝑟_𝑓 be 𝐸_𝑎=𝐸_𝑓=𝜙_𝑎(𝑟_𝑓) for 𝜙_𝑎 and 𝐸_𝑏=𝜙_𝑏(𝑟_𝑓) for 𝜙_𝑏, where Δ=𝐸_𝑏−𝐸_𝑎. Then if 𝑟_𝑓 be the final coordinate due to the 𝜙_𝑏 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 𝜙_𝑏 trajectory calculation? and (ii) can a pseudostochastic potential be imposed from coordinates 𝑟_𝑐 (at virtual time 𝑡_𝑐) to 𝑟_𝑓 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 𝑟_𝑐.

Proof. The external work done 𝛿𝑊 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 𝑖, 𝐅_{𝐞𝐱𝐭,𝐢}⋅Δ𝐫_𝐢=𝛿𝑊_𝑖 where Δ𝐫_𝐢 is the distance increment during at least part of the time step from 𝑟_𝑐 to 𝑟_𝑓. For the nonreacting trajectory over time 𝜆𝛿𝑡 (𝜆≤1) (virtual because it is not the correct path due to the crossover at 𝑟_𝑐), 𝛿𝑊₂+𝛿𝑊₁−(𝜙_𝑏(𝑟_𝑓)−𝜙_𝑏(𝑟_𝑐))=Δ(K.E.), (2.2) where Δ∑(K.E.) is the change of kinetic energy for the (1,2) pair from the First Law between the end points 𝑟_𝑓,𝑟_𝑐. Now over time interval 𝑡_𝑐 to 𝑡_𝑓, for the reactive trajectory, we introduce a "virtual potential" 𝑉^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 𝛿𝑊₂+𝛿𝑊₁−(𝑉^vir(𝑟_𝑓)−𝑉^vir(𝑟_𝑐))=Δ^′(K.E.), (2.3) where Δ∑^′(K.E.) is the change of kinetic energy for the pair with 𝑉^vir turned on and along this trajectory, the change of potential for 𝑉^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, 𝛿𝑉^vir(𝑟)−𝛿𝜙_𝑏(𝑟)=𝛿(Δ(K.E._{𝑥,𝑦,𝑧})−Δ^′(K.E.)_{𝑥,𝑦,𝑧}), (2.4) with momentum conservation, that is, 𝛿𝑉^vir(𝑟_𝑖)=𝛿𝜙_𝑎(𝑟_𝑖) for the variation along the 𝑟_𝑖 coordinate, but 𝛿𝜙_𝑎(𝑟_𝑖)=−𝛿K.E. along internuclear coordinate 𝑟_𝑖 whereas 𝛿𝑉^vir=−K.E. (scaled about all three axes). Continuity of potential implies 𝜙_𝑎(𝑟_𝑓)=𝑉^vir(𝑟_𝑓);𝜙_𝑎(𝑟_𝑐)=𝑉^vir(𝑟_𝑐);𝜙_𝑏(𝑟_𝑐)=𝑉^vir(𝑟_𝑐). (2.5) Subtracting (2.2) from (2.3) and applying b.c.'s (2.5) leads to Δ=𝜙_𝑏(𝑟_𝑓)−𝑉^vir(𝑟_𝑓)=𝜙_𝑏(𝑟_𝑓)−𝜙_𝑎(𝑟_𝑓)=𝐸_𝑏−𝐸_𝑎=Δ^′(K.E.)−Δ(K.E.). (2.6) The above shows that a conservative virtual potential could be said to be operating in the vicinity of the transition (from 𝑡_𝑐to𝑡_𝑎).

Question (i) above leads to the following.

Theorem 2.2. Relative to the velocities at any 𝑟_𝑓 due to the 𝜙_𝑏 potential, the rescaled velocities 𝐯^′ due to the potential difference Δ leading to these final velocities due to the virtual potential can have a form given by 𝐯^′_𝐢=(1+𝛼)𝐯_𝐢+𝜷, (2.7) (where 𝑖=1,2) for a vector 𝜷.

Proof. From the 𝐯 velocities at 𝑟_𝑓 due to 𝜙_𝑏, we compute the 𝐯^′ velocities at 𝑟_𝑓 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 𝐯^′ and 𝐯 in (2.7) yields a definition of 𝜷 (summation from 1 to 2 for what follows, where the mass of particle 𝑖 is 𝑚_𝑖) 𝜷=−𝛼∑𝑚_𝑖v_𝐢∑𝑚_𝑖. (2.8) Defining for any vector 𝐬, 𝐬²=𝐬.𝐬,𝜷²=𝛼²𝑄, where 𝑄=(∑𝑚_𝑖𝐯_𝐢)²(∑𝑚_𝑖)², (2.9) then the rescaled velocities become from (2.7) 𝐯^′_𝐢²=(1+𝛼)²𝐯_𝐢²+2(1+𝛼)𝐯_𝐢⋅𝜷+𝜷². (2.10) With Δ=𝐸_𝑏−𝐸_𝑎, energy conservation implies 12𝑚_𝑖𝐯^′_𝐢²1−2𝑚_𝑖𝐯_𝐢²=Δ. (2.11) The coupling of (2.10) and (2.11) leads, after several steps of algebra to 𝛼Δ=²𝑚₁𝑚₂2(𝑚₁+𝑚₂)[𝐯₁²+𝐯₂²−2𝐯₁⋅𝐯₂]+2𝛼𝑚₁𝑚₂2(𝑚₁+𝑚₂)[𝐯₁²+𝐯₂²−2𝐯₁⋅𝐯₂]. (2.12) Defining 𝑎=(𝐯₁−𝐯₂)², 𝑞=𝑚₂𝑚₁/[2(𝑚₁+𝑚₂)], (𝑞>0,𝑎≥0), then the above is equivalent to the quadratic equation: 𝛼²𝑞𝑎+2𝑞𝑎𝛼−Δ=0, (2.13) and in simulations, only 𝛼 is unknown and can be determined from (2.13) where real solutions exist for Δ/𝑞𝑎≥−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 𝑣^𝑆 has the form 𝑣^𝑆(𝑟_𝑖𝑗)=𝑣(𝑟_𝑖𝑗)−𝑣_𝑐 for 0<𝑟_𝑖𝑗≤𝑟_𝑐, and 𝑣^𝑆=0 otherwise with 𝑣_𝑐=𝑣(𝑟_𝑐). 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 𝑟_𝑓 and 𝑟_𝑏, curves 𝑙1 and 𝑙2 have too large 𝛿𝑡 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). (𝜖_ℎ and 𝜖_𝑙) profiles where for these curves, the (𝜖_ℎ,𝜖_𝑙) values show net heat absorption.

Table 1 

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

Figure 2 

Temperature profile across the cell for different set conditions 𝑎−𝑒 for temperature 𝑇∗ and step time 𝛿𝑡 pairs (𝑇∗,𝛿𝑡) where 𝑎=(8.0,2.0𝑒𝑝−3),𝑏=(8.0,5.0𝑒𝑝−4),𝑐=(8.0,5.0𝑒𝑝−5),𝑑=(12.0,5.0𝑒𝑝−5),𝑒=(16.0,5.0𝑒𝑝−5). The curves {𝑙1,𝑙3,𝑡1,𝑡2,𝑡3} results with the application of the algorithm at 𝑟𝑏 and 𝑟𝑓 with associated conditions 𝑙1⇔𝑎,𝑙3⇔𝑏,𝑡1⇔𝑐,𝑡2⇔𝑑,𝑡3⇔𝑒 while the curves {𝑙2,𝑙4,𝑙5,𝑙6,𝑙7} are for the cases without implementing the algorithm with the associated conditions 𝑙2⇔𝑎, 𝑙4⇔𝑏, 𝑙5⇔𝑐, 𝑙6⇔𝑑, 𝑙7⇔𝑒, where 𝑒𝑝𝑥≡10𝑥.

Figure 3 

Pressure profile across the cell for different runs. The conditions of the runs and the labeling of the curves are exactly as in Figure 2.

The curve at 𝑡1 with 𝛿𝑡=5.0𝑒𝑝−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 (𝑇=12 and 𝑇=16) which was used to determine thermodynamical properties [18]. For this 𝛿𝑡 value and all others, no reasonable stationary equilibrium conditions could be obtained without the application of the algorithm (curves 𝑙2, 𝑙4, 𝑙5, 𝑙6 and 𝑙7). 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.

3. Conclusions

Without difficulty, one can easily construct a reversible system where 𝑟_𝑓 and 𝑟_𝑏 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 𝑛-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.

Acknowledgment

The author is grateful to University of Malaya for a Conference grant to present this algorithm as an Invited Speaker at the Fifth ICDSA (2007, Atlanta) which this communication briefly reviews.