The quest for the best nonpolarizable water model from the adaptive force matching method

Simulation of Linear and Cyclic Alkanes with Second-Order Møller-Plesset Perturbation Theory through Adaptive Force Matching

Predicting ensemble properties, such as density and heat of vaporization, of small hydrocarbons is challenging due to the dispersion-dominated weak interactions between these molecules. With the adaptive force matching (AFM) method, the bonded and short-range nonbonded interactions are fitted to second-order Møller-Plesset perturbation theory (MP2) references computed with the def2-TZVP basis set. The dispersion is modeled using symmetry adapted perturbation theory (SAPT) at MP4 accuracy using the def2-TZVPD basis set. A new charge matrix decomposition technique is described to obtain partial charges in AFM. Although the models developed do not have any empirical parameters, several properties of the resulting models are compared with experiments as validations. The density, heat of vaporization, pressure dependence of density, diffusion constants, and surface tensions all show quantitative agreement with experiments. Although the density shows a very small systematic error, which could be due to missing three-body dispersion, the heat of vaporization agrees with experiments of within 0.5%. The paper shows that AFM can be used as a reliable tool to enable simulations at post-Hartree-Fock quality at the cost of molecular mechanics force fields.

Water Potential from Adaptive Force Matching for Ice and Liquid with Revised Dispersion Predicts Supercooled Liquid Anomalies in Good Agreement with Two Independent Experimental Fits

A revised version of the Water potential from Adaptive force matching for Ice and Liquid (WAIL) was developed by using the previous data set for fitting the WAIL model but with a dispersion term calculated using symmetry adapted perturbation theory (SAPT). The model has no adjustable parameters and relies solely on fitting first-principles information. The new model, named revised WAIL (rWAIL), shows improved predictions of most properties of water when compared to the previously published WAIL model. The rWAIL model also compares favorably to other first-principles-derived water models, such as MB-Pol, at only a fraction of the computational cost. The rWAIL model is used to study the properties of supercooled water. The model shows evidence of a liquid-liquid phase transition (LLPT) in the supercooled regimes with the liquid-liquid critical point (LLCP) at 203 K and 90 MPa. This estimate is in good agreement with a recent polynomial fit to the experimental density of water. Also, the fit to the surface tension of supercooled water based on the rWAIL model shows excellent agreement with the corresponding fit to the experimental data. Consistent with previously published molecular dynamics and experimental data, the surface tension of water exhibits exponential growth in the supercooled regime, which is likely a result of the emergence of a low-density liquid form of water. The simulation thus unites two separate experimental fits with one first-principles-based model, lending strong evidence of an LLPT in real water.

Dipole Cooperativity and Polarization Frustration Determine the Secondary Structure Distribution of Short Alanine Peptides in Water

The physical driving forces for secondary structure preferences of hydrated alanine peptide are investigated with B3LYP-D3(BJ) and the adaptive force matching (AFM) method. The AFM fit to the DFT surface, ALA2022, provides excellent agreement with the nuclear magnetic resonance scalar coupling constants from experiments. In turn, the model is used to gain insight into the physical driving forces behind secondary structure preferences of hydrated peptides. DFT calculations with and without the Conductor-like Screening Model (COSMO) show that the α helix is stabilized by solvent polarization due to dipole cooperativity. The two adjacent amide groups in β strand form a near-planar trapezoid that is not much larger than the size of water molecules. When the finite size of a water molecule is considered, the stabilization from solvent polarization for such a trapezoid is frustrated. Water molecules cannot find orientations to properly stabilize all four polar regions close to each other with such an awkward arrangement. This leads to quite substantial reduction in polarization stabilization. Although the polyproline II (PP-II) conformation is very similar to the β strand, the small twist in the backbone angles allowed much improved polarization stabilization. The improved polarization, when combined with favorable intrapeptide interactions, leads to the PP-II to be lowest in free energy. Other factors, such as the entropic TΔS and the ϕ, ψ coupling terms, are also studied but are found to play only a minor role. The insight shown in this work helps to better understand the structure of globular and intrinsic disordered proteins and facilitate future force field development.

Obtaining QM/MM binding free energies in the SAMPL8 drugs of abuse challenge: indirect approaches

Accurately predicting free energy differences is essential in realizing the full potential of rational drug design. Unfortunately, high levels of accuracy often require computationally expensive QM/MM Hamiltonians. Fortuitously, the cost of employing QM/MM approaches in rigorous free energy simulation can be reduced through the use of the so-called “indirect” approach to QM/MM free energies, in which the need for QM/MM simulations is avoided via a QM/MM “correction” at the classical endpoints of interest. Herein, we focus on the computation of QM/MM binding free energies in the context of the SAMPL8 Drugs of Abuse host–guest challenge. Of the 5 QM/MM correction coupled with force-matching submissions, PM6-D3H4/MM ranked submission proved the best overall QM/MM entry, with an RMSE from experimental results of 2.43 kcal/mol (best in ranked submissions), a Pearson’s correlation of 0.78 (second-best in ranked submissions), and a Kendall \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document}τ correlation of 0.52 (best in ranked submissions).

A comparison of three DFT exchange-correlation functionals and two basis sets for the prediction of the conformation distribution of hydrated polyglycine

