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.     

Computationally efficient canonical molecular dynamics simulations by using a multiple time‐step integrator algorithm combined with the particle mesh Ewald method and with the fast multipole method

An efficient implementation of the canonical molecular dynamics simulation using the reversible reference system propagator algorithm (r‐RESPA) combined with the particle mesh Ewald method (PMEM) and with the macroscopic expansion of the fast multipole method (MEFMM) was examined. The performance of the calculations was evaluated for systems with 3000, 9999, 30,000, 60,000, and 99,840 particles. For a given accuracy, the optimal conditions for minimizing the CPU time for the implementation of the Ewald method, the PMEM, and the MEFMM were first analyzed. Using the optimal conditions, we evaluated the performance and the reliability of the integrated methods. For all the systems examined, the r‐RESPA with the PMEM was about twice as fast as the r‐RESPA with the MEFMM. The difference arose from the difference in the numerical complexities of the fast Fourier transform in the PMEM and from the transformation of the multipole moments into the coefficients of the local field expansion in the MEFMM. Compared with conventional methods, such as the velocity‐verlet algorithm with the Ewald method, significant speedups were obtained by the integrated methods; the speedup of the calculation was a function of system size, and was a factor of 100 for a system with 3000 particles and increased to a factor of 700 for a system with 99,840 particles. These integrated calculations are, therefore, promising for realizing large‐scale molecular dynamics simulations for complex systems. © 2000 John Wiley & Sons, Inc. J Comput Chem 21: 201–217, 2000     

Study on growth rate of TiO2 nanostructured thin films: simulation by molecular dynamics approach and modeling by artificial neural network

Effects of the deposition process parameters on the thickness of TiO₂ nanostructured film were simulated using the molecular dynamics (MD) approach and modeled by the artificial neural network (ANN) and regression method. Accordingly, TiO₂ nanostructured film was prepared experimentally with the sol–gel dip‐coating method. Structural instabilities can be expected, due to short‐ and/or long‐range intermolecular forces, leading to the surface inhomogeneities. In the MD simulation, the Morse potential function was used for the inter‐atomic interactions, and equations of motion for atoms were solved by Verlet algorithm. The effect of the withdrawal velocity, drying temperature and number of deposited layers were studied in order to characterize the film thickness. The results of MD simulations are reasonably consistent with atomic force microscopy, scanning electron microscopy and Dektak surface profiler. Finally, the outputs from experimental data were analyzed by using the ANN in order to investigate the effects of deposition process parameters on the film thickness. In this case, various architectures have been checked using 75% of experimental data for training of the ANN. Among the various architectures, feed‐forward back‐propagation network with trainer training algorithm was found as the best architecture. Based on the R‐squared value, the ANN is better than the regression model in predicting the film thickness. The statistical analysis for those results was then used to verify the fitness of the complex process model. Based on the results, this modeling methodology can explain the characteristics of the TiO₂ nanostructured thin film and growth mechanism varying with process conditions. © 2013 The Authors. Surface and Interface Analysis published by John Wiley & Sons Ltd.     

Coupled‐perturbed Hartree–Fock theory for quasi–one‐dimensional periodic systems: Calculation of static and dynamic nonlinear optical properties of model systems

An alternative method to solve the coupled‐perturbed Hartree–Fock (CPHF) equations for infinite quasi–one‐dimensional systems is presented. The new procedure follows a proposal made by Langhoff, Epstein, and Karplus to obtain perturbed wavefunctions free from arbitrary phase factors in each order of perturbation. It is based on the intermediate orthonormalization of the perturbed wavefunctions (which is different from the usual one) and a corresponding selection of the Lagrangian multipliers. In this way it is possible to incorporate the orthonormalization conditions into the set of CPHF equations. Moreover, a new, advantageous procedure to determine the derivatives of the wavefunction with respect to the quasimomentum k is presented. We report calculations of the dipole moment, the polarizability α, and the first hyperpolarizability β for different polymers (poly‐HF, poly‐H₂O, trans‐polyacetylene, polyyne, and polycarbonitrile) for different frequencies. These results are extensively compared with oligomer calculations. © 2003 Wiley Periodicals, Inc. Int J Quantum Chem 94: 251–268, 2003     

Single‐pass incremental force updates for adaptively restrained molecular dynamics

Adaptively restrained molecular dynamics (ARMD) allows users to perform more integration steps in wall‐clock time by switching on and off positional degrees of freedoms. This article presents new, single‐pass incremental force updates algorithms to efficiently simulate a system using ARMD. We assessed different algorithms for speedup measurements and implemented them in the LAMMPS MD package. We validated the single‐pass incremental force update algorithm on four different benchmarks using diverse pair potentials. The proposed algorithm allows us to perform simulation of a system faster than traditional MD in both NVE and NVT ensembles. Moreover, ARMD using the new single‐pass algorithm speeds up the convergence of observables in wall‐clock time. © 2017 Wiley Periodicals, Inc.     

