High-performance scalable molecular dynamics simulations of a polarizable force field based on classical Drude oscillators in NAMD

Incorporating the influence of induced polarization in large-scale atomistic molecular dynamics (MD) simulations is a critical challenge in the progress toward computations of increased accuracy. One computationally efficient treatment is based on the classical Drude oscillator in which an auxiliary charged particle is attached by a spring to each nucleus. Here, we report the first implementation of this model in the program NAMD. An extended Lagrangian dynamics with a dual-Langevin thermostat scheme applied to the Drude-nucleus pairs is employed to efficiently generate classical dynamic propagation near the self-consistent field limit. Large-scale MD simulations based on the Drude polarizable force field scale very well on massively distributed supercomputing platforms, the computational demand being only about 50-100% higher than for nonpolarizable models. As an illustration, a large-scale 150 mM NaCl aqueous salt solution is simulated, and the calculated ionic conductivity is shown to be in excellent agreement with experiment.     

Molecular dynamics simulations using the drude polarizable force field on GPUs with OpenMM: Implementation,validation, and benchmarks

Presented is the implementation of the Drude force field in the open‐source OpenMM simulation package allowing for access to graphical processing unit (GPU) hardware. In the Drude model, electronic degrees of freedom are represented by negatively charged particles attached to their parent atoms via harmonic springs, such that extra computational overhead comes from these additional particles and virtual sites representing lone pairs on electronegative atoms, as well as the associated thermostat and integration algorithms. This leads to an approximately fourfold increase in computational demand over additive force fields. However, by making the Drude model accessible to consumer‐grade desktop GPU hardware it will be possible to perform simulations of one microsecond or more in less than a month, indicating that the barrier to employ polarizable models has largely been removed such that polarizable simulations with the classical Drude model are readily accessible and practical.     

The structure of aqueous sodium hydroxide solutions: a combined solution x-ray diffraction and simulation study

To determine the structure of aqueous sodium hydroxide solutions, results obtained from x-ray diffraction and computer simulation (molecular dynamics and Car-Parrinello) have been compared. The capabilities and limitations of the methods in describing the solution structure are discussed. For the solutions studied, diffraction methods were found to perform very well in describing the hydration spheres of the sodium ion and yield structural information on the anion's hydration structure. Classical molecular dynamics simulations were not able to correctly describe the bulk structure of these solutions. However, Car-Parrinello simulation proved to be a suitable tool in the detailed interpretation of the hydration sphere of ions and bulk structure of solutions. The results of Car-Parrinello simulations were compared with the findings of diffraction experiments.     

Ab initio based polarizable force field parametrization

Experimental and simulation studies of anion-water systems have pointed out the importance of molecular polarization for many phenomena ranging from hydrogen-bond dynamics to water interfaces structure. The study of such systems at molecular level is usually made with classical molecular dynamics simulations. Structural and dynamical features are deeply influenced by molecular and ionic polarizability, which parametrization in classical force field has been an object of long-standing efforts. Although when classical models are compared to ab initio calculations at condensed phase, it is found that the water dipole moments are underestimated by approximately 30%, while the anion shows an overpolarization at short distances. A model for chloride-water polarizable interaction is parametrized here, making use of Car-Parrinello simulations at condensed phase. The results hint to an innovative approach in polarizable force fields development, based on ab initio simulations, which do not suffer for the mentioned drawbacks. The method is general and can be applied to the modeling of different systems ranging from biomolecular to solid state simulations.     

Hydration force between mica surfaces in aqueous KCl electrolyte solution

