A fast path integral method for polarizable force fields

Path Integral Simulations of Condensed-Phase Vibrational Spectroscopy

Recent theoretical and algorithmic developments have improved the accuracy with which path integral dynamics methods can include nuclear quantum effects in simulations of condensed-phase vibrational spectra. Such methods are now understood to be approximations to the delocalized classical Matsubara dynamics of smooth Feynman paths, which dominate the dynamics of systems such as liquid water at room temperature. Focusing mainly on simulations of liquid water and hexagonal ice, we explain how the recently developed quasicentroid molecular dynamics (QCMD), fast-QCMD, and temperature-elevated path integral coarse-graining simulations (Te PIGS) methods generate classical dynamics on potentials of mean force obtained by averaging over quantum thermal fluctuations. These new methods give very close agreement with one another, and the Te PIGS method has recently yielded excellent agreement with experimentally measured vibrational spectra for liquid water, ice, and the liquid-air interface. We also discuss the limitations of such methods.

Developing machine-learned potentials to simultaneously capture the dynamics of excess protons and hydroxide ions in classical and path integral simulations

The transport of excess protons and hydroxide ions in water underlies numerous important chemical and biological processes. Accurately simulating the associated transport mechanisms ideally requires utilizing ab initio molecular dynamics simulations to model the bond breaking and formation involved in proton transfer and path-integral simulations to model the nuclear quantum effects relevant to light hydrogen atoms. These requirements result in a prohibitive computational cost, especially at the time and length scales needed to converge proton transport properties. Here, we present machine-learned potentials (MLPs) that can model both excess protons and hydroxide ions at the generalized gradient approximation and hybrid density functional theory levels of accuracy and use them to perform multiple nanoseconds of both classical and path-integral proton defect simulations at a fraction of the cost of the corresponding ab initio simulations. We show that the MLPs are able to reproduce ab initio trends and converge properties such as the diffusion coefficients of both excess protons and hydroxide ions. We use our multi-nanosecond simulations, which allow us to monitor large numbers of proton transfer events, to analyze the role of hypercoordination in the transport mechanism of the hydroxide ion and provide further evidence for the asymmetry in diffusion between excess protons and hydroxide ions.

Routine Molecular Dynamics Simulations Including Nuclear Quantum Effects: From Force Fields to Machine Learning Potentials

We report the implementation of a multi-CPU and multi-GPU massively parallel platform dedicated to the explicit inclusion of nuclear quantum effects (NQEs) in the Tinker-HP molecular dynamics (MD) package. The platform, denoted Quantum-HP, exploits two simulation strategies: the Ring-Polymer Molecular Dynamics (RPMD) that provides exact structural properties at the cost of a MD simulation in an extended space of multiple replicas and the adaptive Quantum Thermal Bath (adQTB) that imposes the quantum distribution of energy on a classical system via a generalized Langevin thermostat and provides computationally affordable and accurate (though approximate) NQEs. We discuss some implementation details, efficient numerical schemes, and parallelization strategies and quickly review the GPU acceleration of our code. Our implementation allows an efficient inclusion of NQEs in MD simulations for very large systems, as demonstrated by scaling tests on water boxes with more than 200,000 atoms (simulated using the AMOEBA polarizable force field). We test the compatibility of the approach with Tinker-HP's recently introduced Deep-HP machine learning potentials module by computing water properties using the DeePMD potential with adQTB thermostatting. Finally, we show that the platform is also compatible with the alchemical free energy estimation capabilities of Tinker-HP and fast enough to perform simulations. Therefore, we study how NQEs affect the hydration free energy of small molecules solvated with the recently developed Q-AMOEBA water force field. Overall, the Quantum-HP platform allows users to perform routine quantum MD simulations of large condensed-phase systems and will help to shed new light on the quantum nature of important interactions in biological matter.

Improving Condensed-Phase Water Dynamics with Explicit Nuclear Quantum Effects: The Polarizable Q-AMOEBA Force Field

