Molecular dynamics integration and molecular vibrational theory. II. Simulation of nonlinear molecules

A series of molecular dynamics (MD) simulations of nonlinear molecules has been performed to test the efficiency of newly introduced semianalytical second-order symplectic time-reversible MD integrators that combine MD and the standard theory of molecular vibrations. The simulation results indicate that for the same level of accuracy, the new algorithms allow significantly longer integration time steps than the standard second-order symplectic leap-frog Verlet method. Since the computation cost per integration step using new MD integrators with longer time steps is approximately the same as for the standard method, a significant speed-up in MD simulation is achieved.     

Molecular dynamics integration and molecular vibrational theory. III. The infrared spectrum of water

The new symplectic molecular dynamics (MD) integrators presented in the first paper of this series were applied to perform MD simulations of water. The physical properties of a system of flexible TIP3P water molecules computed by the new integrators, such as diffusion coefficients, orientation correlation times, and infrared (IR) spectra, are in good agreement with results obtained by the standard method. The comparison between the new integrators' and the standard method's integration time step sizes indicates that the resulting algorithm allows a 3.0 fs long integration time step as opposed to the standard leap-frog Verlet method, a sixfold simulation speed-up. The accuracy of the method was confirmed, in particular, by computing the IR spectrum of water in which no blueshifting of the stretching normal mode frequencies is observed as occurs with the standard method.     

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.     

Symplectic integrators and the conservation of angular momentum

In this article we observe that generally symplectic integrators conserve angular momentum exactly, whereas nonsymplectic integrators do not. We show that this observation extends to multiple timesteps and to constrained dynamics. Both of these devices are important for efficient molecular dynamics simulations. © 1995 by John Wiley & Sons, Inc.     

Molecular dynamics with the united-residue model of polypeptide chains. I. Lagrange equations of motion and tests of numerical stability in the microcanonical mode

The Lagrange formalism was implemented to derive the equations of motion for the physics-based united-residue (UNRES) force field developed in our laboratory. The C(alpha)...C(alpha) and C(alpha)...SC (SC denoting a side-chain center) virtual-bond vectors were chosen as variables. The velocity Verlet algorithm was adopted to integrate the equations of motion. Tests on the unblocked Ala(10) polypeptide showed that the algorithm is stable in short periods of time up to the time step of 1.467 fs; however, even with the shorter time step of 0.489 fs, some drift of the total energy occurs because of momentary jumps of the acceleration. These jumps are caused by numerical instability of the forces arising from the U(rot) component of UNRES that describes the energetics of side-chain-rotameric states. Test runs on the Gly(10) sequence (in which U(rot) is not present) and on the Ala(10) sequence with U(rot) replaced by a simple numerically stable harmonic potential confirmed this observation; oscillations of the total energy were observed only up to the time step of 7.335 fs, and some drift in the total energy or instability of the trajectories started to appear in long-time (2 ns and longer) trajectories only for the time step of 9.78 fs. These results demonstrate that the present U(rot) components (which are statistical potentials derived from the Protein Data Bank) must be replaced with more numerically stable functions; this work is under way in our laboratory. For the purpose of our present work, a nonsymplectic variable-time-step algorithm was introduced to reduce the energy drift for regular polypeptide sequences. The algorithm scales down the time step at a given point of a trajectory if the maximum change of acceleration exceeds a selected cutoff value. With this algorithm, the total energy is reasonably conserved up to a time step of 2.445 fs, as tested on the unblocked Ala(10) polypeptide. We also tried a symplectic multiple-time-step reversible RESPA algorithm and achieved satisfactory energy conservation for time steps up to 7.335 fs. However, at present, it appears that the reversible RESPA algorithm is several times more expensive than the variable-time-step algorithm because of the necessity to perform additional matrix multiplications. We also observed that, because Ala(10) folds and unfolds within picoseconds in the microcanonical mode, this suggests that the effective (event-based) time unit in UNRES dynamics is much larger than that of all-atom dynamics because of averaging over the fast-moving degrees of freedom in deriving the UNRES potential.     