Liquid-vapor molecular dynamics simulations are performed to study the interaction forces between two mica surfaces in an aqueous KCl electrolyte solution. Strong repulsive hydration force is obtained within a distance of ~2 nm between the two mica surfaces, which cannot be explained by the continuum theory of double-layer repulsion. We find that this short-range repulsive hydration force is much stronger than the double-layer force between mica surfaces. Whereas the simulation system is much smaller than the surface force measurement system, fundamental mechanisms of repulsive hydration force are revealed. In particular, important features of the step-like force oscillatory behavior during normal compression and force hysteresis during retraction are observed. Detailed analysis of the ionic density distributions shows that the "forced adsorption" of diffusive K(+) ions onto mica surfaces during compression and the subsequent "slow desorption" of the absorbed K(+) ions from mica surfaces upon retraction are responsible for the hysteresis phenomenon. From a mechanics point of view, we attribute the load bearing capacity of the dense electrolyte to the very hard hydration shells of K(+) metal ions under confinement. We find that the hydrated K(+) ions and Cl(-) co-ions remain very diffusive in the aqueous film. Water molecules in the hydration layer are also very fluidic, in the sense that the diffusion constant of water molecules is less than its bulk value by at most 3 orders of magnitude under the extreme confinement.     

Prediction of nitroxide hyperfine coupling constants in solution from combined nanosecond scale simulations and quantum computations

We present a combined theoretical approach based on analyzing molecular dynamics trajectories (at the nanosecond scale) generated by use of classical polarizable force fields and on quantum calculations to compute averaged hyperfine coupling constants. That method is used to estimate the constant of a prototypical nitroxide: the dimethylnitroxide. The molecule is embedded during the simulations in a cubic box containing about 500 water molecules and the molecular dynamics is generated using periodic conditions. Once the trajectories are achieved, the nitroxide and its first hydration shell molecules are extracted, and the coupling constants are computed by considering the latter aggregates by means of quantum computations. However, all the water molecules of the bulk are also accounted for during those computations by means of the electrostatic potential fitted method. Our results exhibit that in order to predict accurate and reliable coupling constants, one needs to describe carefully the out-of-plane motion of the nitroxide nitrogen and to sample trajectories with a time interval of 400 fs at least to generate an uncorrelated large set of nitroxide structures. Compared to Car-Parrinello molecular dynamics techniques, our approach can be used readily to compute hyperfine coupling constants of large systems, such as nitroxides of great size interacting with macromolecules such as proteins or polymers.     

Understanding the dielectric properties of liquid amides from a polarizable force field

The role played by electronic polarization in the dielectric properties of liquid N-methyl acetamide (NMA) is examined using molecular dynamics simulations with a polarizable force field based on classical Drude oscillators. The model presented is the first force field shown to reproduce the anomalously large dielectric constant of liquid NMA. Details of the molecular polarizability are found to be important. For instance, all elements of the polarizability tensor, rather then just the trace, impact on the condensed phase properties. Two factors related to electronic polarizability are found to contribute to this large dielectric constant. First is the significant enhancement of the mean amide molecular dipole magnitude, which is 50% larger in the liquid than in the gas phase. Second is the consequent strong hydrogen bonding between molecular neighbors that enhances the orientational alignment of the molecular dipoles. Polarizable models of amide compounds that have two (acetamide) and zero (N,N-dimethyl acetamide) polar hydrogen-bond donor atoms are also investigated. Experimentally, the neat liquid dielectric constants at 373 K are 100 for NMA, 66 for acetamide and 26 for N,N-dimethyl acetamide. The polarizable models replicate this trend, predicting a dielectric constant of 92+/-5 for NMA, 66+/-3 for acetamide and 23+/-1 for N,N-dimethyl acetamide.     

Water Permeation Through DMPC Lipid Bilayers using Polarizable Charge Equilibration Force Fields

We investigate permeation energetics of water entering a model dimyristoylphosphatidylcholine (DMPC) bilayer via molecular dynamics simulations using polarizable Charge Equilibration (CHEQ) models. Potentials of mean force show 4.5-5.5 kcal/mol barriers for water permeation into bilayers. Barriers are highest when water coordination within the bilayer is prevented, and also when using force fields that accurately reproduce experimental alkane hydration free energies. The magnitude of the average water dipole moment decreases from 2.6 Debye (in bulk) to 1.88 Debye (in membrane interior). This variation correlates with the change in a water molecule's coordination number.     

Thermal and solvent effects on the NMR and UV parameters of some bioreductive drugs