The performance of three density functional theory (DFT) exchange-correlation functionals, namely, Perdew-Burke-Ernzerhof (PBE), BP86, and B3LYP, in predicting conformational distributions of a hydrated glycine peptide is tested with two different basis sets in the framework of adaptive force matching (AFM). The conformational distributions yielded the free energy profiles of the DFT functional and basis set combinations. Unlike traditional validations of potential energy and structural parameters, our approach allows the free energy of DFT to be validated. When compared to experimental distributions, the def2-TZVP basis set provides better agreement than a slightly trimmed aug-cc-pVDZ basis set. B3LYP is shown to be better than BP86 and PBE. The glycine model fitted against B3LYP-D3(BJ) with the def2-TZVP basis set is the most accurate and named the AFM2021 model for glycine. The AFM2021 glycine model provides better agreement with experimental J-coupling constants than C36m and ff14SB, although the margin is very small when compared to C36m. Our previously published alanine model is also refitted with the slightly simplified AFM2021 energy expression. This work shows good promise of AFM for developing force fields for a range of proteinogenic peptides using only DFT as reference.

A replica exchange umbrella sampling (REUS) approach to predict host-guest binding free energies in SAMPL8 challenge

In this study, we report binding free energy calculations of various drugs-of-abuse to Cucurbit-[8]-uril as part of the SAMPL8 blind challenge. Force-field parameters were obtained from force-matching with different quantum mechanical levels of theory. The Replica Exchange Umbrella Sampling (REUS) approach was used with a cylindrical restraint to enhance the sampling of host–guest binding. Binding free energy was calculated by pulling the guest molecule from one side of the symmetric and cylindrical host, then into and through the host, and out the other side (bidirectional) as compared to pulling only to the bound pose inside the cylindrical host (unidirectional). The initial results with force-matched MP2 parameter set led to RMSE of 4.68 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text{kcal}}/{\text{mol}}$$\end{document}kcal/mol from experimental values. However, the follow-up study with CHARMM generalized force field parameters and force-matched PM6-D3H4 parameters resulted in RMSEs from experiment of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.65$$\end{document}2.65 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.72 {\text{kcal}}/{\text{mol}}$$\end{document}1.72kcal/mol, respectively, which demonstrates the potential of REUS for accurate binding free energy calculation given a more suitable description of energetics. Moreover, we compared the free energies for the so called bidirectional and unidirectional free energy approach and found that the binding free energies were highly similar. However, one issue in the bidirectional approach is the asymmetry of profile on the two sides of the host. This is mainly due to the insufficient sampling for these larger systems and can be avoided by longer sampling simulations. Overall REUS shows great promise for binding free energy calculations.

Determining the hydration free energies of selected small molecules with MP2 and local MP2 through adaptive force matching

Force fields for seven small solute molecules, ethanol, 2-methyl-1-propanol, 2-butanol, cyclohexene, tetrahydropyran, 1,4-dioxane, and 1,4-butanediol, in dilute aqueous solutions were created with the adaptive force matching (AFM) method using MP2 or local MP2 as reference. The force fields provide a way to predict the hydration free energies (HFEs) of these molecules with only electronic structure calculations as reference. For six of the seven molecules, the predicted HFEs are in very good agreement with experiments. For 1,4-butanediol, the model created by force matching LMP2 provides a HFE that is too positive. Further investigation suggests that LMP2 may not be sufficiently accurate for computing HFEs for alcohols with AFM. Other properties, such as enthalpy of hydration, diffusion constants, and vibrational spectra, are also computed with the force field developed. The force fields developed by AFM provide a bridge for computing ensemble properties of the reference electronic structure method. With MP2 and LMP2 as reference methods, the computed properties of the small molecular solutes are found to be in good agreement with experiments.

Development and Validation of a DFT-Based Force Field for a Hydrated Homoalanine Polypeptide

A new force field has been created for simulating hydrated alanine polypeptides using the adaptive force matching (AFM) method. Only density functional theory calculations using the Perdew-Burke-Ernzerhof exchange-correlation functional and the D3 dispersion correction were used to fit the force field. The new force field, AFM2020, predicts NMR scalar coupling constants for hydrated homopolymeric alanine in better agreements with experimental data than several other models including those fitted directly to such data. For Ala₇, the new force field shows about 15% helical conformations, 20% conformation in the β basin, and 65% polyproline II. The predicted helical population of short hydrated alanine is higher than previous estimates based on the same experimental data. Gas-phase simulations indicate that the force field developed by AFM solution-phase data is likely to produce a reasonable conformation distribution when hydration water is no longer present, such as the interior of a protein.

Accurate MP2-based force fields predict hydration free energies for simple alkanes and alcohols in good agreement with experiments

Force fields for four small molecules, methane, ethane, methanol, and ethanol, were created by force matching MP2 gradients computed with triple-zeta-quality basis sets using the Adaptive Force Matching method. Without fitting to any experimental properties, the force fields created were able to predict hydration free energies, enthalpies of hydration, and diffusion constants in excellent agreements with experiments. The root mean square error for the predicted hydration free energies is within 1 kJ/mol of experimental measurements of Ben-Naim et al. [J. Chem. Phys. 81(4), 2016-2027 (1984)]. The good prediction of hydration free energies is particularly noteworthy, as it is an important fundamental property. Similar hydration free energies of ethane relative to methane and of ethanol relative to methanol are attributed to a near cancellation of cavitation penalty and favorable contributions from dispersion and Coulombic interactions as a result of the additional methyl group.

Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery

