Protecting High Energy Barriers: A New Equation to Regulate Boost Energy in Accelerated Molecular Dynamics Simulations

Molecular dynamics (MD) is one of the most common tools in computational chemistry. Recently, our group has employed accelerated molecular dynamics (aMD) to improve the conformational sampling over conventional molecular dynamics techniques. In the original aMD implementation, sampling is greatly improved by raising energy wells below a predefined energy level. Recently, our group presented an alternative aMD implementation where simulations are accelerated by lowering energy barriers of the potential energy surface. When coupled with thermodynamic integration simulations, this implementation showed very promising results. However, when applied to large systems, such as proteins, the simulation tends to be biased to high energy regions of the potential landscape. The reason for this behavior lies in the boost equation used since the highest energy barriers are dramatically more affected than the lower ones. To address this issue, in this work, we present a new boost equation that prevents oversampling of unfavorable high energy conformational states. The new boost potential provides not only better recovery of statistics throughout the simulation but also enhanced sampling of statistically relevant regions in explicit solvent MD simulations.     

Enhanced sampling method in molecular simulations using genetic algorithm for biomolecular systems

We propose a molecular simulation method using genetic algorithm (GA) for biomolecular systems to obtain ensemble averages efficiently. In this method, we incorporate the genetic crossover, which is one of the operations of GA, to any simulation method such as conventional molecular dynamics (MD), Monte Carlo, and other simulation methods. The genetic crossover proposes candidate conformations by exchanging parts of conformations of a target molecule between a pair of conformations during the simulation. If the candidate conformations are accepted, the simulation resumes from the accepted ones. While conventional simulations are based on local update of conformations, the genetic crossover introduces global update of conformations. As an example of the present approach, we incorporated genetic crossover to MD simulations. We tested the validity of the method by calculating ensemble averages and the sampling efficiency by using two kinds of peptides, ALA3 and (AAQAA)₃. The results show that for ALA3 system, the distribution probabilities of backbone dihedral angles are in good agreement with those of the conventional MD and replica-exchange MD simulations. In the case of (AAQAA)₃ system, our method showed lower structural correlation of α-helix structures than the other two methods and more flexibility in the backbone ψ angles than the conventional MD simulation. These results suggest that our method gives more efficient conformational sampling than conventional simulation methods based on local update of conformations. © 2018 Wiley Periodicals, Inc.     

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.     

Convergence and sampling efficiency in replica exchange simulations of peptide folding in explicit solvent

Replica exchange methods (REMs) are increasingly used to improve sampling in molecular dynamics (MD) simulations of biomolecular systems. However, despite having been shown to be very effective on model systems, the application of REM in complex systems such as for the simulation of protein and peptide folding in explicit solvent has not been objectively tested in detail. Here we present a comparison of conventional MD and temperature replica exchange MD (T-REMD) simulations of a beta-heptapeptide in explicit solvent. This system has previously been shown to undergo reversible folding on the time scales accessible to MD simulation and thus allows a direct one-to-one comparison of efficiency. The primary properties compared are the free energy of folding and the relative populations of different conformers as a function of temperature. It is found that to achieve a similar degree of precision T-REMD simulations starting from a random set of initial configurations were approximately an order of magnitude more computationally efficient than a single 800 ns conventional MD simulation for this system at the lowest temperature investigated (275 K). However, whereas it was found that T-REMD simulations are more than four times more efficient than multiple independent MD simulations at one temperature (300 K) the actual increase in conformation sampling was only twofold. The overall gain in efficiency using REMD resulted primarily from the ordering of different conformational states over temperature, as opposed to a large increase of conformational sampling. It is also shown that in this system exchanges are accepted primarily based on (random) fluctuations within the solvent and are not strongly correlated with the instantaneous peptide conformation raising questions in regard to the efficiency of T-REMD in larger systems.     

A comparative study of molecular dynamics in Cartesian and in internal coordinates: dynamical instability in the latter caused by nonlinearity of the equations of motion

