Time reversible and symplectic integrators for molecular dynamics simulations of rigid molecules

Molecular dynamics integrators are presented for translational and rotational motion of rigid molecules in microcanonical, canonical, and isothermal-isobaric ensembles. The integrators are all time reversible and are also, in some approaches, symplectic for the microcanonical ensembles. They are developed utilizing the quaternion representation on the basis of the Trotter factorization scheme using a Hamiltonian formalism. The structure is similar to that of the velocity Verlet algorithm. Comparison is made with standard integrators in terms of stability and it is found that a larger time step is stable with the new integrators. The canonical and isothermal-isobaric molecular dynamics simulations are defined by using a chain thermostat approach according to generalized Nosé-Hoover and Andersen methods.     

New open modified trigonometrically-fitted Newton-Cotes type multilayer symplectic integrators for the numerical solution of the Schrödinger equation

In this paper we present new modified open Newton Cotes integrators and we develop a new modified trigonometrically-fitted open Newton-Cotes method. We study the connection between the new proposed schemes, the differential methods and the symplectic integrators. although The research on multistep symplectic integrators is very poor, although, much research has been done on one step symplectic integrators and several of then have obtained based on symplectic geometry. In this paper a new open modified numerical algorithm of Newton-Cotes type is produced. We present the new obtained method as symplectic multilayer integrator. The new obtained symplectic schemes are applied for the solution of the resonance problem of the radial Schrödinger Equation. The results show the efficiency of the new proposed algorithm.     

Explicit symplectic integrators of molecular dynamics algorithms for rigid-body molecules in the canonical, isobaric-isothermal, and related ensembles

The authors propose explicit symplectic integrators of molecular dynamics (MD) algorithms for rigid-body molecules in the canonical and isobaric-isothermal ensembles. They also present a symplectic algorithm in the constant normal pressure and lateral surface area ensemble and that combined with the Parrinello-Rahman algorithm. Employing the symplectic integrators for MD algorithms, there is a conserved quantity which is close to Hamiltonian. Therefore, they can perform a MD simulation more stably than by conventional nonsymplectic algorithms. They applied this algorithm to a TIP3P pure water system at 300 K and compared the time evolution of the Hamiltonian with those by the nonsymplectic algorithms. They found that the Hamiltonian was conserved well by the symplectic algorithm even for a time step of 4 fs. This time step is longer than typical values of 0.5-2 fs which are used by the conventional nonsymplectic algorithms.     

Configurational constant pressure molecular dynamics

We propose two new algorithms for generating isothermal-isobaric molecular dynamics. The algorithms are based on an extended phase space dynamics where two extra degrees of freedom, representing the thermostat and the barostat, are included. These new methods adopt a totally different approach towards molecular dynamics simulation in the isothermal-isobaric ensemble. They are fully configurational in the sense that only the particle positions are required in the control of the system temperature and pressure. Following on from the works of Delhommelle and Evans [Mol. Phys., 99, 1825 (2001)] and of Braga and Travis [J. Chem. Phys., 123, 134101 (2005)] concerning configurational canonical dynamics, these new algorithms can be seen as a natural extension to the isothermal-isobaric ensemble. We have validated both of our new configurational isothermal-isobaric schemes by conducting molecular dynamics simulations of a Lennard-Jones fluid and comparing the static and dynamic properties for a single state point. We find that both schemes generate similar results compared with schemes which use kinetic temperature and pressure control. We have also monitored the response of the system to a series of isothermal compressions and isobaric quenches. We find that the configurational schemes performed at least as well as the kinetic based scheme in bringing the system temperature and pressure into line with the set point values of these variables. These new methods will potentially play a significant role in simulations where the calculation of the kinetic temperature and pressure can be problematic. A well known example resides in the field of nonequilibrium simulations where the kinetic temperature and pressure require a knowledge of the streaming velocity of the fluid in order to calculate the true peculiar velocities (or momenta) that enter into their definitions. These are completely avoided by using our configurational thermostats and barostats, since these are independent of momenta. By extending the analysis of Kusnezov et al. [Ann. Phys., 204, 155 (1990)] in order to derive a set of generalized Nose-Hoover equations of motion which can generate isothermal-isobaric dynamics in a number of different ways, we are able to show that both of our new configurational barostats and Hoover's kinetic isothermal-isobaric scheme are special cases of this more general set of equations. This generalization can be very powerful in generating constant pressure dynamics for a variety of systems.     