15N NMR chemical shifts and n-->pi* electronic transition energy for metronidazole (1) has been calculated and compared with experimental data. A detailed computational study of 1 is presented, with special attention to the performance of various theoretical methods for reproducing spectroscopic parameters in solution. The most sophisticated approach involves density functional based on the Car-Parrinello molecular dynamics simulations of 1 in aqueous solution (BP86 level) and averaging chemical shifts and deltaE(n-->pi*) over snapshots from the trajectory. In the NMR and UV calculations for these snapshots (performed at the B3LYP level), a small number of discrete water molecules are retained, and the remaining bulk solution effects are included via a polarizable continuum model (PCM). A good agreement with experiment is also obtained using static geometry optimization and NMR computation of pristine 1 employing a PCM approach. Further theoretical predictions are also reported for 17O NMR and deltaE(n-->pi*) of three hydroxycinnamic acid derivatives, which suggest that it is essential to incorporate the dynamics and solvent effects for NMR and UV calculations in the condensed phase.     

Implementation of extended Lagrangian dynamics in GROMACS for polarizable simulations using the classical Drude oscillator model

Explicit treatment of electronic polarization in empirical force fields used for molecular dynamics simulations represents an important advancement in simulation methodology. A straightforward means of treating electronic polarization in these simulations is the inclusion of Drude oscillators, which are auxiliary, charge‐carrying particles bonded to the cores of atoms in the system. The additional degrees of freedom make these simulations more computationally expensive relative to simulations using traditional fixed‐charge (additive) force fields. Thus, efficient tools are needed for conducting these simulations. Here, we present the implementation of highly scalable algorithms in the GROMACS simulation package that allow for the simulation of polarizable systems using extended Lagrangian dynamics with a dual Nosé–Hoover thermostat as well as simulations using a full self‐consistent field treatment of polarization. The performance of systems of varying size is evaluated, showing that the present code parallelizes efficiently and is the fastest implementation of the extended Lagrangian methods currently available for simulations using the Drude polarizable force field.     

The ice-vapor interface and the melting point of ice I(h) for the polarizable POL3 water model

We use molecular dynamics simulations to determine the melting point of ice I(h) for the polarizable POL3 water force field (Dang, L. X. J. Chem. Phys.1992, 97, 2659). Simulations are performed on a slab of ice I(h) with two free surfaces at several different temperatures. The analysis of the time evolution of the total energy in the course of the simulations at the set of temperatures yields the melting point of the POL3 model to be T(m) = 180 ± 10 K. Moreover, the results of the simulations show that the degree of hydrogen-bond disorder occurring in the bulk of POL3 ice is larger (at the corresponding degree of undercooling) than in ice modeled by nonpolarizable water models. These results demonstrate that the POL3 water force field is rather a poor model for studying ice and ice-liquid or ice-vapor interfaces. While a number of polarizable water models have been developed over the past years, little is known about their performance in simulations of supercooled water and ice. This study thus highlights the need for testing of the existing polarizable water models over a broad range of temperatures, pressures, and phases, and developing a new polarizable water force field, reliable over larger areas of the phase diagram.     

Probing NMR parameters, structure and dynamics of 5-nitroimidazole derivatives. Density functional study of prototypical radiosensitizers

The 15N chemical shifts of metronidazole (1), secnidazole (2), nimorazole (3) and tinidazole (4), radiosensitizers based on the 5-nitroimidazole motif, are reported. A detailed computational study of 1 is presented, calling special attention to the performance of various theoretical methods in reproducing the 13C and 15N data observed in solution. The most sophisticated approach involves density functional-based Car-Parrinello molecular dynamics simulations (CPMD) of 1 in aqueous solution (BP86 level) and averaging chemical shifts over snapshots from the trajectory. In the NMR calculations for these snapshots (performed at the B3LYP level), a small number of discrete water molecules are retained, and the remaining bulk solution effects are included via a polarizable continuum model (PCM). A similarly good accord with experiment is obtained from much less involved, static geometry optimization and NMR computation of pristine 1 employing a PCM approach. Solvent effects on delta(15N), which are of the order of up to 20 ppm, are not due to changes in geometric parameters upon solvation, but arise from the direct response of the electronic wavefunction to the presence of the solvent, which can be represented by discrete molecules and/or the dielectric bulk.     

The onset of calcium carbonate nucleation: a density functional theory molecular dynamics and hybrid microsolvation/continuum study