Predicting protein-ligand binding affinities and the associated thermodynamics of biomolecular recognition is a primary objective of structure-based drug design. Alchemical free energy simulations offer a highly accurate and computationally efficient route to achieving this goal. While the AMBER molecular dynamics package has successfully been used for alchemical free energy simulations in academic research groups for decades, widespread impact in industrial drug discovery settings has been minimal because of the previous limitations within the AMBER alchemical code, coupled with challenges in system setup and postprocessing workflows. Through a close academia-industry collaboration we have addressed many of the previous limitations with an aim to improve accuracy, efficiency, and robustness of alchemical binding free energy simulations in industrial drug discovery applications. Here, we highlight some of the recent advances in AMBER20 with a focus on alchemical binding free energy (BFE) calculations, which are less computationally intensive than alternative binding free energy methods where full binding/unbinding paths are explored. In addition to scientific and technical advances in AMBER20, we also describe the essential practical aspects associated with running relative alchemical BFE calculations, along with recommendations for best practices, highlighting the importance not only of the alchemical simulation code but also the auxiliary functionalities and expertise required to obtain accurate and reliable results. This work is intended to provide a contemporary overview of the scientific, technical, and practical issues associated with running relative BFE simulations in AMBER20, with a focus on real-world drug discovery applications.

Altering the Solubility of the Antibiotic Candidate Nisin-A Computational Study

The growing bacterial resistance to available antibiotics makes it necessary to look for new drug candidates. An example is the lanthionine-containing nisin, which has a broad spectrum of antimicrobial activity. While nisin is widely utilized as a food preservative, its poor solubility and low stability at physiological pH hinder its use as an antibiotic. As the solubility of nisin is controlled by the residues of the hinge region, we have performed molecular dynamics simulations of various mutants and studied their effects on nisin's solubility. These simulations are complicated by the presence of two uncommon residues (dehydroalanine and dehydrobutyrine) in the peptide. The primary goal of the present study is to derive rules for designing new mutants that will be more soluble at physiological pH and, therefore, may serve as a basis for the future antibiotic design. Another aim of our study is to evaluate whether existing force fields can model the solubility of these amino acids accurately in order to motivate further developments of force fields to account for solubility information.

Exploring optimization strategies for improving explicit water models: Rigid n-point model and polarizable model based on Drude oscillator

Rigid n-point water models are widely used in atomistic simulations, but have known accuracy drawbacks. Increasing the number of point charges, as well as adding electronic polarizability, are two common strategies for accuracy improvements. Both strategies come at considerable computational cost, which weighs heavily against modest possible accuracy improvements in practical simulations. In an effort to provide guidance for model development, here we have explored the limiting accuracy of "electrostatically globally optimal" n-point water models in terms of their ability to reproduce properties of water dimer-a mimic of the condensed state of water. For a given n, each model is built upon a set of reference multipole moments (e.g. ab initio) and then optimized to reproduce water dimer total dipole moment. The models are then evaluated with respect to the accuracy of reproducing the geometry of the water dimer. We find that global optimization of the charge distribution alone can deliver high accuracy of the water model: for n = 4 or n = 5, the geometry of the resulting water dimer can be almost within 50 of the ab initio reference, which is half that of the experimental error margin. Thus, global optimization of the charge distribution of classical n-point water models can lead to high accuracy models. We also find that while the accuracy improvement in going from n = 3 to n = 4 is substantial, the additional accuracy increase in going from n = 4 to n = 5 is marginal. Next, we have explored accuracy limitations of the standard practice of adding electronic polarizability (via a Drude particle) to a "rigid base"-pre-optimization rigid n-point water model. The resulting model (n = 3) shows a relatively small improvement in accuracy, suggesting that the strategy of merely adding the polarizability to an inferior accuracy water model used as the base cannot fix the defects of the latter. An alternative strategy in which the parameters of the rigid base model are globally optimized along with the polarizability parameter is much more promising: the resulting 3-point polarizable model out-performs even the 5-point optimal rigid model by a large margin. We suggest that future development efforts consider 3- and 4-point polarizable models where global optimization of the "rigid base" is coupled to optimization of the polarizability to deliver globally optimal solutions.

Development of a Robust Indirect Approach for MM → QM Free Energy Calculations That Combines Force-Matched Reference Potential and Bennett's Acceptance Ratio Methods