Field response and switching times in biaxial nematics

We study by means of virtual molecular dynamics computer experiments the response of a bulk biaxial nematic to an applied external field and, in particular, the relative speed of reorientation of the principal director axis and of the secondary one, typical of these new materials, upon a pi2 field switch. We perform the simulations setting up and integrating the equations of motion for biaxial Gay-Berne particles using quaternions and a suitable time reversible symplectic integrator. We find that switching of the secondary axis is up to an order of magnitude faster than that of the principal axis, and that under fields above a certain strength a reorganization of local domains, temporarily disrupting the nematic and biaxial ordering, rather than a collective concerted reorientation occurs.     

High‐order symplectic integration in quasi‐classical trajectory simulation: Case study for O(1D) + H2

Classical trajectory calculations for the O(¹D) + H₂ reaction system are employed to assess the effectiveness of the symplectic integrators. The sixth‐order symplectic integrator has been found to be the most suitable method for the quasi‐classical trajectory calculation of a long‐lived complex‐forming reaction system. In comparison with the traditional fourth‐order Runge–Kutta initialized fourth‐order Admas–Moulton–Hamming predictor‐corrector integrator (RK4‐AMH4), the sixth‐order symplectic integrator is six times less computationally expensive and exhibits better energy conservation. © 2006 Wiley Periodicals, Inc. Int J Quantum Chem, 2006     

Ab initio Ehrenfest dynamics

We present an ab initio direct Ehrenfest dynamics scheme using a three time-step integrator. The three different time steps are implemented with nuclear velocity Verlet, nuclear-position-coupled midpoint Fock integrator, and time-dependent Hartree-Fock with a modified midpoint and unitary transformation algorithm. The computational cost of the ab initio direct Ehrenfest dynamics presented here is found to be only a factor of 2-4 larger than that of Born-Oppenheimer (BO) dynamics. As an example, we compute the vibration of the NaCl molecule and the intramolecular torsional motion of H2C=NH2+ by Ehrenfest dynamics compared with BO dynamics. For the vibration of NaCl with an initial kinetic energy of 1.16 eV, Ehrenfest dynamics converges to BO dynamics with the same vibrational frequency. The intramolecular rotation of H2C=NH2+ produces significant electronic excitation in the Ehrenfest trajectory. The amount of nonadiabaticity, suggested by the amplitude of the coherent progression of the excited and ground electronic states, is observed to be directly related to the strength of the electron-nuclear coupling. Such nonadiabaticity is seen to have a significant effect on the dynamics compared with the adiabatic approximation.     

Geometric integration in Born-Oppenheimer molecular dynamics

Geometric integration schemes for extended Lagrangian self-consistent Born-Oppenheimer molecular dynamics, including a weak dissipation to remove numerical noise, are developed and analyzed. The extended Lagrangian framework enables the geometric integration of both the nuclear and electronic degrees of freedom. This provides highly efficient simulations that are stable and energy conserving even under incomplete and approximate self-consistent field (SCF) convergence. We investigate three different geometric integration schemes: (1) regular time reversible Verlet, (2) second order optimal symplectic, and (3) third order optimal symplectic. We look at energy conservation, accuracy, and stability as a function of dissipation, integration time step, and SCF convergence. We find that the inclusion of dissipation in the symplectic integration methods gives an efficient damping of numerical noise or perturbations that otherwise may accumulate from finite arithmetics in a perfect reversible dynamics.     

Langevin dynamics for rigid bodies of arbitrary shape

We present an algorithm for carrying out Langevin dynamics simulations on complex rigid bodies by incorporating the hydrodynamic resistance tensors for arbitrary shapes into an advanced rotational integration scheme. The integrator gives quantitative agreement with both analytic and approximate hydrodynamic theories for a number of model rigid bodies and works well at reproducing the solute dynamical properties (diffusion constants and orientational relaxation times) obtained from explicitly solvated simulations.     

Rigid-body dynamics in the isothermal-isobaric ensemble: a test on the accuracy and computational efficiency