Atomic forces for geometry‐dependent point multipole and Gaussian multipole models

In standard treatments of atomic multipole models, interaction energies, total molecular forces, and total molecular torques are given for multipolar interactions between rigid molecules. However, if the molecules are assumed to be flexible, two additional multipolar atomic forces arise because of (1) the transfer of torque between neighboring atoms and (2) the dependence of multipole moment on internal geometry (bond lengths, bond angles, etc.) for geometry‐dependent multipole models. In this study, atomic force expressions for geometry‐dependent multipoles are presented for use in simulations of flexible molecules. The atomic forces are derived by first proposing a new general expression for Wigner function derivatives . The force equations can be applied to electrostatic models based on atomic point multipoles or Gaussian multipole charge density. Hydrogen‐bonded dimers are used to test the intermolecular electrostatic energies and atomic forces calculated by geometry‐dependent multipoles fit to the ab initio electrostatic potential. The electrostatic energies and forces are compared with their reference ab initio values. It is shown that both static and geometry‐dependent multipole models are able to reproduce total molecular forces and torques with respect to ab initio, whereas geometry‐dependent multipoles are needed to reproduce ab initio atomic forces. The expressions for atomic force can be used in simulations of flexible molecules with atomic multipoles. In addition, the results presented in this work should lead to further development of next generation force fields composed of geometry‐dependent multipole models. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Molecular dynamics in the isothermal-isobaric ensemble: the requirement of a "shell" molecule. II. Simulation results

The results of a series of constant pressure and temperature molecular-dynamics (MD) simulation studies based on the rigorous shell particle formulation of the isothermal-isobaric (NpT) ensemble are presented. These MD simulations validate the newly proposed constant pressure equations of motion in which a "shell" particle is used to define uniquely the volume of the system [M. J. Uline and D. S. Corti, J. Chem. Phys. (to be published), preceding paper]. Ensemble averages obtained with the new MD NpT algorithm match the ensemble averages obtained using the previously derived shell particle Monte Carlo NpT method [D. S. Corti, Mol. Phys. 100, 1887 (2002)]. In addition, we also verify that the Hoover NpT MD algorithm [W. G. Hoover, Phys. Rev. A 31, 1695 (1985); 34, 2499 (1986)] generates the correct ensemble averages, though only when periodic boundary conditions are employed. The extension of the shell particle MD algorithm to multicomponent systems is also discussed, in which we show for equilibrium properties that the identity of the shell particle is completely arbitrary when periodic boundary conditions are applied. Self-diffusion coefficients determined with the shell particle equations of motion are also identical to those obtained in other ensembles. 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 effects of the shell particle on the dynamics of the system. Overall, the shell particle MD algorithm is an effective simulation method for studying systems exposed to a constant external pressure and may provide an advantage over other existing constant pressure approaches when developing nonequilibrium MD methods.     

Simulations of vibrational spectra from classical trajectories: calibration with ab initio force fields

An algorithm allowing simulating vibrational spectra from classical time-dependent trajectories was applied for infrared absorption, vibrational circular dichroism, Raman, and Raman optical activity of model harmonic systems. The implementation of the theory within the TINKER molecular dynamics (MD) program package was tested with ab initio harmonic force fields in order to determine the feasibility for more extended MD simulations. The results suggest that sufficiently accurate frequencies can be simulated with integration time steps shorter than about 0.5 fs. For a given integration time step, lower vibrational frequencies ( approximately 0-2000 cm(-1)) could be reproduced with a higher accuracy than higher-frequency vibrational modes (e.g., O-H and C-H stretching). In principle, the algorithm also provides correct intensities for ideal systems. In applied simulations, however, the intensity profiles are affected by an unrealistic energy distribution between normal modes and a slow energy relaxation. Additionally, the energy fluctuations may cause weakening of the intensities on average. For ab initio force fields, these obstacles could be overcome by an arbitrary normal mode energy correction. For general MD simulations, averaging of many shorter MD trajectories started with randomly distributed atomic velocities provided the best spectral shapes. alpha-pinene, D-gluconic acid, formaldehyde dimer, and the acetylprolineamide molecule were used in the tests.     

New algorithm for tensor contractions on multi‐core CPUs,GPUs, and accelerators enables CCSD and EOM‐CCSD calculations with over 1000 basis functions on a single compute node

A new hardware‐agnostic contraction algorithm for tensors of arbitrary symmetry and sparsity is presented. The algorithm is implemented as a stand‐alone open‐source code libxm . This code is also integrated with general tensor library libtensor and with the Q‐Chem quantum‐chemistry package. An overview of the algorithm, its implementation, and benchmarks are presented. Similarly to other tensor software, the algorithm exploits efficient matrix multiplication libraries and assumes that tensors are stored in a block‐tensor form. The distinguishing features of the algorithm are: (i) efficient repackaging of the individual blocks into large matrices and back, which affords efficient graphics processing unit (GPU)‐enabled calculations without modifications of higher‐level codes; (ii) fully asynchronous data transfer between disk storage and fast memory. The algorithm enables canonical all‐electron coupled‐cluster and equation‐of‐motion coupled‐cluster calculations with single and double substitutions (CCSD and EOM‐CCSD) with over 1000 basis functions on a single quad‐GPU machine. We show that the algorithm exhibits predicted theoretical scaling for canonical CCSD calculations, O (N ⁶), irrespective of the data size on disk. © 2017 Wiley Periodicals, Inc.     