We use the PBE0/6-31G* density functional method to perform ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations under periodic boundary conditions with rigorous electrostatics using the ambient potential composite Ewald method in order to test the convergence of MM → QM/MM free energy corrections for the prediction of 17 small-molecule solvation free energies and eight ligand binding free energies to T4 lysozyme. The "indirect" thermodynamic cycle for calculating free energies is used to explore whether a series of reference potentials improve the statistical quality of the predictions. Specifically, we construct a series of reference potentials that optimize a molecular mechanical (MM) force field's parameters to reproduce the ab initio QM/MM forces from a QM/MM simulation. The optimizations form a systematic progression of successively expanded parameters that include bond, angle, dihedral, and charge parameters. For each reference potential, we calculate benchmark quality reference values for the MM → QM/MM correction by performing the mixed MM and QM/MM Hamiltonians at 11 intermediate states, each for 200 ps. We then compare forward and reverse application of Zwanzig's relation, thermodynamic integration (TI), and Bennett's acceptance ratio (BAR) methods as a function of reference potential, simulation time, and the number of simulated intermediate states. We find that Zwanzig's equation is inadequate unless a large number of intermediate states are explicitly simulated. The TI and BAR mean signed errors are very small even when only the end-state simulations are considered, and the standard deviations of the TI and BAR errors are decreased by choosing a reference potential that optimizes the bond and angle parameters. We find a robust approach for the data sets of fairly rigid molecules considered here is to use bond + angle reference potential together with the end-state-only BAR analysis. This requires QM/MM simulations to be performed in order to generate reference data to parametrize the bond + angle reference potential, and then this same simulation serves a dual purpose as the full QM/MM end state. The convergence of the results with respect to time suggests that computational resources may be used more efficiently by running multiple simulations for no more than 50 ps, rather than running one long simulation.

Surface Penetration without Enrichment: Simulations Show Ion Surface Propensities Consistent with Both Elevated Surface Tension and Surface Sensitive Spectroscopy

Molecular dynamics was used to investigate ion surface propensities in NaCl, KBr, and CsI solutions with an MP2-based force field. Although NaCl is found to be strongly repelled from the liquid-vapor interface, softer ions, such as I-, penetrate closely to the interface. Despite the surface penetration, the concentration of CsI near the interface is still lower than that in the bulk, thus leading to no surface enrichment. The salt concentration is found to affect relative surface propensities of the ions. More significant surface penetration is observed at higher salt concentrations. Softer ions at higher concentrations form a complex multilayer arrangement that can not be characterized as a simple surface bilayer. The simulated ion distributions explain the spectroscopic evidence of surface perturbation by soft ions with a negative surface excess consistent with an increased surface tension of salt solutions.

The Good, the Bad, and the Ugly: "HiPen", a New Dataset for Validating (S)QM/MM Free Energy Simulations

Indirect (S)QM/MM free energy simulations (FES) are vital to efficiently incorporating sufficient sampling and accurate (QM) energetic evaluations when estimating free energies of practical/experimental interest. Connecting between levels of theory, i.e., calculating Δ A l o w → h i g h , remains to be the most challenging step within an indirect FES protocol. To improve calculations of Δ A l o w → h i g h , we must: (1) compare the performance of all FES methods currently available; and (2) compile and maintain datasets of Δ A l o w → h i g h calculated for a wide-variety of molecules so that future practitioners may replicate or improve upon the current state-of-the-art. Towards these two aims, we introduce a new dataset, "HiPen", which tabulates Δ A g a s M M → 3 o b (the free energy associated with switching from an M M to an S C C - D F T B molecular description using the 3ob parameter set in gas phase), calculated for 22 drug-like small molecules. We compare the calculation of this value using free energy perturbation, Bennett's acceptance ratio, Jarzynski's equation, and Crooks' equation. We also predict the reliability of each calculated Δ A g a s M M → 3 o b by evaluating several convergence criteria including sample size hysteresis, overlap statistics, and bias metric ( Π ). Within the total dataset, three distinct categories of molecules emerge: the "good" molecules, for which we can obtain converged Δ A g a s M M → 3 o b using Jarzynski's equation; "bad" molecules which require Crooks' equation to obtain a converged Δ A g a s M M → 3 o b ; and "ugly" molecules for which we cannot obtain reliably converged Δ A g a s M M → 3 o b with either Jarzynski's or Crooks' equations. We discuss, in depth, results from several example molecules in each of these categories and describe how dihedral discrepancies between levels of theory cause convergence failures even for these gas phase free energy simulations.

Accelerating QM/MM Free Energy Computations via Intramolecular Force Matching

The calculation of free energy differences between levels of theory has numerous potential pitfalls. Chief among them is the lack of overlap, i.e., ensembles generated at one level of theory (e.g., "low") not being good approximations of ensembles at the other (e.g., "high"). Numerous strategies have been devised to mitigate this issue. However, the most straightforward approach is to ensure that the "low" level ensemble more closely resembles that of the "high". Ideally, this is done without increasing computational cost. Herein, we demonstrate that by reparametrizing classical intramolecular potentials to reproduce high level forces (i.e., force matching) configurational overlap between a "low" (i.e., classical) and "high" (i.e., quantum) level can be significantly improved. This procedure is validated on two test cases and results in vastly improved convergence of free energy simulations.

Force matching as a stepping stone to QM/MM CB[8] host/guest binding free energies: a SAMPL6 cautionary tale