The stability of a general molecular dynamics (MD) integration scheme is examined for simulations in generalized (internal plus external) coordinates (GCs). An analytic expression is derived for the local error in energy during each integration time step. This shows that the explicit dependence of the mass-matrix on GCs, which makes the system's Lagrange equations of motion nonlinear, causes MD simulations in GCs to be less stable than those in Cartesian coordinates (CCs). In terms of CCs, the corresponding mass-matrix depends only on atomic masses and thus atomistic motion is subject to the linear Newton equations, which makes the system more stable. Also investigated are two MD methods in GCs that utilize nonzero elements of the vibrational spectroscopic B-matrices. One updates positions and velocities in GCs that are iteratively adjusted so as to conform to the velocity Verlet equivalent in GCs. The other updates positions in GCs and velocities in CCs that are adjusted to satisfy the internal constraints of the new constrained WIGGLE MD scheme. The proposed methods are applied to an isolated n-octane molecule and their performances are compared with those of several CCMD schemes. The simulation results are found to be consistent with the analytic stability analysis. Finally, a method is presented for computing nonzero elements of B-matrices for external rotations without imposing the Casimir-Eckart conditions.     

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.     

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.     

Investigation of ligand exchange reactions in aqueous uranyl carbonate complexes using computational approaches

Carbonate anion exchange reactions with water in the uranyl-carbonate and calcium-uranyl-carbonate aqueous systems have been investigated using computational methods. Classical molecular dynamics (MD) simulations with the umbrella sampling technique were employed to determine potentials of mean force for the exchange reactions of water and carbonate. The presence of calcium counter-ions is predicted to increase the stability of the uranyl-carbonate species in accordance with previous experimental observations. However, the free energy barrier to carbonate exchange with water is found to be comparable both in the presence and absence of calcium cations. Possible implications of these results for uranyl adsorption on mineral surfaces are discussed. Density functional theory (DFT) calculations were also used to confirm the trends observed in classical molecular dynamics simulations and to corroborate the validity of the potential parameters employed in the MD scheme.     

FAMUSAMM: An algorithm for rapid evaluation of electrostatic interactions in molecular dynamics simulations

Within molecular dynamics simulations of protein–solvent systems the exact evaluation of long-range Coulomb interactions is computationally demanding and becomes prohibitive for large systems. Conventional truncation methods circumvent that computational problem, but are hampered by serious artifacts concerning structure and dynamics of the simulated systems. To avoid these artifacts we have developed an efficient and yet sufficiently accurate approximation scheme which combines the structure-adapted multipole method (SAMM) [C. Niedermeier and P. Tavan, J. Chem. Phys., 101 , 734 (1994)] with a multiple-time-step method. The computational effort for MD simulations required within our fast multiple-time-step structure-adapted multipole method (FAMUSAMM) scales linearly with the number of particles. For a system with 36,000 atoms we achieve a computational speed-up by a factor of 60 as compared with the exact evaluation of the Coulomb forces. Extended test simulations show that the applied approximations do not seriously affect structural or dynamical properties of the simulated systems. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1729–1749, 1997     

Simple algorithm for isothermal-isobaric molecular dynamics

We review principles of non-Hamiltonian statistical mechanics and present a new set of equations and integration algorithm for isothermal-isobaric dynamics. The chief advantage of the present scheme is that it is somewhat simpler than previous methods. We perform numerical simulations to test the accuracy of the algorithm and compare its stability to that of a "gold standard," a symplectic integrator for Hamiltonian dynamics of the same system. The stability of the isothermal-isobaric algorithm is comparable to the stability of the symplectic integrator.     

The performance of ANI-ML potentials for ligand-n(H2O) interaction energies and estimation of hydration free energies from end-point MD simulations