We have developed a time-reversible rigid-body (rRB) molecular dynamics algorithm in the isothermal-isobaric (NPT) ensemble. The algorithm is an extension of rigid-body dynamics [Matubayasi and Nakahara, J Chem Phys 1999, 110, 3291] to the NPT ensemble on the basis of non-Hamiltonian statistical mechanics [Martyna, G. J. et al., J Chem Phys 1994, 101, 4177]. A series of MD simulations of water as well as fully hydrated lipid bilayer systems have been undertaken to investigate the accuracy and efficiency of the algorithm. The rRB algorithm was shown to be superior to the state-of-the-art constraint-dynamics algorithm SHAKE/RATTLE/ROLL, with respect to computational efficiency. However, it was revealed that both algorithms produced accurate trajectories of molecules in the NPT as well as NVT ensembles, as long as a reasonably short time step was used. A couple of multiple time-step (MTS) integration schemes were also examined. The advantage of the rRB algorithm for computational efficiency increased when the MD simulation was carried out using MTS on parallel processing computer systems; total computer time for MTS-MD of a lipid bilayer using 64 processors was reduced by about 40% using rRB instead of SHAKE/RATTLE/ROLL.     

Inherent speedup limitations in multiple time step/particle mesh Ewald algorithms

Multiple time step (MTS) algorithms present an effective integration approach to reduce the computational cost of dynamics simulations. By using force splitting to allow larger time steps for the more slowly varying force components, computational savings can be realized. The Particle-Mesh-Ewald (PME) method has been independently devised to provide an effective and efficient treatment of the long-range electrostatics interactions. Here we examine the performance of a combined MTS/PME algorithm previously developed for AMBER on a large polymerase beta/DNA complex containing 40,673 atoms. Our goal is to carefully combine the robust features of the Langevin/MTS (LN) methodology implemented in CHARMM-which uses position rather than velocity Verlet with stochasticity to make possible outer time steps of 150 fs-with the PME formulation. The developed MTS/PME integrator removes fast terms from the reciprocal-space Ewald component by using switch functions. We analyze the advantages and limitations of the resulting scheme by comparing performance to the single time step leapfrog Verlet integrator currently used in AMBER by evaluating different time-step protocols using three assessors for accuracy, speedup, and stability, all applied to long (i.e., nanosecond) simulations to ensure proper energy conservation. We also examine the performance of the algorithm on a parallel, distributed shared-memory computer (SGI Origin 2000 with 8 300-MHz R12000 processors). Good energy conservation and stability behavior can be demonstrated, for Newtonian protocols with outer time steps of up to 8 fs and Langevin protocols with outer time steps of up to 16 fs. Still, we emphasize the inherent limitations imposed by the incorporation of MTS methods into the PME formulation that may not be widely appreciated. Namely, the limiting factor on the largest outer time-step size, and hence speedup, is an intramolecular cancellation error inherent to PME. This error stems from the excluded-nonbonded correction term contained in the reciprocal-space component. This cancellation error varies in time and introduces artificial frequencies to the governing dynamics motion. Unfortunately, we find that this numerical PME error cannot be easily eliminated by refining the PME parameters (grid resolution and/or order of interpolating polynomial). We suggest that methods other than PME for fast electrostatics may allow users to reap the full advantages from MTS algorithms.     

Molecular dynamics in the isothermal-isobaric ensemble: the requirement of a "shell" molecule. III. Discontinuous potentials

Based on the approach of Gruhn and Monson [Phys. Rev. E 63, 061106 (2001)], we present a new method for deriving the collisions dynamics for particles that interact via discontinuous potentials. By invoking the conservation of the extended Hamiltonian, we generate molecular dynamics (MD) algorithms for simulating the hard-sphere and square-well fluids within the isothermal-isobaric (NpT) ensemble. Consistent with the recent rigorous reformulation of the NpT ensemble partition function, the equations of motion impose a constant external pressure via the introduction of a shell particle of known mass [M. J. Uline and D. S. Corti, J. Chem. Phys. 123, 164101 (2005); 123, 164102 (2005)], which serves to define uniquely the volume of the system. The particles are also connected to a temperature reservoir through the use of a chain of Nose-Hoover thermostats, the properties of which are not affected by a hard-sphere or square-well collision. By using the Liouville operator formalism and the Trotter expansion theorem to integrate the equations of motion, the update of the thermostat variables can be decoupled from the update of the positions of the particles and the momentum changes upon a collision. Hence, once the appropriate collision dynamics for the isobaric-isenthalpic (NpH) equations of motion is known, the adaptation of the algorithm to the NpT ensemble is straightforward. Results of MD simulations for the pure component square-well fluid are presented and serve to validate our algorithm. Finally, since the mass of the shell particle is known, the system itself, and not a piston of arbitrary mass, controls the time scales for internal pressure and volume fluctuations. We therefore consider the influence of the shell particle algorithm on the dynamics of the square-well fluid.     