We introduce a new parametrization of the AMOEBA polarizable force field for water denoted Q-AMOEBA, for use in simulations that explicitly account for nuclear quantum effects (NQEs). This study is made possible thanks to the recently introduced adaptive Quantum Thermal Bath (adQTB) simulation technique which computational cost is comparable to classical molecular dynamics. The flexible Q-AMOEBA model conserves the initial AMOEBA functional form, with an intermolecular potential including an atomic multipole description of electrostatic interactions (up to quadrupole), a polarization contribution based on the Thole interaction model and a buffered 14-7 potential to model van der Waals interactions. It has been obtained by using a ForceBalance fitting strategy including high-level quantum chemistry reference energies and selected condensed-phase properties targets. The final Q-AMOEBA model is shown to accurately reproduce both gas-phase and condensed-phase properties, notably improving the original AMOEBA water model. This development allows the fine study of NQEs on water liquid phase properties such as the average H-O-H angle compared to its gas-phase equilibrium value, isotope effects, and so on. Q-AMOEBA also provides improved infrared spectroscopy prediction capabilities compared to AMOEBA03. Overall, we show that the impact of NQEs depends on the underlying model functional form and on the associated strength of hydrogen bonds. Since adQTB simulations can be performed at near classical computational cost using the Tinker-HP package, Q-AMOEBA can be extended to organic molecules, proteins, and nucleic acids opening the possibility for the large-scale study of the importance of NQEs in biophysics.

Generalized Moment Correction for Long-Ranged Electrostatics

Describing long-ranged electrostatics using short-ranged pair potentials is appealing because the computational complexity scales linearly with the number of particles. The foundation of the approach presented here is to mimic the long-ranged medium response by cancelling electric multipoles within a small cutoff sphere. We propose a rigorous and formally exact new method that cancels up to infinitely many multipole moments and is free of operational damping parameters often required in existing theories. Using molecular dynamics simulations of water with and without added salt, we discuss radial distribution functions, Kirkwood–Buff integrals, dielectrics, diffusion coefficients, and angular correlations in relation to existing electrostatic models. We find that the proposed method is an efficient and accurate alternative for handling long-ranged electrostatics as compared to Ewald summation schemes. The methodology and proposed parameterization are applicable also for dipole–dipole interactions.

On short-ranged pair-potentials for long-range electrostatics

Fast and accurate summation of long-range electrostatics by using a short-ranged pair-potential that ensures moment cancellation in the cutoff sphere.

Using a monomer potential energy surface to perform approximate path integral molecular dynamics simulation of ab initio water at near-zero added cost

It is now established that nuclear quantum motion plays an important role in determining water's hydrogen bonding, structure, and dynamics. Such effects are important to include in density functional theory (DFT) based molecular dynamics simulation of water. The standard way of treating nuclear quantum effects, path integral molecular dynamics (PIMD), multiplies the number of energy/force calculations by the number of beads required. In this work we introduce a method whereby PIMD can be incorporated into a DFT simulation with little extra cost and little loss in accuracy. The method is based on the many body expansion of the energy and has the benefit of including a monomer level correction to the DFT energy. Our method calculates intramolecular forces using the highly accurate monomer potential energy surface developed by Partridge-Schwenke, which is cheap to evaluate. Intermolecular forces and energies are calculated with DFT only once per timestep using the centroid positions. We show how our method may be used in conjunction with a multiple time step algorithm for an additional speedup and how it relates to ring polymer contraction and other schemes that have been introduced recently to speed up PIMD simulations. We show that our method, which we call "monomer PIMD", correctly captures changes in the structure of water found in a full PIMD simulation but at much lower computational cost.

Quantum Dynamics and Spectroscopy of Ab Initio Liquid Water: The Interplay of Nuclear and Electronic Quantum Effects

Understanding the reactivity and spectroscopy of aqueous solutions at the atomistic level is crucial for the elucidation and design of chemical processes. However, the simulation of these systems requires addressing the formidable challenges of treating the quantum nature of both the electrons and nuclei. Exploiting our recently developed methods that provide acceleration by up to 2 orders of magnitude, we combine path integral simulations with on-the-fly evaluation of the electronic structure at the hybrid density functional theory level to capture the interplay between nuclear quantum effects and the electronic surface. Here we show that this combination provides accurate structure and dynamics, including the full infrared and Raman spectra of liquid water. This allows us to demonstrate and explain the failings of lower-level density functionals for dynamics and vibrational spectroscopy when the nuclei are treated quantum mechanically. These insights thus provide a foundation for the reliable investigation of spectroscopy and reactivity in aqueous environments.

Path integral-GC-AdResS simulation of a large hydrophobic solute in water: a tool to investigate the interplay between local microscopic structures and quantum delocalization of atoms in space

We perform large scale quantum (path integral) molecular dynamics simulations of a C₆₀ -like molecule in water.

Force-field functor theory: classical force-fields which reproduce equilibrium quantum distributions

