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.     

Atomic multipoles: Electrostatic potential fit,local reference axis systems,and conformational dependence

Currently, all standard force fields for biomolecular simulations use point charges to model intermolecular electrostatic interactions. This is a fast and simple approach but has deficiencies when the electrostatic potential (ESP) is compared to that from ab initio methods. Here, we show how atomic multipoles can be rigorously implemented into common biomolecular force fields. For this, a comprehensive set of local reference axis systems is introduced, which represents a universal solution for treating atom‐centered multipoles for all small organic molecules and proteins. Furthermore, we introduce a new method for fitting atomic multipole moments to the quantum mechanically derived ESP. This methods yields a 50–90% error reduction compared to both point charges fit to the ESP and multipoles directly calculated from the ab initio electron density. It is shown that it is necessary to directly fit the multipole moments of conformational ensembles to the ESP. Ignoring the conformational dependence or averaging over parameters from different conformations dramatically deteriorates the results obtained with atomic multipole moments, rendering multipoles worse than partial charges. © 2012 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     

Force field dependence of phospholipid headgroup and acyl chain properties: Comparative molecular dynamics simulations of DMPC bilayers

The reliability of molecular simulations largely depends on the quality of the empirical force field parameters. Force fields used in lipid simulations continue to be improved to enhance the agreement with experiments for a number of different properties. In this work, we have carried out molecular dynamics simulations of neat DMPC bilayers using united‐atom Berger force field and three versions of all‐atom CHARMM force fields. Three different systems consisting of 48, 72, and 96 lipids were studied. Both particle mesh Ewald (PME) and spherical cut‐off schemes were used to evaluate the long‐range electrostatic interactions. In total, 21 simulations were carried out and analyzed to find out the dependence of lipid properties on force fields, system size, and schemes to calculate long‐range interactions. The acyl chain order parameters calculated from Berger and the recent versions of CHARMM simulations have shown generally good agreement with the experimental results. However, both sets of force fields deviate significantly from the experimentally observed P‐C dipolar coupling values for the carbon atoms that link the choline and glycerol groups with the phosphate groups. Significant differences are also observed in several headgroup parameters between CHARMM and Berger simulations. Our results demonstrate that when changes were introduced to improve CHARMM force field using PME scheme, all the headgroup parameters have not been reoptimized. The headgroup properties are likely to play a significant role in lipid–lipid, protein–lipid, and ligand–lipid interactions and hence headgroup parameters in phospholipids require refinement for both Berger and CHARMM force fields. © 2009 Wiley Periodicals, Inc.J Comput Chem, 2010     

An N log N approximation based on the natural organization of biomolecules for speeding up the computation of long range interactions

Presented here is a method, the hierarchical charge partitioning (HCP) approximation, for speeding up computation of pairwise electrostatic interactions in biomolecular systems. The approximation is based on multiple levels of natural partitioning of biomolecular structures into a hierarchical set of its constituent structural components. The charge distribution in each component is systematically approximated by a small number of point charges, which, for the highest level component, are much fewer than the number of atoms in the component. For short distances from the point of interest, the HCP uses the full set of atomic charges available. For long‐distance interactions, the approximate charge distributions with smaller sets of charges are used instead. For a structure consisting of N charges, the computational cost of computing the pairwise interactions via the HCP scales as O(N log N), under assumptions about the structural organization of biomolecular structures generally consistent with reality. A proof‐of‐concept implementation of the HCP shows that for large structures it can lead to speed‐up factors of up to several orders of magnitude relative to the exact pairwise O(N²) all‐atom computation used as a reference. For structures with more than 2000–3000 atoms the relative accuracy of the HCP (relative root‐mean‐square force error per atom), approaches the accuracy of the particle mesh Ewald (PME) method with parameter settings typical for biomolecular simulations. When averaged over a set of 600 representative biomolecular structures, the relative accuracies of the two methods are roughly equal. The HCP is also significantly more accurate than the spherical cutoff method. The HCP has been implemented in the freely available nucleic acids builder (NAB) molecular dynamics (MD) package in Amber tools. A 10 ns simulation of a small protein indicates that the HCP based MD simulation is stable, and that it can be faster than the spherical cutoff method. A critical benefit of the HCP approximation is that it is algorithmically very simple, and unlike the PME, the HCP is straightforward to use with implicit solvent models. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Effective multipoles and Yukawa electrostatics in dressed molecule theory

