Polymer dynamics in torsion space

The large scale motion of proteins, or covalently bonded polymers in general, is governed by the dynamics of the torsion angles, with bond lengths and bond angles kept approximately constant. In the present work, the Lagrangian equations of torsion motion are derived for a general macromolecule. The dynamics is implemented numerically for a test protein, using the velocity Verlet method as the integrator. The results indicate time steps of up to about 30 fs can be used for short time (up to at least 20 ps) simulations, before the dynamics and energy start to differ significantly from results obtained with smaller time steps. For longer time simulations, up to 1000 ps, a time step of 10 fs is relatively safe. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

LIN: A new algorithm to simulate the dynamics of biomolecules by combining implicit-integration and normal mode techniques

A central goal in molecular dynamics simulations is increasing the integration time-step to allow the capturing of biomolecular motion on biochemically interesting time frames. We previously made a step in that direction by developing the Langevin/implicit–Euler scheme. Here, we present a modified Langevin/implicit–Euler formulation for molecular dynamics. The new method still maintains the major advantage of the original scheme, namely, stability over a wide range of time-steps. However, it substantially reduces the damping effect of the high-frequency modes inherent in the original implicit scheme. The new formulation involves separation of the solution into two components, one of which is solved exactly using normal-mode techniques, the other of which is solved by implicit numerical integration. In this way, the high-frequency and fast-varying components are well resolved in the analytic solution component, while the remaining components of the motion are obtained by a large time-step integration phase. Full details of the new scheme are presented, accompanied by illustrative examples for a simple pendulum system. An application to liquid butane demonstrates stability of the simulations at time-steps up to 50 fs, still with activation of the high-frequency modes. © John Wiley & Sons, Inc.     

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.     

Molecular dynamics with the united-residue model of polypeptide chains. II. Langevin and Berendsen-bath dynamics and tests on model alpha-helical systems

The implementation of molecular dynamics (MD) with our physics-based protein united-residue (UNRES) force field, described in the accompanying paper, was extended to Langevin dynamics. The equations of motion are integrated by using a simplified stochastic velocity Verlet algorithm. To compare the results to those with all-atom simulations with implicit solvent in which no explicit stochastic and friction forces are present, we alternatively introduced the Berendsen thermostat. Test simulations on the Ala(10) polypeptide demonstrated that the average kinetic energy is stable with about a 5 fs time step. To determine the correspondence between the UNRES time step and the time step of all-atom molecular dynamics, all-atom simulations with the AMBER 99 force field and explicit solvent and also with implicit solvent taken into account within the framework of the generalized Born/surface area (GBSA) model were carried out on the unblocked Ala(10) polypeptide. We found that the UNRES time scale is 4 times longer than that of all-atom MD simulations because the degrees of freedom corresponding to the fastest motions in UNRES are averaged out. When the reduction of the computational cost for evaluation of the UNRES energy function is also taken into account, UNRES (with hydration included implicitly in the side chain-side chain interaction potential) offers about at least a 4000-fold speed up of computations relative to all-atom simulations with explicit solvent and at least a 65-fold speed up relative to all-atom simulations with implicit solvent. To carry out an initial full-blown test of the UNRES/MD approach, we ran Berendsen-bath and Langevin dynamics simulations of the 46-residue B-domain of staphylococcal protein A. We were able to determine the folding temperature at which all trajectories converged to nativelike structures with both approaches. For comparison, we carried out ab initio folding simulations of this protein at the AMBER 99/GBSA level. The average CPU time for folding protein A by UNRES molecular dynamics was 30 min with a single Alpha processor, compared to about 152 h for all-atom simulations with implicit solvent. It can be concluded that the UNRES/MD approach will enable us to carry out microsecond and, possibly, millisecond simulations of protein folding and, consequently, of the folding process of proteins in real time.     

Balancing an accurate representation of the molecular surface in generalized born formalisms with integrator stability in molecular dynamics simulations