Partial rigid-body dynamics in NPT, NPAT and NPgammaT ensembles for proteins and membranes

A partial rigid-body method of molecular dynamics simulations for proteins and membranes is presented. In this method, the symplectic integrator for rigid bodies is combined with the equations of motion for the NPT ensemble. The standard NPT ensemble is extended to the membrane-specific ensembles, the NPAT (constant normal pressure and lateral surface area of membranes and constant temperature) and NPgammaT (constant normal pressure and lateral surface tension of membranes and constant temperature) ensembles. By more than 30-ns simulations of aqueous proteins and hydrated lipid bilayers, the results of the partial rigid-body method demonstrated excellent conservation of total energy and consistent behavior with the traditional constraint method in terms of structural distribution and fluctuation of proteins and lipids. The efficient implementation of the partial rigid-body method in parallel computation is presented, which is shown to work well in large-scale molecular dynamics simulations.     

A_2B模型分子经典轨迹的辛算法计算

采用辛算法计算了Ａ２Ｂ模型分子的经典轨迹并与传统Ｒｕｎｇｅ－Ｋｕｔｔａ　（Ｒ－Ｋ）算法进行了比较．结果表明，在微观反应动力学研究所应考虑的时间范围内，辛算法的结果与理论分析一致，Ｒ－Ｋ法的结果则面目全非．因此，用辛算法替代传统数值方法有可能克服目前经典轨迹计算存在的困难，从根本上改进微观反应动力学研究的经典轨迹方法．     

Molecular dynamics integration meets standard theory of molecular vibrations

An iterative SISM (split integration symplectic method) for molecular dynamics (MD) integration is described. This work explores an alternative for the internal coordinate system prediction in the SISM introduced by JaneZic et al. (J. Chem. Phys. 2005, 122, 174101). The SISM, which employs a standard theory of molecular vibrations, analytically resolves the internal high-frequency molecular vibrations. This is accomplished by introducing a translating and rotating internal coordinate system of a molecule and calculating normal modes of an isolated molecule only. The Eckart frame, which is usually used in the standard theory of molecular vibrations as an internal coordinate system of a molecule, is adopted to be used within the framework of the second order generalized leapfrog scheme. In the presented MD integrator the internal coordinate frame at the end of the integration step is predicted halfway through the integration step using a predictor-corrector type iterative approach thus ensuring the method's time reversibility. The iterative SISM, which is applicable to any system of molecules with one equilibrium configuration, was applied here to perform all-atom MD simulations of liquid CO2 and SO2. The simulation results indicate that for the same level of accuracy, this algorithm allows significantly longer integration time steps than the standard second-order leapfrog Verlet (LFV) method.     

Direct evaluation of multicomponent phase equilibria using flat-histogram methods

We present a method for directly locating density-driven phase transitions in multicomponent systems. Phase coexistence conditions are determined through manipulation of a total density probability distribution evaluated over a density range that includes both coexisting phases. Saturation quantities are determined through appropriate averaging of density-dependent mean values of a given property of interest. We discuss how to implement the method in both the grand-canonical and isothermal-isobaric semigrand ensembles. Calculations can be conducted using any of the recently introduced flat-histogram techniques. Here, we combine the general algorithm with a transition-matrix approach to produce an efficient self-adaptive technique for determining multicomponent phase equilibrium properties. To assess the performance of the new method, we generate phase diagrams for a number of binary and ternary Lennard-Jones mixtures.     

A cutoff phenomenon in accelerated stochastic simulations of chemical kinetics via flow averaging (FLAVOR-SSA)