In this paper we derive the multipolar expansion of the screened Coulomb potential in electrolyte solutions with molecular solvent. The solute and solvent molecules can have arbitrary sizes, shapes, and internal charge distributions. We use the exact statistical mechanical definition of renormalized charge distributions coming from "dressed molecule theory" to determine the effective multipoles of a molecule immersed in an electrolyte. The effects of many-body correlations are fully included in our formally exact theory. We restrict ourselves to sufficiently dilute solutions so the screened Coulomb potential decays for large distances like a Yukawa function, exp(-kappa r)/r, where r is the distance and 1/kappa is the decay length (it is normally different from the Debye length). The resulting "Yukawa electrostatics" differ in many respects from ordinary, unscreened electrostatics. The "Yukawa charge" of a molecule (the lowest order moment in the multipolar expansion) is in general not equal to its Coulombic charge and it is not the integral of the renormalized charge distribution of the molecule. Moreover, as shown in this paper, the multipolar expansion of the Yukawa potential does not correspond, contrary to the case of the Coulomb potential, to its asymptotic expansion for large r. As a consequence, the charge term in the multipolar expansion is not the leading term in the asymptotic expansion. Instead, for large r values, multipoles of all orders contribute to the leading asymptotic term. Thus, the electrostatic potential from, for example, an electroneutral solvent molecule in an electrolyte solution has generally the same range as that from an ion. The proper asymptotic expansion for electrostatic interactions in electrolytes is derived. It is briefly shown how the multipole expansion formalism can also be applied in the Poisson-Boltzmann approximation for primitive model electrolytes.     

Force related atomic multipoles in planar molecules. Derivation of atomic quadrupole and octupole moments

Atomic multipoles as defined by current methods generally do not account for forces in molecules that arise from external electrostatic fields. It is pointed out that such forces and the electrostatic potential that the molecule itself generates are both determined by the molecular multipolar tensors. The latter constitute therefore the fundamental molecular constants that determine the molecular electrostatics apart from polarization. In general the multipolar tensors include contributions from the atomic multipoles and their fluxes. In planar molecules, however, the perpendicular charge flux is zero by symmetry. This gives rise to a (previously introduced) formalism that extracts analytical, force-related, atomic multipoles from the molecular multipolar tensors. This formalism is extended in this work to include force-related (FR) atomic quadrupoles and octupoles in planar molecules. The properties of the FR atomic multipoles, including their perpendicular fluxes, are discussed and some formal theoretical and computational advantages that characterize them are indicated. As an example, the electrostatics of OCS, including the molecular electrostatic potential and the forces on the nuclei due to an external point charge, is discussed.     

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

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

Scalable improvement of SPME multipolar electrostatics in anisotropic polarizable molecular mechanics using a general short‐range penetration correction up to quadrupoles

We propose a general coupling of the Smooth Particle Mesh Ewald SPME approach for distributed multipoles to a short‐range charge penetration correction modifying the charge‐charge, charge‐dipole and charge‐quadrupole energies. Such an approach significantly improves electrostatics when compared to ab initio values and has been calibrated on Symmetry‐Adapted Perturbation Theory reference data. Various neutral molecular dimers have been tested and results on the complexes of mono‐ and divalent cations with a water ligand are also provided. Transferability of the correction is adressed in the context of the implementation of the AMOEBA and SIBFA polarizable force fields in the TINKER‐HP software. As the choices of the multipolar distribution are discussed, conclusions are drawn for the future penetration‐corrected polarizable force fields highlighting the mandatory need of non‐spurious procedures for the obtention of well balanced and physically meaningful distributed moments. Finally, scalability and parallelism of the short‐range corrected SPME approach are addressed, demonstrating that the damping function is computationally affordable and accurate for molecular dynamics simulations of complex bio‐ or bioinorganic systems in periodic boundary conditions. © 2016 Wiley Periodicals, Inc.     

Implicit solvent models: Combining an analytical formulation of continuum electrostatics with simple models of the hydrophobic effect

The recent development of approximate analytical formulations of continuum electrostatics opens the possibility of efficient and accurate implicit solvent models for biomolecular simulations. One such formulation (ACE, Schaefer & Karplus, J. Phys. Chem., 1996, 100:1578) is used to compute the electrostatic contribution to solvation and conformational free energies of a set of small solutes and three proteins. Results are compared to finite-difference solutions of the Poisson equation (FDPB) and explicit solvent simulations and experimental data where available. Small molecule solvation free energies agree with FDPB within 1–1.5 kcal/mol, which is comparable to differences in FDPB due to different surface treatments or different force field parameterizations. Side chain conformation free energies of aspartate and asparagine are in qualitative agreement with explicit solvent simulations, while 74 conformations of a surface loop in the protein Ras are accurately ranked compared to FDPB. Preliminary results for solvation free energies of small alkane and polar solutes suggest that a recent Gaussian model could be used in combination with analytical continuum electrostatics to treat nonpolar interactions. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 322–335, 1999     

A convergent multipole expansion for 1,3 and 1,4 Coulomb interactions