Different integrator time steps in NVT and NVE simulations of protein and nucleic acid systems are tested with the GBMV (Generalized Born using Molecular Volume) and GBSW (Generalized Born with simple SWitching) methods. The simulation stability and energy conservation is investigated in relation to the agreement with the Poisson theory. It is found that very close agreement between generalized Born methods and the Poisson theory based on the commonly used sharp molecular surface definition results in energy drift and simulation artifacts in molecular dynamics simulation protocols with standard 2-fs time steps. New parameters are proposed for the GBMV method, which maintains very good agreement with the Poisson theory while providing energy conservation and stable simulations at time steps of 1 to 1.5 fs.     

Multiple time scale molecular dynamics for fluids with orientational degrees of freedom. II. Canonical and isokinetic ensembles

We have developed several multiple time stepping techniques to overcome the limitations on efficiency of molecular dynamics simulations of complex fluids. They include the modified canonical and isokinetic schemes, as well as the extended isokinetic Nosé-Hoover chain approach. The latter generalizes the method of Minary, Tuckerman, and Martyna for translational motion [Phys. Rev. Lett. 93, 150201 (2004)] to systems with both translational and orientational degrees of freedom. Although the microcanonical integrators are restricted to relatively small outer time steps of order of 16 fs, we show on the basis of molecular dynamics simulations of ambient water that in the canonical and isokinetic thermostats the size of these steps can be increased to 50 and 75 fs, respectively (at the same inner time step of 4 fs). Within the generalized isokinetic Nosé-Hoover chain algorithm we have derived, huge outer time steps of order of 500 fs can be used without losing numerical stability and affecting equilibrium properties.     

Evaluation of atomic pressure in the multiple time‐step integration algorithm

In molecular dynamics (MD) calculations, reduction in calculation time per MD loop is essential. A multiple time‐step (MTS) integration algorithm, the RESPA (Tuckerman and Berne, J. Chem. Phys. 1992, 97, 1990–2001), enables reductions in calculation time by decreasing the frequency of time‐consuming long‐range interaction calculations. However, the RESPA MTS algorithm involves uncertainties in evaluating the atomic interaction‐based pressure (i.e., atomic pressure) of systems with and without holonomic constraints. It is not clear which intermediate forces and constraint forces in the MTS integration procedure should be used to calculate the atomic pressure. In this article, we propose a series of equations to evaluate the atomic pressure in the RESPA MTS integration procedure on the basis of its equivalence to the Velocity‐Verlet integration procedure with a single time step (STS). The equations guarantee time‐reversibility even for the system with holonomic constrants. Furthermore, we generalize the equations to both (i) arbitrary number of inner time steps and (ii) arbitrary number of force components (RESPA levels). The atomic pressure calculated by our equations with the MTS integration shows excellent agreement with the reference value with the STS, whereas pressures calculated using the conventional ad hoc equations deviated from it. Our equations can be extended straightforwardly to the MTS integration algorithm for the isothermal NVT and isothermal–isobaric NPT ensembles. © 2017 Wiley Periodicals, Inc.     

Application of Molecular Dynamics Simulations in Molecular Property Prediction I: Density and Heat of Vaporization

Molecular mechanical force field (FF) methods are useful in studying condensed phase properties. They are complementary to experiment and can often go beyond experiment in atomic details. Even a FF is specific for studying structures, dynamics and functions of biomolecules, it is still important for the FF to accurately reproduce the experimental liquid properties of small molecules that represent the chemical moieties of biomolecules. Otherwise, the force field may not describe the structures and energies of macromolecules in aqueous solutions properly. In this work, we have carried out a systematic study to evaluate the General AMBER Force Field (GAFF) in studying densities and heats of vaporization for a large set of organic molecules that covers the most common chemical functional groups. The latest techniques, such as the particle mesh Ewald (PME) for calculating electrostatic energies, and Langevin dynamics for scaling temperatures, have been applied in the molecular dynamics (MD) simulations. For density, the average percent error (APE) of 71 organic compounds is 4.43% when compared to the experimental values. More encouragingly, the APE drops to 3.43% after the exclusion of two outliers and four other compounds for which the experimental densities have been measured with pressures higher than 1.0 atm. For heat of vaporization, several protocols have been investigated and the best one, P4/ntt0, achieves an average unsigned error (AUE) and a root-mean-square error (RMSE) of 0.93 and 1.20 kcal/mol, respectively. How to reduce the prediction errors through proper van der Waals (vdW) parameterization has been discussed. An encouraging finding in vdW parameterization is that both densities and heats of vaporization approach their "ideal" values in a synchronous fashion when vdW parameters are tuned. The following hydration free energy calculation using thermodynamic integration further justifies the vdW refinement. We conclude that simple vdW parameterization can significantly reduce the prediction errors. We believe that GAFF can greatly improve its performance in predicting liquid properties of organic molecules after a systematic vdW parameterization, which will be reported in a separate paper.     