Density functional theory (Perdew-Burke-Ernzerhof) based methods have been used to study the structure and hydration environment of the building blocks of CaCO 3 in aqueous solutions of calcium bicarbonate and calcium carbonate. Car-Parrinello molecular dynamics simulations of Ca(2+)/CO3(2-) and Ca (2+)/HCO3(-) in explicit water were performed to investigate the formation of CaCO3 and the hydration shell of the solvated hetero-ion pair. Our simulations show that the formation of the monomer of CaCO3 occurs with an associative mechanism and that the dominant building block of calcium (bi)carbonate in aqueous solution is Ca[eta(1)-(H)CO3](H2O)5, i.e., the preferred hydration number is five, while the (bi)carbonate is coordinated to the calcium in a monodentate mode. This result agrees with static calculations, where a hybrid approach using a combination of explicit solvent molecules and a polarizable continuum model has been applied to compute the solvation free energies of calcium bicarbonate species. Furthermore, the discrete-continuum calculations predict that the Ca(HCO3)2 and Ca(HCO3)3(-) species are stable in an aqueous environment preferentially as Ca(HCO3)2(H2O)4 and Ca(HCO3)3(H2O)2(-), respectively.     

Long-time correlations and hydrophobe-modified hydrogen-bonding dynamics in hydrophobic hydration

The physical mechanisms behind hydrophobic hydration have been debated for over 65 years. Spectroscopic techniques have the ability to probe the dynamics of water in increasing detail, but many fundamental issues remain controversial. We have performed systematic first-principles ab initio Car-Parrinello molecular dynamics simulations over a broad temperature range and provide a detailed microscopic view on the dynamics of hydration water around a hydrophobic molecule, tetramethylurea. Our simulations provide a unifying view and resolve some of the controversies concerning femtosecond-infrared, THz-GHz dielectric relaxation, and nuclear magnetic resonance experiments and classical molecular dynamics simulations. Our computational results are in good quantitative agreement with experiments, and we provide a physical picture of the long-debated "iceberg" model; we show that the slow, long-time component is present within the hydration shell and that molecular jumps and over-coordination play important roles. We show that the structure and dynamics of hydration water around an organic molecule are non-uniform.     

Assessing atomistic and coarse-grained force fields for protein-lipid interactions: the formidable challenge of an ionizable side chain in a membrane

Ionizable amino acid side chains play important roles in membrane protein structure and function, including the activation of voltage-gated ion channels, where it has been previously suggested that charged side chains may move through the hydrocarbon core of the membrane. However, all-atom molecular dynamics simulations have demonstrated large free energy barriers for such lipid-exposed motions. These simulations have also revealed that the membrane will deform due to the presence of a charged side chain, leading to a complex solvation microenvironment for which empirical force fields were not specifically parametrized. We have tested the ability of the all-atom CHARMM, Drude polarizable CHARMM, and a recent implementation of a coarse-grained force field to measure the thermodynamics of arginine-membrane interactions as a function of protonation state. We have employed model systems to attempt to match experimental bulk partitioning and quantum mechanical interactions within the membrane and found that free energy profiles from nonpolarizable and polarizable CHARMM simulations are accurate to within 1-2 kcal/mol. In contrast, the coarse-grained simulations failed to reproduce the same membrane deformations, exhibit interactions that are an order of magnitude too small, and thus, have incorrect free energy profiles. These results illustrate the need for careful parametrization of coarse-grained force fields and demonstrate the utility of atomistic molecular dynamics for providing quantitative thermodynamic and mechanistic analysis of protein-lipid interactions.     

Effective force fields for condensed phase systems from ab initio molecular dynamics simulation: a new method for force-matching