Implementation of a symplectic multiple-time-step molecular dynamics algorithm, based on the united-residue mesoscopic potential energy function

A symplectic multiple-time-step (MTS) algorithm has been developed for the united-residue (UNRES) force field. In this algorithm, the slow-varying forces (which contain most of the long-range interactions and are, therefore, expensive to compute) are integrated with a larger time step, termed the basic time step, and the fast-varying forces are integrated with a shorter time step, which is an integral fraction of the basic time step. Based on the split operator formalism, the equations of motion were derived. Separation of the fast- and slow-varying forces leads to stable molecular dynamics with longer time steps. The algorithms were tested with the Ala(10) polypeptide chain and two versions of the UNRES force field: the current one in which the energy components accounting for the energetics of side-chain rotamers (U(rot)) can lead to numerically unstable forces and a modified one in which the the present U(rot) was replaced by a numerically stable expression which, at present, is parametrized only for polyalanine chains. With the modified UNRES potential, stable trajectories were obtained even when extending the basic time step to 15 fs and, with the original UNRES potentials, the basic time step is 1 fs. An adaptive multiple-time-step (A-MTS) algorithm is proposed to handle instabilities in the forces; in this method, the number of substeps in the basic time step varies depending on the change of the magnitude of the acceleration. With this algorithm, the basic time step is 1 fs but the number of substeps and, consequently, the computational cost are reduced with respect to the MTS algorithm. The use of the UNRES mesoscopic energy function and the algorithms derived in this work enables one to increase the simulation time period by several orders of magnitude compared to conventional atomic-resolution molecular dynamics approaches and, consequently, such an approach appears applicable to simulating protein-folding pathways, protein functional dynamics in a real molecular environment, and dynamical molecular recognition processes.     

Langevin stabilization of molecular-dynamics simulations of polymers by means of quasisymplectic algorithms

Algorithms for the numerical integration of Langevin equations are compared in detail from the point of view of their accuracy, numerical efficiency, and stability to assess them as potential candidates for molecular-dynamics simulations of polymeric systems. Some algorithms are symplectic in the deterministic frictionless limit and prove to stabilize long time-step integrators. They are tested against other popular algorithms. The optimal algorithm depends on the main goal: accuracy or efficiency. The former depends on the observable of interest. A recently developed quasisymplectic algorithm with great accuracy in the position evaluation exhibits better overall accuracy and stability than the other ones. On the other hand, the well-known BrunGer-Brooks-Karplus [Chem. Phys. Lett. 105, 495 (1982)] algorithm is found to be faster with limited accuracy loss but less stable. It is also found that using higher-order algorithms does not necessarily improve the accuracy. Moreover, they usually require more force evaluations per single step, thus leading to poorer performances.     

Symplectic splitting operator methods for the time-dependent Schrodinger equation

We present a family of symplectic splitting methods especially tailored to solve numerically the time-dependent Schrodinger equation. When discretized in time, this equation can be recast in the form of a classical Hamiltonian system with a Hamiltonian function corresponding to a generalized high-dimensional separable harmonic oscillator. The structure of the system allows us to build highly efficient symplectic integrators at any order. The new methods are accurate, easy to implement, and very stable in comparison with other standard symplectic integrators.     

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.     

Molecular dynamics integration and molecular vibrational theory. I. New symplectic integrators

New symplectic integrators have been developed by combining molecular dynamics integration with the standard theory of molecular vibrations to solve the Hamiltonian equations of motion. The presented integrators analytically resolve the internal high-frequency molecular vibrations by introducing a translating and rotating internal coordinate system of a molecule and calculating normal modes of an isolated molecule only. The translation and rotation of a molecule are treated as vibrational motions with the vibrational frequency zero. All types of motion are thus described in terms of the normal coordinates. The method's time reversibility requirement was used to determine the equations of motion for internal coordinate system of a molecule. The calculation of long-range forces is performed numerically within the generalized second-order leap-frog scheme, in the same way as in standard second-order symplectic methods. The new methods for integrating classical equations of motion using normal mode analysis allow us to use a long integration step and are applicable to any system of molecules with one equilibrium configuration.     