Multiple time scale molecular dynamics for fluids with orientational degrees of freedom. I. Microcanonical ensemble

We propose a new approach to eliminate the resonance instabilities inherent in multiple time step molecular dynamics simulations. The approach is developed within the microcanonical ensemble on the basis of an energy-constrained technique in the presence of orientational degrees of freedom. While the single and standard multiscale methods are restricted to small time steps of 5 and 8 fs, respectively, it is shown in simulations of water that the algorithms we have derived postpone the appearance of the instabilities to larger steps of about 16 fs. Such steps are close to the upper theoretical limit of 20 fs peculiar to the microcanonical ensemble and can be used without affecting static and dynamical properties.     

New ways to boost molecular dynamics simulations

We describe a set of algorithms that allow to simulate dihydrofolate reductase (DHFR, a common benchmark) with the AMBER all‐atom force field at 160 nanoseconds/day on a single Intel Core i7 5960X CPU (no graphics processing unit (GPU), 23,786 atoms, particle mesh Ewald (PME), 8.0 Å cutoff, correct atom masses, reproducible trajectory, CPU with 3.6 GHz, no turbo boost, 8 AVX registers). The new features include a mixed multiple time‐step algorithm (reaching 5 fs), a tuned version of LINCS to constrain bond angles, the fusion of pair list creation and force calculation, pressure coupling with a “densostat,” and exploitation of new CPU instruction sets like AVX2. The impact of Intel's new transactional memory, atomic instructions, and sloppy pair lists is also analyzed. The algorithms map well to GPUs and can automatically handle most Protein Data Bank (PDB) files including ligands. An implementation is available as part of the YASARA molecular modeling and simulation program from www.YASARA.org . © 2015 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

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.     

Difficulties with multiple time stepping and fast multipole algorithm in molecular dynamics

Numerical experiments are performed on a 36,000-atom protein–DNA–water simulation to ascertain the effectiveness of two devices for reducing the time spent computing long-range electrostatics interactions. It is shown for Verlet-I/r-RESPA multiple time stepping, which is based on approximating long-range forces as widely separated impulses, that a long time step of 5 fs results in a dramatic energy drift and that this is reduced by using an even larger long time step. It is also shown that the use of as many as six terms in a fast multipole algorithm approximation to long-range electrostatics still fails to prevent significant energy drift even though four digits of accuracy is obtained. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1785–1791, 1997     

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.     

Optimization of parameters for molecular dynamics simulation using smooth particle-mesh Ewald in GROMACS 4.5

Based on our critique of requirements for performing an efficient molecular dynamics simulation with the particle-mesh Ewald (PME) implementation in GROMACS 4.5, we present a computational tool to enable the discovery of parameters that produce a given accuracy in the PME approximation of the full electrostatics. Calculations on two parallel computers with different processor and communication structures showed that a given accuracy can be attained over a range of parameter space, and that the attributes of the hardware and simulation system control which parameter sets are optimal. This information can be used to find the fastest available PME parameter sets that achieve a given accuracy. We hope that this tool will stimulate future work to assess the impact of the quality of the PME approximation on simulation outcomes, particularly with regard to the trade-off between cost and scientific reliability in biomolecular applications.     

Electron and nuclear dynamics of molecular clusters in ultraintense laser fields. II. Electron dynamics and outer ionization of the nanoplasma