Here, we investigate the performance of “Accurate NeurAl networK engINe for Molecular Energies” (ANI), trained on small organic compounds, on bulk systems including non-covalent interactions and applicability to estimate solvation (hydration) free energies using the interaction between the ligand and explicit solvent (water) from single-step MD simulations. The method is adopted from ANI using the Atomic Simulation Environment (ASE) and predicts the non-covalent interaction energies at the accuracy of wb97x/6-31G(d) level by a simple linear scaling for the conformations sampled by molecular dynamics (MD) simulations of ligand-n(H₂O) systems. For the first time, we test ANI potentials' abilities to reproduce solvation free energies using linear interaction energy (LIE) formulism by modifying the original LIE equation. Our results on ~250 different complexes show that the method can be accurate and have a correlation of R² = 0.88–0.89 (MAE <1.0 kcal/mol) to the experimental solvation free energies, outperforming current end-state methods. Moreover, it is competitive to other conventional free energy methods such as FEP and BAR with 15-20 × fold reduced computational cost.     

Essential dynamics for the study of microstructures in liquids

Essential Dynamics (ED) is a powerful tool for analyzing molecular dynamics (MD) simulations and it is widely adopted for conformational analysis of large molecular systems such as, for example, proteins and nucleic acids. In this study, we extend the use of ED to the study of clusters of arbitrary size constituted by weakly interacting particles, for example, atomic clusters and supramolecular systems. The key feature of the method we present is the identification of the relevant atomic‐molecular clusters to be analyzed by ED for extracting the information of interest. The application of this computational approach allows a straightforward and unbiased conformational study of the local microstructures in liquids, as emerged from semiclassical MD simulations. The good performance of the method is demonstrated by calculating typical observables of liquid water, that is, NMR, NEXAFS O1s, and IR spectra, known to be rather sensitive both to the presence and to the conformational features of hydrogen‐bonded clusters. © 2014 Wiley Periodicals, Inc.     

Increasing the time step with mass scaling in Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations

Born‐Oppenheimer ab initio QM/MM molecular dynamics simulation with umbrella sampling is a state‐of‐the‐art approach to calculate free energy profiles of chemical reactions in complex systems. To further improve its computational efficiency, a mass‐scaling method with the increased time step in MD simulations has been explored and tested. It is found that by increasing the hydrogen mass to 10 amu, a time step of 3 fs can be employed in ab initio QM/MM MD simulations. In all our three test cases, including two solution reactions and one enzyme reaction, the resulted reaction free energy profiles with 3 fs time step and mass scaling are found to be in excellent agreement with the corresponding simulation results using 1 fs time step and the normal mass. These results indicate that for Born‐Oppenheimer ab initio QM/MM molecular dynamics simulations with umbrella sampling, the mass‐scaling method can significantly reduce its computational cost while has little effect on the calculated free energy profiles. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2009     

Hybrid particle‐field molecular dynamics simulations: Parallelization and benchmarks

The parallel implementation of a recently developed hybrid scheme for molecular dynamics (MD) simulations (Milano and Kawakatsu, J Chem Phys 2009, 130, 214106) where self‐consistent field theory (SCF) and particle models are combined is described. Because of the peculiar formulation of the hybrid method, considering single particles interacting with density fields, the most computationally expensive part of the hybrid particle‐field MD simulation can be efficiently parallelized using a straightforward particle decomposition algorithm. Benchmarks of simulations, including comparisons of serial MD and MD‐SCF program profiles, serial MD‐SCF and parallel MD‐SCF program profiles, and parallel benchmarks compared with efficient MD program GROMACS 4.5.4 are tested and reported. The results of benchmarks indicate that the proposed parallelization scheme is very efficient and opens the way to molecular simulations of large scale systems with reasonable computational costs. © 2012 Wiley Periodicals, Inc.     

Phase separation in binary mixtures containing polymers: A quantitative comparison of single-chain-in-mean-field simulations and computer simulations of the corresponding multichain systems