Monte Carlo vs molecular dynamics for all-atom polypeptide folding simulations

An efficient Monte Carlo (MC) algorithm including concerted rotations is directly compared to molecular dynamics (MD) in all-atom statistical mechanics folding simulations of small polypeptides. The previously reported algorithm "concerted rotations with flexible bond angles" (CRA) has been shown to successfully locate the native state of small polypeptides. In this study, the folding of three small polypeptides (trpzip2/H1/Trp-cage) is investigated using MC and MD, for a combined sampling time of approximately 10(11) MC configurations and 8 micros, respectively. Both methods successfully locate the experimentally determined native states of the three systems, but they do so at different speed, with 2-2.5 times faster folding of the MC runs. The comparison reveals that thermodynamic and dynamic properties can reliably be obtained by both and that results from folding simulations do not depend on the algorithm used. Similar to previous comparisons of MC and MD, it is found that one MD integration step of 2 fs corresponds to one MC scan, revealing the good sampling of MC. The simplicity and efficiency of the MC method will enable its future use in folding studies involving larger systems and the combination with replica exchange algorithms.     

High order closed Newton-Cotes exponentially and trigonometrically fitted formulae as multilayer symplectic integrators and their application to the radial Schr?dinger equation

In this paper we study the connection between: (i) closed Newton-Cotes formulae of high order, (ii) trigonometrically-fitted and exponentially-fitted differential methods, (iii) symplectic integrators. Several one step symplectic integrators have been produced based on symplectic geometry during the last decades (see the relevant literature and the references here). However, the study of multistep symplectic integrators is very poor. In this paper we investigate the High Order Closed Newton-Cotes Formulae and we write them as symplectic multilayer structures. We develop trigonometrically-fitted and exponentially-fitted symplectic methods which are based on the closed Newton-Cotes formulae. We apply the symplectic schemes in order to solve the resonance problem of the radial Schr?dinger equation. Based on the theoretical and numerical results, conclusions on the efficiency of the new obtained methods are given.     

Implementations of Nosé-Hoover and Nosé-Poincaré thermostats in mesoscopic dynamic simulations with the united-residue model of a polypeptide chain

Molecular dynamics (MD) simulations generate a canonical ensemble only when integration of the equations of motion is coupled to a thermostat. Three extended phase space thermostats, one version of Nose-Hoover and two versions of Nose-Poincare, are compared with each other and with the Berendsen thermostat and Langevin stochastic dynamics. Implementation of extended phase space thermostats was first tested on a model Lennard-Jones fluid system; subsequently, they were implemented with our physics-based protein united-residue (UNRES) force field MD. The thermostats were also implemented and tested for the multiple-time-step reversible reference system propagator (RESPA). The velocity and temperature distributions were analyzed to confirm that the proper canonical distribution is generated by each simulation. The value of the artificial mass constant, Q, of the thermostat has a large influence on the distribution of the temperatures sampled during UNRES simulations (the velocity distributions were affected only slightly). The numerical stabilities of all three algorithms were compared with each other and with that of microcanonical MD. Both Nose-Poincare thermostats, which are symplectic, were not very stable for both the Lennard-Jones fluid and UNRES MD simulations started from nonequilibrated structures which implies major changes of the potential energy throughout a trajectory. Even though the Nose-Hoover thermostat does not have a canonical symplectic structure, it is the most stable algorithm for UNRES MD simulations. For UNRES with RESPA, the "extended system inside-reference system propagator algorithm" of the RESPA implementation of the Nose-Hoover thermostat was the only stable algorithm, and enabled us to increase the integration time step.     

Trigonometrically fitted and exponentially fitted symplectic methods for the numerical integration of the Schrödinger equation

The solution of the one-dimensional time-independent Schr?dinger equation is considered by exponentially fitted symplectic integrators. The Schr?dinger equation is first transformed into a Hamiltonian canonical equation. Numerical results are obtained for the one-dimensional harmonic oscillator and the doubly anharmonic oscillator.     

Closed Newton–Cotes trigonometrically-fitted formulae of high order for the numerical integration of the Schrödinger equation