Use of quantum mechanical/molecular mechanical (QM/MM) methods in binding free energy calculations, particularly in the SAMPL challenge, often fail to achieve improvement over standard additive (MM) force fields. Frequently, the implementation is through use of reference potentials, or the so-called "indirect approach", and inherently relies on sufficient overlap existing between MM and QM/MM configurational spaces. This overlap is generally poor, particularly for the use of free energy perturbation to perform the MM to QM/MM free energy correction at the end states of interest (e.g., bound and unbound states). However, by utilizing MM parameters that best reproduce forces obtained at the desired QM level of theory, it is possible to lessen the configurational disparity between MM and QM/MM. To this end, we sought to use force matching to generate MM parameters for the SAMPL6 CB[8] host-guest binding challenge, classically compute binding free energies, and apply energetic end state corrections to obtain QM/MM binding free energy differences. For the standard set of 11 molecules and the bonus set (including three additional challenge molecules), error statistics, such as the root mean square deviation (RMSE) were moderately poor (5.5 and 5.4 kcal/mol). Correlation statistics, however, were in the top two for both standard and bonus set submissions ([Formula: see text] of 0.42 and 0.26, [Formula: see text] of 0.64 and 0.47 respectively). High RMSE and moderate correlation strongly indicated the presence of systematic error. Identifiable issues were ameliorated for two of the guest molecules, resulting in a reduction of error and pointing to strong prospects for the future use of this methodology.

Development of reactive force fields using ab initio molecular dynamics simulation minimally biased to experimental data

Incorporation of quantum mechanical electronic structure data is necessary to properly capture the physics of many chemical processes. Proton hopping in water, which involves rearrangement of chemical and hydrogen bonds, is one such example of an inherently quantum mechanical process. Standard ab initio molecular dynamics (AIMD) methods, however, do not yet accurately predict the structure of water and are therefore less than optimal for developing force fields. We have instead utilized a recently developed method which minimally biases AIMD simulations to match limited experimental data to develop novel multiscale reactive molecular dynamics (MS-RMD) force fields by using relative entropy minimization. In this paper, we present two new MS-RMD models using such a parameterization: one which employs water with harmonic internal vibrations and another which uses anharmonic water. We show that the newly developed MS-RMD models very closely reproduce the solvation structure of the hydrated excess proton in the target AIMD data. We also find that the use of anharmonic water increases proton hopping, thereby increasing the proton diffusion constant.

Water graphene contact surface investigated by pairwise potentials from force-matching PAW-PBE with dispersion correction

A pairwise additive atomistic potential was developed for modeling liquid water on graphene. The graphene-water interaction terms were fit to map the PAW-PBE-D3 potential energy surface using the adaptive force matching method. Through condensed phase force matching, the potential developed implicitly considers the many-body effects of water. With this potential, the graphene-water contact angle was determined to be 86° in good agreement with a recent experimental measurement of 85° ± 5° on fully suspended graphene. Furthermore, the PAW-PBE-D3 based model was used to study contact line hysteresis. It was found that the advancing and receding contact angles of water do agree on pristine graphene, however a long simulation time was required to reach the equilibrium contact angle. For water on suspended graphene, sharp peaks in the water density profile disappear when the flexibility of graphene was explicitly considered. The water droplet induces graphene to wrap around it leading to a slightly concave contact interface.

The strengths and limitations of effective centroid force models explored by studying isotopic effects in liquid water

The development of effective centroid potentials (ECPs) is explored with both the constrained-centroid and quasi-adiabatic force matching using liquid water as a test system. A trajectory integrated with the ECP is free of statistical noises that would be introduced when the centroid potential is approximated on the fly with a finite number of beads. With the reduced cost of ECP, challenging experimental properties can be studied in the spirit of centroid molecular dynamics. The experimental number density of H2O is 0.38% higher than that of D2O. With the ECP, the H2O number density is predicted to be 0.42% higher, when the dispersion term is not refit. After correction of finite size effects, the diffusion constant of H2O is found to be 21% higher than that of D2O, which is in good agreement with the 29.9% higher diffusivity for H2O observed experimentally. Although the ECP is also able to capture the redshifts of both the OH and OD stretching modes in liquid water, there are a number of properties that a classical simulation with the ECP will not be able to recover. For example, the heat capacities of H2O and D2O are predicted to be almost identical and higher than the experimental values. Such a failure is simply a result of not properly treating quantized vibrational energy levels when the trajectory is propagated with classical mechanics. Several limitations of the ECP based approach without bead population reconstruction are discussed.

Fibril-Barrel Transitions in Cylindrin Amyloids

We introduce Replica-Exchange-with-Tunneling (RET) simulations as a tool for studies of the conversion between polymorphic amyloids. For the 11-residue amyloid-forming cylindrin peptide we show that this technique allows for a more efficient sampling of the formation and interconversion between fibril-like and barrel-like assemblies. We describe a protocol for optimized analysis of RET simulations that allows us to propose a mechanism for formation and interconversion between various cylindrin assemblies. Especially, we show that an interchain salt bridge between residues K3 and D7 is crucial for formation of the barrel structure.

Accurate Prediction of the Hydration Free Energies of 20 Salts through Adaptive Force Matching and the Proper Comparison with Experimental References