We discuss a phenomenological, coarse-grained simulation scheme, single-chain-in-mean-field (SCMF) simulation, for investigating the kinetics of phase separation in dense polymer blends and mixtures of polymers and solvents. In the spirit of self-consistent-field calculations, we approximate the interacting multichain problem by that of a single chain in an external field, which, in turn, depends on the local densities of the components. To study the time evolution of the mixture, we perform an explicit Monte Carlo (MC) simulation of an ensemble of independent chains in the external field and periodically calculate the average densities and update the external field. Unlike dynamic self-consistent-field theory, these SCMF simulations do not assume that the chain conformations relax much more quickly than the density and incorporate the single-chain dynamics explicitly rather than via an Onsager coefficient. This allows us to study systems with large spatial inhomogeneities and dynamic asymmetries. To assess the accuracy and limitations of the simulation scheme, we compare the results of SCMF simulations using a discretized Edwards Hamiltonian with computer simulations of the corresponding multichain system for (1) the early stages of spinodal decomposition of a symmetric binary polymer blend in response to a quench from χN = 0.314 to χN = 5 (where χ is the Flory–Huggins parameter and N is the number of segments), for which the growth rate of composition fluctuations is compared with MC simulations of the bond fluctuation model and alternative dynamic self-consistent-field calculations, and (2) the evaporation of a solvent from a low-molecular-weight thin polymer film, for which a comparison is made with molecular dynamics (MD) simulations of a bead-necklace model with a monomeric solvent. In the latter case, the polymer conformations are extracted from MD simulations and modeled in the SCMF simulations by a discretized Edwards Hamiltonian augmented by a chain-bending potential. From the MD simulations of thin polymer films in equilibrium with its vapor, phase coexistence has been determined, and the second- and third-order virial coefficients in the SCMF simulations have been adjusted accordingly. Finally, MD simulations of bulk solutions of a polymer and a solvent over a range of compositions, as well as the pure solvent at various densities, have been performed to determine self-diffusion coefficients that enter the SCMF simulations in the form of density-dependent segmental mobilities. A comparison of the polymer and solvent profiles in a thin film as a function of time and the fraction of the solvent evaporating from a solvent-swollen film, as obtained from MD simulations and parameterized SCMF simulations, shows satisfactory agreement for this simple mapping procedure. © 2005 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys 43: 934–958, 2005     

Atomic mobility in strained glassy polymers: The role of fold catastrophes on the potential energy surface

Deformation is known to enhance the atomic mobility in disordered systems such as polymer materials. To elucidate the origin of this phenomenon, we carry out two types of simulations: molecular dynamics (MD) simulations, which determine the atomic trajectories at finite temperature, and quasi-static simulations, which determine the atomic trajectories in the limit of zero temperature (and in the limit of zero shear rate). The quasi-static simulations show discontinuous changes in properties, such as system energy and atomic mobility. We use a normal mode analysis to show that these discontinuous changes arise from fold catastrophes of the potential energy landscape, in which energy minima flatten out and the heights of energy barriers decrease to zero; this was demonstrated by normal mode frequencies following a power law with an exponent of 0.5 as the discontinuous change is approached. After the fold catastrophe, the system relaxes to a different energy minimum, giving rise to atomic displacements. These fold catastrophes are the only mechanism for diffusive atomic displacements in the quasi-static simulations, where there is no thermal energy. We compared the mean-squared displacements as a function of strain from the quasi-static simulations to those from MD simulations (which do include thermal effects)—the similarity of the values of the mean-squared displacements in these two types of simulations demonstrates that the fold catastrophes underlie the enhanced dynamics in strained polymer systems even at finite temperature. © 2012 Wiley Periodicals, Inc. J Polym Sci Part B: Polym Phys, 2012     

Finite element interpolation for combined classical/quantum mechanical molecular dynamics simulations

A method is presented to interpolate the potential energy function for a part of a system consisting of a few degrees of freedom, such as a molecule in solution. The method is based on a modified finite element (FE) interpolation scheme. The aim is to save computer time when expensive methods such as quantum-chemical calculations are used to determine the potential energy function. The expensive calculations are only carried out if the molecule explores new unknown regions of the conformation space. If the molecule resides in regions previously explored, a cheap interpolation is performed instead of an expensive calculation, using known neighboring points. We report the interpolation techniques for the energies and the forces of the molecule, the handling of the FE mesh, and an application to a simple test example in molecular dynamics (MD) simulations. Good performance of the method was obtained (especially for MD simulations with a preceding Monte Carlo mesh generation) without losing accuracy. © 1997 John Wiley & Sons, Inc. J Comput Chem 18 : 1484–1495, 1997     

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.