In this paper we investigate the connection between (i) closed Newton–Cotes formulae, (ii) trigonometrically-fitted differential methods, (iii) symplectic integrators and (iv) efficient solution of the Schrödinger equation. In the last decades several one step symplectic integrators have been produced based on symplectic geometry, (see the relevant literature and the references here). However, the study of multistep symplectic integrators is very poor. In this paper we investigate the closed Newton–Cotes formulae and we write them as symplectic multilayer structures. We also develop trigonometrically-fitted symplectic methods which are based on the closed Newton–Cotes formulae. We apply the symplectic schemes to the well known radial Schrödinger equation in order to investigate the efficiency of the proposed method to these type of problems.     

Symplectic Methods for the Numerical Solution of the Radial Shrödinger Equation

In this paper new symplectic-schemes for the numerical solution of the radial Shrödinger equation are proposed. In particular, symplectic integrators for Hamiltonian systems have been developed. Based on this approach, second- and third-order methods are proposed. These methods are more accurate than the existing ones. We compare these methods not only with the existing symplectic methods, but also with a classical Runge–Kutta–Nyström method.     

Event-driven Brownian dynamics for hard spheres

Brownian dynamics algorithms integrate Langevin equations numerically and allow to probe long time scales in simulations. A common requirement for such algorithms is that interactions in the system should vary little during an integration time step; therefore, computational efficiency worsens as the interactions become steeper. In the extreme case of hard-body interactions, standard numerical integrators become ill defined. Several approximate schemes have been invented to handle such cases, but little emphasis has been placed on testing the correctness of the integration scheme. Starting from the two-body Smoluchowski equation, the authors discuss a general method for the overdamped Brownian dynamics of hard spheres, recently developed by one of the authors. They test the accuracy of the algorithm and demonstrate its convergence for a number of analytically tractable test cases.     

New Second-order Exponentially and Trigonometrically Fitted Symplectic Integrators for the Numerical Solution of the Time-independent Schrödinger Equation

The solution of the one-dimensional time-independent Schr?dinger equation is considered by trigonometrically and exponentially fitted symplectic integrators. The Schr?dinger equation is first transformed into a Hamiltonian canonical equation. Numerical results are obtained for the one-dimensional harmonic oscillator, doubly anharmonic oscillator and the exponential potential.     

Ewald summation and multiple time step methods for molecular dynamics simulation of biological molecules

Methods by which to determine conditions for a molecular dynamics (MD) simulation of biological molecules were investigated. Derivation of the optimal parameters of the Ewald summation was described so as to give same precision to the real space, the reciprocal space summations and the van der Waals interaction. Later, the procedure by which to determine the condition of the multiple time step method by RESPA (REference System Propagator Algorithm; Tuckerman et al., 1992, J. Chem. Phys., 97, 1990) was described as exemplified by MD simulations of a solvated β-sheet peptide. The conservation of the total energy in a microcanonical ensemble was measured to investigate the stability of the simulation conditions. The most feasible respective combinations of the time steps were: 0.25 fs for bond, angle and torsion interactions; 2 fs for van der Waals interaction and Ewald real-space summation; and 4 fs for Ewald reciprocal-space summation. Though it retained an acceptable accuracy, this condition accelerated the simulation ten-fold compared to that in which a simple velocity-Verlet method with a time step of 0.25 fs was used. The update of the correction term due to excluded neighbors was then investigated. Better results were obtained when the correction was updated with the real-space than when it was updated with the reciprocal-space summation. Finally, an MD simulation as long as 50 ps performed under the optimal Ewald and RESPA parameters was thus determined. The trajectory showed a good stability, indicating the feasibility of the parameters.     

Symplectic Methods of Fifth Order for the Numerical Solution of the Radial Shrödinger Equation

In this paper the numerical solution of the radial Shrödinger equation via new proposed symplectic-schemes is investigated. In particular, the radial Schödinger equation is transformed into Hamiltonian canonical form and is solved via symplectic integrators. Based on this approach, fifth-order methods are proposed. We compare these methods with well-known existing symplectic methods. The numerical results show the efficiency of the proposed method.