Simple pairwise potentials for five alkali ions and four halide ions were developed by only fitting to ab initio MP2 forces with the adaptive force matching (AFM) method. Without fitting to any experimental information, the AFM models predict the hydration free energies of all 10 fluoride and chloride salts formed by these ions within 1.5% of experimental references. The predicted hydration free energies for the 10 bromide and iodide salts are within 5-6% of experimental references with the larger error likely due to the neglect of explicit treatment of polarization and charge transfer. An inconsistency in the treatment of the gas phase entropy term between experimental and theoretical approaches is discussed. A new simplified hydration free energy for the ions is reported for use as a more appropriate experimental reference for further theoretical studies. The simulations show different dipole alignments for the hydration waters of cations and anions. While hydration waters of small cations tend to align their molecular dipole toward the ion, the dipole of one of the water OH bonds is aligned with the field of an anion.

Solvation Dynamics of CO₂(g) by Monoethanolamine at the Gas-Liquid Interface: A Molecular Mechanics Approach

A classical force field approach was used to characterize the solvation dynamics of high-density CO₂(g) by monoethanolamine (MEA) at the air–liquid interface. Intra- and intermolecular CO₂ and MEA potentials were parameterized according to the energetics calculated at the MP2 and BLYP-D2 levels of theory. The thermodynamic properties of CO₂ and MEA, such as heat capacity and melting point, were consistently predicted using this classical potential. An approximate interfacial simulation for CO₂(g)/MEA(l) was performed to monitor the depletion of the CO₂(g) phase, which was influenced by amino and hydroxyl groups of MEA. There are more intramolecular hydrogen bond interactions notably identified in the interfacial simulation than the case of bulk MEA(l) simulation. The hydroxyl group of MEA was found to more actively approach CO₂ and overpower the amino group to interact with CO₂ at the air–liquid interface. With artificially reducing the dipole moment of the hydroxyl group, CO₂–amino group interaction was enhanced and suppressed CO₂(g) depletion. The hydroxyl group of MEA was concluded to play dual but contradictory roles for CO₂ capture.

Possible Evidence for a New Form of Liquid Buried in the Surface Tension of Supercooled Water

Contrary to the historical data, several recent experiments indicate that the surface tension of supercooled water follows a smooth extrapolation of the IAPWS equation in the supercooled regime. It can be seen, however, that a small deviation from the IAPWS equation is present in the recent experimental measurements. It is shown with simulations using the WAIL water potential that the small deviation in the experimental data is consistent with the tail of an exponential growth in surface tension as temperature decreases. The emergence temperature, Te, of a substantial deviation from the IAPWS equation is shown to be 227 K for the WAIL water and 235 K for real water. Since the 227 K Te is close to the Widom line in WAIL water, we argue that real water at 235 K approaches a similar crossover line at one atmospheric pressure.

The effects of replacing the water model while decoupling water-water and water-solute interactions on computed properties of simple salts

The effects of decoupling the water-water and water-solute interactions are studied with selected mono-valent ions as the solute. Using the ion-water cross terms developed for the BLYPSP-4F water model, we replaced the water potential with WAIL, TIP4P, and TIP3P without changing the ion-water parameters. When the adaptive force matching (AFM) derived BLYPSP-4F model is replaced by the other AFM derived WAIL model, the difference in ion properties, such as hydration free energies, radial distribution functions, relative diffusion constants, is negligible, demonstrating the feasibility for combining AFM parameters from different sources. Interestingly, when the AFM-derived ion-water cross-terms are used with a non-AFM based water model, only small changes in the ion properties are observed. The final combined models with TIP3P or TIP4P water reproduce the salt hydration free energies within 6% of experiments. The feasibility of combining AFM models with other non-AFM models is of significance since such combinations allow more complex systems to be studied without specific parameterization. In addition, the study suggests an interesting prospect of reusing the cross-terms when a part of a general force field is replaced with a different model. The prevailing practice, which is to re-derive all cross-terms with combining rules, may not have been optimal.

Computational investigation of structure, dynamics and nucleation kinetics of a family of modified Stillinger-Weber model fluids in bulk and free-standing thin films

In recent years, computer simulations have found increasingly widespread use as powerful tools for studying phase transitions in wide variety of systems. In the particular and very important case of aqueous systems, the commonly used force-fields tend to offer quite different predictions with respect to a wide range of thermodynamic and kinetic properties, including the ease of ice nucleation, the propensity to freeze at a vapor-liquid interface, and the existence of a liquid-liquid phase transition. It is thus of fundamental and practical interest to understand how different features of a given water model affect its thermodynamic and kinetic properties. In this work, we use the forward-flux sampling technique to study the crystallization kinetics of a family of modified Stillinger-Weber (SW) potentials with energy (ε) and length (σ) scales taken from the monoatomic water (mW) model, but with different tetrahedrality parameters (λ). By increasing λ from 21 to 24, we observe the nucleation rate increases by 48 orders of magnitude at a supercooling of ζ = T/Tm = 0.845. Using classical nucleation theory, we are able to demonstrate that this change can largely be accounted for by the increase in |Δμ|, the thermodynamic driving force. We also perform rate calculations in freestanding thin films of the supercooled liquid, and observe a crossover from surface-enhanced crystallization at λ = 21 to bulk-dominated crystallization for λ ≥ 22.