Numerical Simulation of the Pressure Denaturation of a Helical β‐Peptide Heptamer Solvated in Methanol

The effect of elevated pressure on the conformational behavior of a β‐peptide heptamer ( 1 ) in MeOH solution was considered. The response of the peptide to elevated pressure was probed by means of molecular dynamics (MD) simulations, and described in atomic terms. The most‐striking features of the response are that the region of the ‘unfolded’ state of the peptide accessible at elevated pressure is narrow, and that thermal and pressure denaturation produce similar ‘unfolded’ states in the case of 1 .     

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.     

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.     

Improvement of methods for molecular dynamics integration

An application of symplectic implicit Runge–Kutta (RK ) integration schemes, the s-stage Gauss–Legendre Runge–Kutta (GLRK ) methods of order 2s, for the numerical solution of molecular dynamics (MD ) equation is described. The two-stage fourth-order GLRK method, the implicit midpoint rule, and the three-stage diagonally implicit RK method of order four are studied. The fixed-point iteraction was used for solving the resulting nonlinear system of equations. The algorithms were applied to a complex system of N particles interacting through a Lennard-Jones potential. The proposed symplectic methods for MD integration permit a wide range of time steps, are highly accurate and stable, and are thus suitable for the MD integration. © 1994 John Wiley & Sons, Inc.     

Reduced variable molecular dynamics

This article describes an extension to previously developed constraint techniques. These enhanced constraint methods will enable the study of large computational chemistry problems that cannot be easily handled with current constrained molecular dynamics (MD) methods. These methods are based on an O(N) solution to the constrained equations of motion. The benefits of this approach are that (1) the system constraints are solved exactly at each time step, (2) the solution algorithm is noniterative, (3) the algorithm is recursive and scales as O(N), (4) the algorithm is numerically stable, (5) the algorithm is highly amenable to parallel processing, and (6) potentially greater integration step sizes are possible. It is anticipated that application of this methodology will provide a 10- to 100-improvement in the speed of a large molecular trajectory as compared with the time required to run a conventional atomistic unconstrained simulation. It is, therefore, anticipated that this methodology will provide an enabling capacity for pursuing the drug discovery process for large molecular systems. © 1995 John Wiley & Sons, Inc.     

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.     

Reconstruction of atomistic details from coarse‐grained structures

We present an algorithm to reconstruct atomistic structures from their corresponding coarse‐grained (CG) representations and its implementation into the freely available molecular dynamics (MD) program package GROMACS. The central part of the algorithm is a simulated annealing MD simulation in which the CG and atomistic structures are coupled via restraints. A number of examples demonstrate the application of the reconstruction procedure to obtain low‐energy atomistic structural ensembles from their CG counterparts. We reconstructed individual molecules in vacuo (NCQ tripeptide, dipalmitoylphosphatidylcholine, and cholesterol), bulk water, and a WALP transmembrane peptide embedded in a solvated lipid bilayer. The first examples serve to optimize the parameters for the reconstruction procedure, whereas the latter examples illustrate the applicability to condensed‐phase biomolecular systems. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2010     

Exact and efficient calculation of Lagrange multipliers in biological polymers with constrained bond lengths and bond angles: proteins and nucleic acids as example cases

To accelerate molecular dynamics simulations, it is common to impose holonomic constraints on the hardest degrees of freedom. In this way, the time step used to integrate the equations of motion can be increased, thereby allowing longer total simulation times. The imposition of such constraints results in an aditional set of N_c equations (the equations of constraint) and unknowns (their associated Lagrange multipliers), whose solution is closely related to any algorithm implementing the constraints in Euclidean coordinates. In this work, it is shown that, due to the essentially linear structure of typical biological polymers the algebraic equations that need to be solved involve a matrix which is not only sparse, but also banded if the constraints are indexed in a skilful way. This allows the Lagrange multipliers to be obtained through a noniterative procedure, which can be considered exact up to machine precision, and which takes O(N_c) operations, instead of the usual O(N) for generic molecular systems. We develop the formalism, and describe the appropriate indexing for a number of model molecules. Finally, we provide a numerical example of the technique in a series of polyalanine peptides of different lengths. Although a use of the Lagrange multipliers without any modification in the solution of the underlying ordinary differential equations yields unstable integration algorithms, the central role of these quantities makes their efficient calculation useful for the improvement of methods that correctly enforce the exact satisfaction of the constraints at each time step. We provide several examples of this. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011