Traditionally force fields express 1,3 and 1,4 interactions as bonded terms via potentials that involve valence and torsion angles, respectively. These interactions are not modeled by point charge terms, which are confined to electrostatic interactions between more distant atoms (1,n where n>4). Here we show that both 1,3 and 1,4 interactions can be described on the same footing as 1,n (n>4) interactions by a convergent multipole expansion of the Coulomb energy of the participating atom pairs. The atomic multipole moments are generated by the theory of quantum chemical topology. The procedure to make the multipole expansion convergent is based on a "shift procedure" described in earlier work [L. Joubert and P. L. A. Popelier, Molec. Phys. 100, 3357 (2002)].     

Fast evaluation of polarizable forces

Polarizability is considered to be the single most significant development in the next generation of force fields for biomolecular simulations. However, the self-consistent computation of induced atomic dipoles in a polarizable force field is expensive due to the cost of solving a large dense linear system at each step of a simulation. This article introduces methods that reduce the cost of computing the electrostatic energy and force of a polarizable model from about 7.5 times the cost of computing those of a nonpolarizable model to less than twice the cost. This is probably sufficient for the routine use of polarizable forces in biomolecular simulations. The reduction in computing time is achieved by an efficient implementation of the particle-mesh Ewald method, an accurate and robust predictor based on least-squares fitting, and non-stationary iterative methods whose fast convergence is accelerated by a simple preconditioner. Furthermore, with these methods, the self-consistent approach with a larger timestep is shown to be faster than the extended Lagrangian approach. The use of dipole moments from previous timesteps to calculate an accurate initial guess for iterative methods leads to an energy drift, which can be made acceptably small. The use of a zero initial guess does not lead to perceptible energy drift if a reasonably strict convergence criterion for the iteration is imposed.     

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     

Long-range Lennard-Jones and electrostatic interactions in interfaces: application of the isotropic periodic sum method

Molecular dynamics (MD) simulations of heptane/vapor, hexadecane/vapor, water/vapor, hexadecane/water, and dipalmitoylphosphatidylcholine (DPPC) bilayers and monolayers are analyzed to determine the accuracy of treating long-range interactions in interfaces with the isotropic periodic sum (IPS) method. The method and cutoff (rc) dependences of surface tensions, density profiles, water dipole orientation, and electrostatic potential profiles are used as metrics. The water/vapor, heptane/vapor, and hexadecane/vapor interfaces are accurately and efficiently calculated with 2D IPS (rc=10 A). It is demonstrated that 3D IPS is not practical for any of the interfacial systems studied. However, the hybrid method PME/IPS [Particle Mesh Ewald for electrostatics and 3D IPS for Lennard-Jones (LJ) interactions] provides an efficient way to include both types of long-range forces in simulations of large liquid/vacuum and all liquid/liquid interfaces, including lipid monolayers and bilayers. A previously published pressure-based long-range LJ correction yields results similar to those of PME/IPS for liquid/liquid interfaces. The contributions to surface tension of LJ terms arising from interactions beyond 10 A range from 13 dyn/cm for the hexadecane/vapor interface to approximately 3 dyn/cm for hexadecane/water and DPPC bilayers and monolayers. Surface tensions of alkane/vapor, hexadecane/water, and DPPC monolayers based on the CHARMM lipid force fields agree very well with experiment, whereas surface tensions of the TIP3P and TIP4P-Ew water models underestimate experiment by 16 and 11 dyn/cm, respectively. Dipole potential drops (DeltaPsi) are less sensitive to long-range LJ interactions than surface tensions. However, DeltaPsi for the DPPC bilayer (845+/-3 mV proceeding from water to lipid) and water (547+/-2 mV for TIP4P-Ew and 521+/-3 mV for TIP3P) overestimate experiment by factors of 3 and 5, respectively, and represent expected deficiencies in nonpolarizable force fields.     

Simulation of Ca2+ and Mg2+ solvation using polarizable atomic multipole potential

The alkaline earth metals calcium and magnesium are critically involved in many biomolecular processes. To understand the hydration thermodynamics of these ions, we have performed molecular dynamics simulations using a polarizable potential. Particle-mesh Ewald for point multipoles has been applied to the calculation of electrostatic interactions. The parameters in this model have been determined from an ab initio quantum mechanical calculation of dimer interactions between ions and water. Two methods for ion solvation free energy calculation, free energy perturbation, and the Bennett acceptance ratio have been compared. Both predict results consistent with other theoretical estimations while the Bennett approach leads to a much smaller statistical error. Based on the Born theory and the ion-oxygen radial distribution functions, we estimate the effective size of the ions in solution, concluding that K(+) > Na(+) congruent with Ca(2+) > Mg(2+). There appears to be much stronger perturbation in water structure, dynamics, and dipole moment around the divalent cations than the monovalent K(+) and Na(+). The average water coordination numbers for Ca(2+) and Mg(2+) are 7.3 and 6, respectively. The lifetime of water molecules in the first solvation shell of Mg(2+) is on the order of hundreds of picoseconds, in contrast to only few picoseconds for Ca(2+), K(+), or Na(+).