Pairwise-additive force fields for selected aqueous monovalent ions from adaptive force matching

Simple non-polarizable potentials were developed for Na(+), K(+), Cl(-), and Br(-) using the adaptive force matching (AFM) method with ab initio MP2 method as reference. Our MP2-AFM force field predicts the solvation free energies of the four salts formed by the ions with an error of no more than 5%. Other properties such as the ion-water radial distribution functions, first solvation shell water tilt angle distributions, ion diffusion constants, concentration dependent diffusion constant of water, and concentration dependent surface tension of the solutions were calculated with this potential. Very good agreement was achieved for these properties. In particular, the diffusion constants of the ions are within 6% of experimental measurements. The model predicts bromide to be enriched at the interface in the 1.6M KBr solution but predicts the ion to be repelled for the surface at lower concentration.

The liquid-vapor equilibria of TIP4P/2005 and BLYPSP-4F water models determined through direct simulations of the liquid-vapor interface

In this paper, the surface tension and critical properties for the TIP4P/2005 and BLYPSP-4F models are reported. A clear dependence of surface tension on the van der Waals cutoff radius (rvdw) is shown when van der Waals interactions are modeled with a simple cutoff scheme. A linear extrapolation formula is proposed that can be used to determine the infinite rvdw surface tension through a few simulations with finite rvdw. A procedure for determining liquid and vapor densities is proposed that does not require fitting to a profile function. Although the critical temperature of water is also found to depend on the choice of rvdw, the dependence is weaker. We argue that a rvdw of 1.75 nm is a good compromise for water simulations when long-range van der Waals correction is not applied. Since the majority of computational programs do not support rigorous treatment of long-range dispersion, the establishment of a minimal acceptable rvdw is important for the simulation of a variety of inhomogeneous systems, such as water bubbles, and water in confined environments. The BLYPSP-4F model predicts room temperature surface tension marginally better than TIP4P/2005 but overestimates the critical temperature. This is expected since only liquid configurations were fit during the development of the BLYPSP-4F potential. The potential is expected to underestimate the stability of vapor and thus overestimate the region of stability for the liquid.

Continuous and Discontinuous Dynamic Crossover in Supercooled Water in Computer Simulations

The dynamic crossover behavior of supercooled water as described by the first-principle based WAIL potential was investigated. Below the second liquid-liquid critical point, the viscosity shows a discontinuous jump consistent with a first-order phase transition between the high density liquid and the low density liquid. Above the critical point, a continuous transition occurs with only the first derivative of viscosity being discontinuous, and the dynamic crossover temperature is about 8 K below the thermodynamic switchover temperature. The 8 K shift can be explained by a delay in dynamic crossover, which does not occur until the more viscous liquid starts to dominate the population and jams the flow. On the basis of finite-size effects observed in our simulations, we believe that dynamic discontinuity may be observable above the critical point in confined water when the confinement is on a length scale shorter than the spatial correlation.

The adaptive buffered force QM/MM method in the CP2K and AMBER software packages

The implementation and validation of the adaptive buffered force (AdBF) quantum-mechanics/molecular-mechanics (QM/MM) method in two popular packages, CP2K and AMBER are presented. The implementations build on the existing QM/MM functionality in each code, extending it to allow for redefinition of the QM and MM regions during the simulation and reducing QM-MM interface errors by discarding forces near the boundary according to the buffered force-mixing approach. New adaptive thermostats, needed by force-mixing methods, are also implemented. Different variants of the method are benchmarked by simulating the structure of bulk water, water autoprotolysis in the presence of zinc and dimethyl-phosphate hydrolysis using various semiempirical Hamiltonians and density functional theory as the QM model. It is shown that with suitable parameters, based on force convergence tests, the AdBF QM/MM scheme can provide an accurate approximation of the structure in the dynamical QM region matching the corresponding fully QM simulations, as well as reproducing the correct energetics in all cases. Adaptive unbuffered force-mixing and adaptive conventional QM/MM methods also provide reasonable results for some systems, but are more likely to suffer from instabilities and inaccuracies. © 2015 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Building Water Models: A Different Approach

Simplified classical water models are currently an indispensable component in practical atomistic simulations. Yet, despite several decades of intense research, these models are still far from perfect. Presented here is an alternative approach to constructing widely used point charge water models. In contrast to the conventional approach, we do not impose any geometry constraints on the model other than the symmetry. Instead, we optimize the distribution of point charges to best describe the "electrostatics" of the water molecule. The resulting "optimal" 3-charge, 4-point rigid water model (OPC) reproduces a comprehensive set of bulk properties significantly more accurately than commonly used rigid models: average error relative to experiment is 0.76%. Close agreement with experiment holds over a wide range of temperatures. The improvements in the proposed model extend beyond bulk properties: compared to common rigid models, predicted hydration free energies of small molecules using OPC are uniformly closer to experiment, with root-mean-square error <1 kcal/mol.

Graphene: a partially ordered non-periodic solid