We explore electron dynamics in molecular (CD4)(1061) clusters and elemental Xen (n=249-2171) clusters, responding to ultraintense (intensity I=10(16)-10(19) W cm(-2)) laser fields. Molecular dynamics simulations (including magnetic field and relativistic effects) and analyses of high-energy electron dynamics and nuclear ion dynamics in a cluster interacting with a Gaussian shaped laser field (frequency 0.35 fs(-1), photon energy 1.44 eV, phase 0, temporal width 25 fs) elucidated the time dependence of inner ionization, the formation of a nanoplasma of unbound electrons within the cluster or its vicinity, and of outer ionization. We determined the cluster size and the laser intensity dependence of these three sequential-parallel electronic processes. The characteristic times for cluster inner ionization (tau(ii)) and for outer ionization (tau(oi)) fall in the femtosecond time domain, i.e., tau(ii)=2-9 fs and tau(oi)=4-15 fs for (CD4)(1061), tau(ii)=7-30 fs and tau(oi)=5-13 fs for Xe(n) (n=479,1061), with both tau(ii) and tau(oi) decreasing with increasing I, in accord with the barrier suppression ionization mechanism for inner ionization of the constituents and the cluster barrier suppression ionization mechanism for outer ionization. The positive delay times Deltatau(OI) between outer and inner ionization (e.g., Deltatau(OI)=6.5 fs for Xen at I=10(16) W cm(-2) and Deltatau(OI)=0.2 fs for (CD4)(1061) at I=10(19) W cm(-2)) demonstrate that the outer/inner ionization processes are sequential. For (CD4)(1061), tau(ii)tau(oi), reflecting on the energetic hierarchy in the ionization of the Xe atoms. Quasiresonance contributions to the outer ionization of the nanoplasma were established, as manifested in the temporal oscillations in the inner/outer ionization levels, and in the center of mass of the nanoplasma electrons. The formation characteristics, dynamics, and response of the nanoplasma in molecular or elemental clusters were addressed. The nanoplasma is positively charged, with a high-average electron density [rho(P)=(2-3)10(22) cm(-3)], being characterized by high-average electron energies epsilon(av) (e.g., in Xe(1061) clusters epsilon(av)=54 eV at I=10(16) W cm(-2) and epsilon(av)=0.56-0.37 keV at I=10(18) W cm(-2), with epsilon(av) proportional, variant I(1/2)). Beyond the cluster boundary the average electron energy markedly increases, reaching electron energies in the range of 1.2-40 keV for outer ionization of Xe(n) (n=249-2171) clusters. The nanoplasma exhibits spatial inhomogeneity and angular anisotropy induced by the laser field. Femtosecond time scales are predicted for the nanoplasma production (rise times 7-3 fs), for the decay (decay times approximately 5 fs), and for the persistence time (30-10 fs) of a transient nanoplasma at I=10(17)-10(18) W cm(-2). At lower intensities of I=10(16) W cm(-2) a persistent nanoplasma with a "long" lifetime of > 50 fs will prevail.     

Fast lookup tables for interatomic interactions

We describe novel lookup tables for the rapid calculation of interatomic interactions. The tables have nonuniform distributions of bin widths tailored to minimize numerical error and maximize computational speed. Since interaction energies are precalculated, computer time requirements are essentially independent of the form of the potential function used. In test calculations using the AMBER force field and an internal coordinate Monte Carlo algorithm, the lookup table runs 15% faster than direct calculation of nonbonded interactions. The method is more advantageous for more complicated energy functions. As an example of a more complicated potential function, we have tested a pairwise approximation to accessible surface area. In this case, the use of the lookup table results in a speedup of a factor of two. The method is straightforward to implement and should be widely applicable. © 1996 by John Wiley & Sons, Inc.     

Stiffness and energy conservation in molecular dynamics: An improved integrator

Molecular dynamics is the integration of a set of coupled differential equations describing the motion of atoms over time. These equations exhibit the unfortunate property of stiffness, that is, terms of the equations (the forces on the atoms) are defined on several scales—ranging from tens of kcal/mol/Å to thousands of kcal/mol/Å. Additional nonconservative and stiff effects occur when a distance cutoff is used for the electrostatics and nonbonded potentials. Because the first derivative at the cutoff is essentially infinite, small variations in positions will cause large variations in energy and violate conservation of energy. The effects are demonstrated in a small system of 125 isolated water molecules. It is possible to greatly reduce and nearly eliminate the stiff integration effects with an improved integrator. The nonconservative effects of the distance cutoff cannot be removed by changing the integrator. © John Wiley & Sons, Inc.     

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.     

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.