A novel least-squares fitting approach is presented to obtain classical force fields from trajectory and force databases produced by ab initio (e.g., Car-Parrinello) molecular dynamics (MD) simulations. The method was applied to derive effective nonpolarizable three-site force fields for liquid water at ambient conditions from Car-Parrinello MD simulations in the Becke-Lee-Yang-Parr approximation to the electronic density functional theory. The force-matching procedure includes a fit of short-ranged nonbonded forces, bonded forces, and atomic partial charges. The various parameterizations of the water force field differ by an enforced smooth cut-off applied to the short-ranged interaction term. These were obtained by fitting to the trajectory and force data produced by Car-Parrinello MD simulations of systems of 32 and 64 H(2)O molecules. The new water force fields were developed assuming both flexible or rigid molecular geometry. The simulated structural and self-diffusion properties of liquid water using the fitted force fields are in close agreement with those observed in the underlying Car-Parrinello MD simulations. The resulting empirical models compare to experiment much better than many conventional simple point charge (SPC) models. The fitted potential is also shown to combine well with more sophisticated intramolecular potentials. Importantly, the computational cost of the new models is comparable to that for SPC-like potentials.     

Effect of hydration on coordination properties of uranyl(VI) complexes. A first-principles molecular dynamics study

Results from Car-Parrinello molecular dynamics simulations are reported for [UO2(OH2)5]2+, UO2(NO3)2(OH2)2, and UO2(NO3)2(eta2-tmma) (tmma = tetramethylmalonamide) in the gas phase and in aqueous solution. The distances between uranyl and neutral ligands such as water and tmma are decreased by up to 0.2 angstroms upon hydration, whereas those between uranyl and the nitrate ion are increased by up to 0.08 angstroms. According to pointwise thermodynamic integration involving constrained molecular dynamics simulations, solvation facilitates the transition of the chelating nitrate ligand to a eta1-bonding mode: the free energy of UO2(eta2-NO3)(eta1-NO3)(OH2)2 relative to the bis-chelating minimum drops from 3.9 kcal/mol in vacuo to 1.4 kcal/mol in water. Optimizations in a polarizable continuum (specifically, the conductor-like screening model in conjunction with the zero-order regular approximation and triple-zeta Slater basis sets) can qualitatively reproduce the geometrical changes from explicit hydration.     

Solvation properties of N-acetyl-β-glucosamine: molecular dynamics study incorporating electrostatic polarization

N-Acetyl-β-glucosamine (NAG) is an important moiety of glycoproteins and is involved in many biological functions. However, conformational and dynamical properties of NAG molecules in aqueous solution, the most common biological environment, remain ambiguous due to limitations of experimental methods. Increasing efforts are made to probe structural properties of NAG and NAG-containing macromolecules, like peptidoglycans and polymeric chitin, at the atomic level using molecular dynamics simulations. In this work, we develop a polarizable carbohydrate force field for NAG and contrast simulation results of various properties using this novel force field and an analogous nonpolarizable (fixed charge) model. Aqueous solutions of NAG and its oligomers are investigated; we explore conformational properties (rotatable bond geometry), electrostatic properties (dipole moment distribution), dynamical properties (self-diffusion coefficient), hydrogen bonding (water bridge structure and dynamics), and free energy of hydration. The fixed-charge carbohydrate force field exhibits deviations from the gas phase relative rotation energy of exocyclic hydroxymethyl side chain and of chair/boat ring distortion. The polarizable force field predicts conformational properties in agreement with corresponding first-principles results. NAG-water hydrogen bonding pattern is studied through radial distribution functions (RDFs) and correlation functions. Intermolecular hydrogen bonding between solute and solvent is found to stabilize NAG solution structures while intramolecular hydrogen bonds define glycosidic linkage geometry of NAG oligomers. The electrostatic component of hydration free energy is highly dependent on force field atomic partial charges, influencing a more favorable free energy of hydration in the fixed-charge model compared to the polarizable model.     

Hydration structure of Y3+ and La3+ compared: an application of metadynamics

We studied the hydration structures of Y(3+) and La(3+) in aqueous solutions by applying the metadynamics method recently introduced as a tool to explore reaction pathways based on the Car-Parrinello molecular dynamics. By employing the metal-water oxygen coordination number as a collective variable of the metadynamics a couple of aqua and chloroaqua complexes are successfully generated within the time scale of 10 ps. The reconstructed free-energy surface captures the characteristics of the hydration of the light and heavy trivalent rare-earth ions. The present study demonstrates that the metadynamics based on the Car-Parrinello molecular dynamics is a promising tool for exploring the free-energy surface of complicated systems such as solutions.