We present a simple algorithm for the simulation of stiff, discrete-space, continuous-time Markov processes. The algorithm is based on the concept of flow averaging for the integration of stiff ordinary and stochastic differential equations and ultimately leads to a straightforward variation of the the well-known stochastic simulation algorithm (SSA). The speedup that can be achieved by the present algorithm [flow averaging integrator SSA (FLAVOR-SSA)] over the classical SSA comes naturally at the expense of its accuracy. The error of the proposed method exhibits a cutoff phenomenon as a function of its speed-up, allowing for optimal tuning. Two numerical examples from chemical kinetics are provided to illustrate the efficiency of the method.     

Molecular dynamics comparative study of methane–nitrogen and methane–nitrogen–ethane systems

This work concerns the site–site interaction study of 256 particles using the Buckingham potential model. We have calculated the new parameters of the Buckingham potential using an iterative algorithm with a mean square method. This adapted model allows determining the characteristics for each state point. We have applied this model to study the liquefied natural gas LNG properties for methane-nitrogen and methane–nitrogen–ethane mixtures by molecular dynamics. We have calculated the thermodynamic, dynamic and structural properties for both the microcanonical NVT and the isothermal-isobaric NPT ensembles of binary and ternary systems from the SP₁ to SP₉ points. Then, we have compared the results between binary and ternary systems. We have obtained a good prediction on transport properties. From the calculated values of self-diffusion coefficient and viscosity, we have confirmed the liquid state of the liquefied natural gas LNG system.     

Enthalpy-entropy contributions to the potential of mean force of nanoscopic hydrophobic solutes

Entropic and enthalpic contributions to the hydrophobic interaction between nanoscopic hydrophobic solutes, modeled as graphene plates in water, have been calculated using molecular dynamics simulations in the isothermal-isobaric (NPT) ensemble with free energy perturbation methodology. We find the stabilizing contribution to the free energy of association (contact pair formation) to be the favorable entropic part, the enthalpic contribution being highly unfavorable. The desolvation barrier is dominated by the unfavorable enthalpic contribution, despite a fairly large favorable entropic compensation. The enthalpic contributions, incorporating the Lennard-Jones solute-solvent terms, largely determine the stability of the solvent separated configuration. We decompose the enthalpy into a direct solute-solute term, the solute-solvent interactions, and the remainder that contains pressure-volume work as well as contributions due to solvent reorganization. The enthalpic contribution due to changes in water-water interactions arising from solvent reorganization around the solute molecules is shown to have major contribution in the solvent induced enthalpy change.     

Molecular simulation study of cavity-generated instabilities in the superheated Lennard-Jones liquid

Previous equilibrium-based density-functional theory (DFT) analyses of cavity formation in the pure component superheated Lennard-Jones (LJ) liquid [S. Punnathanam and D. S. Corti, J. Chem. Phys. 119, 10224 (2003); M. J. Uline and D. S. Corti, Phys. Rev. Lett. 99, 076102 (2007)] revealed that a thermodynamic limit of stability appears in which no liquidlike density profile can develop for cavity radii greater than some critical size (being a function of temperature and bulk density). The existence of these stability limits was also verified using isothermal-isobaric Monte Carlo (MC) simulations. To test the possible relevance of these limits of stability to a dynamically evolving system, one that may be important for homogeneous bubble nucleation, we perform isothermal-isobaric molecular dynamics (MD) simulations in which cavities of different sizes are placed within the superheated LJ liquid. When the impermeable boundary utilized to generate a cavity is removed, the MD simulations show that the cavity collapses and the overall density of the system remains liquidlike, i.e., the system is stable, when the initial cavity radius is below some certain value. On the other hand, when the initial radius is large enough, the cavity expands and the overall density of the system rapidly decreases toward vaporlike densities, i.e., the system is unstable. Unlike the DFT predictions, however, the transition between stability and instability is not infinitely sharp. The fraction of initial configurations that generate an instability (or a phase separation) increases from zero to unity as the initial cavity radius increases over a relatively narrow range of values, which spans the predicted stability limit obtained from equilibrium MC simulations. The simulation results presented here provide initial evidence that the equilibrium-based stability limits predicted in the previous DFT and MC simulation studies may play some role, yet to be fully determined, in the homogeneous nucleation and growth of embryos within metastable fluids.