Feynman and Hibbs were the first to variationally determine an effective potential whose associated classical canonical ensemble approximates the exact quantum partition function. We examine the existence of a map between the local potential and an effective classical potential which matches the exact quantum equilibrium density and partition function. The usefulness of such a mapping rests in its ability to readily improve Born-Oppenheimer potentials for use with classical sampling. We show that such a map is unique and must exist. To explore the feasibility of using this result to improve classical molecular mechanics, we numerically produce a map from a library of randomly generated one-dimensional potential/effective potential pairs then evaluate its performance on independent test problems. We also apply the map to simulate liquid para-hydrogen, finding that the resulting radial pair distribution functions agree well with path integral Monte Carlo simulations. The surprising accessibility and transferability of the technique suggest a quantitative route to adapting Born-Oppenheimer potentials, with a motivation similar in spirit to the powerful ideas and approximations of density functional theory.

Systematic Parametrization of Polarizable Force Fields from Quantum Chemistry Data

We introduce ForceBalance, a method and free software package for systematic force field optimization with the ability to parametrize a wide variety of functional forms using flexible combinations of reference data. We outline several important challenges in force field development and how they are addressed in ForceBalance, and present an example calculation where these methods are applied to develop a highly accurate polarizable water model. ForceBalance is available for free download at https://simtk.org/home/forcebalance.

Isotope effects in water as investigated by neutron diffraction and path integral molecular dynamics

The structures of heavy and light water at 300 K were investigated by using a joint approach in which the method of neutron diffraction with oxygen isotope substitution was complemented by path integral molecular dynamics simulations. The diffraction results, which give intra-molecular O-D and O-H bond distances of 0.985(5) and 0.990(5) Å, were found to be in best agreement with those obtained by using the flexible anharmonic TTM3-F water model. Both techniques show a difference of  ≃ 0.5% between the O-D and O-H intra-molecular bond lengths, and the results support a competing quantum effects model for water in which its structural and dynamical properties are governed by an offset between intra-molecular and inter-molecular quantum contributions. Further consideration of the O-O correlations is needed in order to improve agreement with experiment.

Oxygen as a site specific probe of the structure of water and oxide materials

The method of oxygen isotope substitution in neutron diffraction is introduced as a site specific structural probe. It is employed to measure the structure of light versus heavy water, thus circumventing the assumption of isomorphism between H and D as used in more traditional neutron diffraction methods. The intramolecular and intermolecular O-H and O-D pair correlations are in excellent agreement with path integral molecular dynamics simulations, both techniques showing a difference of ≃0.5% between the O-H and O-D intramolecular bond distances. The results support the validity of a competing quantum effects model for water in which its structural and dynamical properties are governed by an offset between intramolecular and intermolecular quantum contributions.

The polarizable point dipoles method with electrostatic damping: implementation on a model system

Recently, the use of polarizable force fields in Molecular Dynamics simulations has been gaining importance, since they allow a better description of heterogeneous systems compared to simple point charges force fields. Among the various techniques developed in the last years the one based on polarizable point dipoles represents one of the most used. In this paper, we review the basic technical issues of the method, illustrating the way to implement intramolecular and intermolecular damping of the electrostatic interactions, either with and without the Ewald summation method. We also show how to reduce the computational overhead for evaluating the dipoles, introducing to the state-of-the-art methods: the extended Lagrangian method and the always stable predictor corrector method. Finally we discuss the importance of screening the electrostatic interactions at short range, defending this technique against simpler approximations usually made. We compare results of density functional theory and classical force field-based Molecular Dynamics simulations of chloride in water.

Properties of liquid water from a systematic refinement of a high-rank multipolar electrostatic potential

We build on previous work [S. Y. Liem and P. L. A. Popelier, J. Chem. Theory Comput. 4, 353 (2008)], where for the first time, a high-rank multipolar electrostatic potential was used in molecular dynamics simulations of liquid water at a wide range of pressures and temperatures, and using a multipolar Ewald summation. Water is represented as a rigid body, with atomic multipole moments defined by quantum chemical topology partitioning its gas phase electron density. The effect of the level of theory on the local structure of liquid water is systematically addressed. Values for Lennard-Jones (LJ) parameters are optimized, for both oxygen and hydrogen atoms, against bulk properties. The best LJ parameters were then used in a set of simulations at 30 different temperatures (1 atm) and another set at 11 different pressures (at 298 K). Inclusion of the hydrogen LJ parameters significantly increases the self-diffusion coefficient. The behavior of bulk properties was studied and the local water structure analyzed by both radial and spatial distribution functions. Comparisons with familiar point-charge potentials, such as TIP3P, TIP4P, TIP5P, and simple point charge, show the benefits of multipole moments.