Molecular dynamics simulations were performed to study the structural features of graphene over a wide range of temperatures from 50 to 4000 K using the PPBE-G potential [D. Wei, Y. Song, and F. Wang, J. Chem. Phys. 134, 184704 (2011)]. This potential was developed by force matching the Perdew-Burke-Ernzerhof (PBE) exchange correlation functional and has been validated previously to provide accurate potential energy surface for graphene at temperatures as high as 3000 K. Simulations with the PPBE‑G potential are the best available approximation to a direct Car-Parrinello Molecular Dynamics study of graphene. One advantage of the PBE-G potential is to allow large simulation boxes to be modeled efficiently so that properties showing strong finite size effects can be studied. Our simulation box contains more than 600,000 C atoms and is one of the largest graphene boxes ever modeled. With the PPBE-G potential, the thermal-expansion coefficient is negative up to 4000 K. With a large box and an accurate potential, the critical exponent for the scaling properties associated with the normal-normal and height-height correlation functions was confirmed to be 0.85. This exponent remains constant up to 4000 K suggesting graphene to be in the deeply cooled regime even close to the experimental melting temperature. The reduced peak heights in the radial distribution function of graphene show an inverse power law dependence to distance, which indicates that a macroscopic graphene sheet will lose long-range crystalline order as predicted by the Mermin-Wagner instability. Although graphene loses long-range translational order, it retains long range orientational order as indicated by its orientational correlation function; graphene is thus partially ordered but not periodic.

Systematic improvement of a classical molecular model of water

We report the iAMOEBA ("inexpensive AMOEBA") classical polarizable water model. The iAMOEBA model uses a direct approximation to describe electronic polarizability, in which the induced dipoles are determined directly from the permanent multipole electric fields and do not interact with one another. The direct approximation reduces the computational cost relative to a fully self-consistent polarizable model such as AMOEBA. The model is parameterized using ForceBalance, a systematic optimization method that simultaneously utilizes training data from experimental measurements and high-level ab initio calculations. We show that iAMOEBA is a highly accurate model for water in the solid, liquid, and gas phases, with the ability to fully capture the effects of electronic polarization and predict a comprehensive set of water properties beyond the training data set including the phase diagram. The increased accuracy of iAMOEBA over the fully polarizable AMOEBA model demonstrates ForceBalance as a method that allows the researcher to systematically improve empirical models by efficiently utilizing the available data.

Liquid-liquid transition in supercooled water suggested by microsecond simulations

The putative liquid-liquid phase transition in supercooled water has been used to explain many anomalous behaviors of water. However, no direct experimental verification of such a phase transition has been accomplished, and theoretical studies from different simulations contradict each other. We investigated the putative liquid-liquid phase transition using the Water potential from Adaptive Force Matching for Ice and Liquid (WAIL). The simulation reveals a first-order phase transition in the supercooled regime with the critical point at ~207 K and 50 MPa. Normal water is high-density liquid (HDL). Low-density liquid (LDL) emerges at lower temperatures. The LDL phase has a density only slightly larger than that of the ice-Ih and shows more long-range order than HDL. However, the transformation from LDL to HDL is spontaneous across the first-order phase transition line, suggesting the LDL configuration is not poorly formed nanocrystalline ice. It has been demonstrated in the past that the WAIL potential provides reliable predictions of water properties such as melting temperature and temperature of maximum density. Compared with other simple water potentials, WAIL is not biased by fitting to experimental properties, and simulation with this potential reflects the prediction of a high-quality first-principle potential energy surface.

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.

Static dielectric constants and molecular dipole distributions of liquid water and ice-Ih investigated by the PAW-PBE exchange-correlation functional

The static dielectric constants, ε(s), of ice-Ih and liquid water were calculated using density functional theory with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional and the projector-augmented-wave (PAW) approach. Proton disordered ice configurations and uncorrelated liquid configurations were sampled with the electrostatic switching method using force fields specially designed to facilitate the ab initio free energy perturbation calculations. Our results indicate that PAW-PBE underestimates the ε(s) of both ice-Ih and liquid water but predicts the ratio of ice and water ε(s) in good agreement with experimental measurements. PAW-PBE gives average water dipole moments of 2.50 D in ice-Ih and 2.48 D in the liquid. Our results show that the fixed-charge water models developed by adaptive force matching can reproduce the PAW-PBE dipole moments with an error of approximately 5%. The ice and liquid models created in this work have polarizabilities of 1.32 Å(3) and 1.30 Å(3), respectively, along the HOH bisector direction.

Predicting the melting temperature of ice-Ih with only electronic structure information as input

The melting temperature of ice-Ih was calculated with only electronic structure information as input by creating a problem-specific force field. The force field, Water model by AFM for Ice and Liquid (WAIL), was developed with the adaptive force matching (AFM) method by fitting to post-Hartree-Fock quality forces obtained in quantum mechanics∕molecular mechanics calculations. WAIL predicts the ice-Ih melting temperature to be 270 K. The model also predicts the densities of ice and water, the temperature of maximum density of water, the heat of vaporizations, and the radial distribution functions for both ice and water in good agreement with experimental measurements. The non-dissociative WAIL model is very similar to a flexible version of the popular TIP4P potential and has comparable computational cost. By customizing to problem-specific configurations with the AFM approach, the resulting model is remarkably more accurate than any variants of TIP4P for simulating ice-Ih and water in the temperature range from 253 K and 293 K under ambient pressure.