Simulations of water and hydrophobic hydration using a neural network potential
SSimulations of water and hydrophobic hydration using a neural network potential
Alexander S. Lyons, Jr. and Steven W. Rick ∗ Department of Chemistry, University of New Orleans, New Orleans, LA, 70148
Using a neural network potential (ANI-1ccx) generated from quantum data on a large data setof molecules and pairs of molecules, isothermal, constant volume simulations demonstrate that themodel can be as accurate as ab initio molecular dynamics for simulations of pure liquid water andthe aqueous solvation of a methane molecule. No theoretical or experimental data for the liquidphase is used to train the model, suggesting that the ANI-1ccx approach is an effective method tolink high level ab initio methods to potentials for large scale simulations.
I. INTRODUCTION
Conventional molecular simulations require a model, orforce field, which can involve developing potential func-tions that capture the essential physics of the systemand optimizing a large set of parameters. The param-eterization often uses a combination of quantum calcula-tions and experimental data. The resulting models canin many cases be accurate and efficient, but the time re-quired in their development and the need for experimen-tal data can be limiting, especially for the simulation ofnew materials. Simulations without models can be car-ried out using ab initio molecular dynamics (AIMD),[1, 2]however these are limited in system size and time scalesand can require using relatively low level density func-tional theory (DFT) methods.Machine learning (ML) methods offer a means to carryout simulations with high accuracy without a potentialfunction or experimental data.[3–10] These methods use ab initio calculations to train a neural network represen-tation of a potential energy surface. The goal of thesemethods is to develop potentials that can be used just as ab initio methods, so they can treat new materials and re-acting systems, but are much faster than AIMD. Centralto these methods are descriptors, or generalized coordi-nates, that transfer the high dimensional set of coordi-nates of the system to a set of inputs, the descriptors, forthe ML algorithms. Several descriptors have been devel-oped that satisfy the necessary physical conditions of be-ing invariant with respect to translation and rotation, aswell as exchange of identical atoms.[3, 5–7, 11–14] Otherapproaches, which do not build in these symmetries, havebeen successful as well.[15] For molecular dynamics sim-ulations, the descriptors should also be continuous anddifferentiable so they can be used to calculate forces.[16]Using these descriptors, artificial neural networks (NN)with multiple layers learn to reproduce ab initio energiesfor a selected[12, 16] reference set of molecules. Whilethe resulting models are one to two orders of magnitudeslower than conventional force fields, they can be sev-eral orders of magnitude faster to compute than DFTmethods.[16] ∗ Electronic address: [email protected]
One of the most promising features of ML in generalis their applicability to systems for which they have notbeen trained. In this work, we will use a NN potential tosimulate liquid water and a methane molecule solvatedby water. Water and hydrophobic solvation were chosenfor their significance and the abundance of experimentaldata for comparison. NN potentials for water have beendeveloped and have shown high accuracy in reproduc-ing the properties of liquid water and ice over a range ofstate points.[17–22] These models were trained using wa-ter clusters or liquid water, specifically for water. Someuse a purely machine-learning approach[17, 18, 20, 21]while others add combine ML with physical potentials,such as Coulombic interactions.[19, 22] Other approachesuse ML to optimize parameters for a classical forcefield.[23, 24] These and other ML methods for molec-ular simulations,[4, 7, 11, 13, 14] have recently beencompared.[25]. The point of the present work is not todevelop or compare water models, but to test a generalpurpose NN potential. Can a model not trained for wa-ter, or the condensed phases, be effective for liquid water?The NN potential of Roitberg and co-workers.[9, 10, 12]was chosen for its generality and ease of use. This NNpotential, ANI-1ccx,[10] was trained for molecules madeup of H, C, N, and O atoms on a large, diverse data setto reproduce energies at the coupled cluster level. Thetraining dataset contains structures of single moleculesand dimers only.[9] We use the model for simulations inthe liquid phase, as a test of the transferability of themodel to the condensed phase.
II. METHODS
The ANI-1ccx potential.
Full details of the modelhave been published by Roitberg and co-workers;[9]here we aim to provide a summary of the main ideas.This method builds on the NN potential of Behler andParrinello.[5] The total energy of the system, E, is givenby a sum over a contribution from each atom, i , E = (cid:88) i E i (1)where E i is the atomic contribution to the energy. TheNN is trained to reproduce E i for each atom, with a dif-ferent NN for each element. Each element is treated by a r X i v : . [ phy s i c s . c h e m - ph ] J a n the same set of NN, independent of its chemical environ-ment. ANI-1ccx has one NN for each element with fivelayers, with 3 hidden layers between the input and outputlayers, each with a different number of nodes. The NNtakes as input descriptors, termed atom centered sym-metry functions[5] or atomic environment variables[12]which are determined from the three dimensional struc-ture. There are two types of symmetry functions, radialand angular. The radial functions are given by G Ri = (cid:88) j (cid:54) = i e − η ( r ij − R s ) f c ( r ij ) (2)where r ij is the distance between atoms i and j, η andR s are fixed parameters for a particular node, and f c isa cut-off function that smoothly goes to zero at a lengthR C . In the ANI-1ccx potential, the η parameter is fixedand multiple R S values are used, equally spaced out tothe cut-off length of 5.2 ˚A.The angular function is G Ai = 2 − ζ (cid:88) j,k (cid:54) = i (1 + cos( θ ijk − θ S )) ζ exp (cid:104) − η (( r ij + r ik ) / − R s ) (cid:105) f c ( r ij ) f c ( r ik ) (3)where θ ijk is the angle between atoms i , j , and k , with i at the vertex. The parameters θ S and R S center thefunction at a particular angle and distance, with ζ and η determining the widths of the functions. This form is dif-ferent from that of Behler and Parrinello and allows forgreater spacial and angular resolution of an atom’s en-vironment. The angular function establishes the angularpotentials between both bonded and non-bonded atoms,allowing for many-body polarization to be treated at thethree-body level.It is worth emphasizing that these potentials are notusing ML to optimize a potential model, although thatcan be done.[24]. Equations 2 and 3 are fairly agnosticfunctions, acting as windows to represent an atom’s envi-ronment. The two main approximations that go into themodels are that they do not go beyond a cut-off lengthof R C and non-additive effects are only included at thethree-body level.The ANI-1ccx model was trained using 5.2 million con-formations of 64 thousand molecules containing C,N,O,and H atoms with energies evaluated at the DFT level.Of those, 500 thousand were selected for retraining usingtransfer learning[26]. The retraining was done using ener-gies calculated approaching the accuracy of the coupledcluster with full treatment of singles and doubles, andtriples calculated with perturbation theory extrapolatedto the complete basis let limit (CCSD (T)/CBS),[10]using the domain localized, local pair natural orbitalDPLNO-CCSD(T) method.[27] The conformations in-cluded dimers, energy minimized structures, and struc-tures generated from normal mode sampling to generatestructures away from energy minima.[9, 12] Simulation details.
Simulations were done usingthe Atomic Simulation Environment (ASE) package.[28]The ANI-1ccx potential implemented in ASE wasdownloaded from GitHub (https: // github.com / isayev / ASE ANI). The liquid water simulations used 128molecules and the water/methane simulations used 128water plus one methane molecule. Simulations were runin the NVT ensemble at a density of 1.0 g cm − , us-ing Langevin thermostatting, with a time step of 0.1 fs,and at a temperature of 300 K. Initial structures weretaken from equilibrated simulations using the SPC/Emodel.[29] The simulations were equilibrated for 5 ps anddata was collected over an additional 10 ps for water.The methane/water simulations collected data over 10ps. The simulation times are similar to those used inAIMD simulations of water, which can range from 5 to50 ps.[30–35] In our hands, with a small 4-node GPU-enabled (NVIDIA K-80; Intel Haswell) compute cluster,1 ps of simulation time takes about 5 days. The relativelylong wall clock times are not due to inherent computa-tional cost of the ANI-1ccx potential, but in its imple-mentation. The ASE program, written in python, inte-grates well with the ANI ML potential, but it is not asfast as other simulation packages.The dynamical properties, the diffusion constant andthe Nuclear Magnetic Resonance (NMR) relaxation time, τ NMR , were calculated using 6 NVE simulations of 2 pseach. The diffusion constant is calculated form the mean-square displacement of a molecule’s center of mass usingthe Einstein relation, fit to the range from 0.5 to 1.2 psof each 2.0 ps trajectory. Corrections for the finite size ofthe system were made for the diffusion constant, accord-ing to D=D(L) + 2.837 k B T/(6 πη L), where D(L) is thediffusion constant for simulation with side length L, k B is Boltzmann’s constant, and η is the viscosity.[36] Weuse the experimental value for the viscosity of 8.9 × − kg m − s − .[37] The NMR relaxation time, τ NMR , wasfound from the second order Legendre polynomial of thetime correlation function for rotations around the axisconnecting the hydrogen atoms.[38] We fit this to an ex-ponential function over the range from 1.0 to 1.5 ps.Voronoi tessellation was used to estimate the volumeoccupied by solvent and solute molecules. The Voronoimethod of assigning a region of space to a specific atomassumes that all atoms are of equal size. AlternativeVoronoi volumes can be defined which take into ac-count molecular size by assigning a molecule a radius,r A .[39, 40] In the surface-based (S-cell) method, a regionof space is assigned to an atom using the distance to thesurface of a sphere of radius R A around the atom. Ra-dial distribution functions were used to assign r A . Thewater radial distribution function (RDF) reaches one ata distance of about 2.6 ˚A, giving a water r A of 1.3 ˚A.The water-methane RDF reaches one at a distance of 3.5˚A, giving a methane r A of 2.2 ˚A, after subtracting thewater radius. Monte Carlo integration was used to cal-culate both the normal and the S-cell Voronoi volumesusing 10 random points per configuration. III. RESULTS
Water.
The minimized water monomer has a OH bondlength of 0.958 ˚A and a bond angle of 104.6 ◦ , very close tothe experimental values of 0.957 ˚A and 104.474 ◦ ,[41] andan energy of -2078.50 eV. At 300 K, the average monomerenergy is -2078.43 eV. The water dimer energies for theANI-1ccx model are compared to the ab initio results,at the CCSD(T) level in the CBS limit,[42] as shown inFigure 1(A).[43] Note that the range of the ANI-1ccx -7-6-5-4-3-2-10 d i m e r ene r g y ( kc a l / m o l ) pa i r ene r g y ( kc a l / m o l ) FIG. 1: (A) Minimized energy of the water pair as a functionof oxygen separation comparing the ANI-1ccx model to abinitio results[42]. (B) Minimized dimer energy (solid greenline) and average pair energy for the liquid at 300 K (dashedred line) for the TIP4P/2005 model.[44]. model is shorter than the ab initio results, going to zeroabove 5.2 ˚A, the cut-off of the potential. Results fromliquid phase simulations using pair potentials assess theimportance of the long-ranged interactions. Figure 1(B)
TABLE I: Cluster binding energies (in kcal/mol) forthe ANI-1ccx model, CCSD(T)/aug-cc-pVDZ[46], andTIP4P/2005[47].N ANI-1ccx CCSD(T) TIP4P/20052 -4.57 -5.29 -6.813 -15.64 -16.33 -18.404 -27.47 -28.54 -30.665 -36.22 -37.43 -39.996 prism -46.06 -48.82 -51.036 cage -45.87 -48.38 -52.116 book -45.87 -47.71 -50.786 cyclic -45.41 -46.17 -48.82 shows the dimer energy for TIP4P/2005 model[44] andlike other non-polarizable models it has a deeper mini-mum, since it is optimized to treat liquid phase inter-actions. Since the TIP4P/2005 model is pair-wise addi-tive, the interaction energies between pairs in the liquidphase can be determined. These results, shown by the reddashed line in Figure 1(B), show that beyond the nearest-neighbor distances the interaction between pairs is nearzero. At distances beyond 3.5 ˚A, the geometries betweenpairs are not optimal—liquid structure is largely deter-mined by optimizing nearest neighbors interactions—sothe average interactions are small. The longer-rangedinteractions are not as significant as the dimer energiessuggest.The water trimer was used to assess the size of non-additive contributions to the energy. For the minimumenergy cyclic structure, with three hydrogen bonds, theANI-1ccx potential gives a binding energy (the energy ofthe timer minus three times the energy of the monomer)of -15.64 kcal/mol, or -5.2 kcal/mol per molecule. Thedimer minimum is -4.57 kcal/mol, so the trimer formshydrogen bonds − ab initio methods find that thethree-body contribution to the energy is around 11%to 14%.[45] The minimum energies for clusters of up tosix molecules are given in Table I. For the N=6 clus-ter, energies for 4 local minima structures[46] are alsogiven. These are compared to ab initio results at theCCSD(T)/aug-cc-pVDZ level[46] and the TIP4P/2005model, using previously published results for N=2-5[47]and our own calculations for the hexamer structures. TheCCSD(T)/aug-cc-pVDZ results give a three-body contri-bution to the hydrogen bond energy ( δ E=(E(N=3) − − δ E is +0.68 kcal/mol.It cannot be negative, without many-body terms, andthe positive value reflects that in the trimer, the rigidTIP4P/2005 model cannot form completely optimal hy-drogen bonds between each pair. For the larger clusters,the agreement between ANI-1ccx and the CCSD(T) cal-culations is pretty good, with the NN model consistentlygiving higher energies. It correctly ranks the prism struc-ture as the lowest energy hexamer structure, although it
TABLE II: Properties of liquid water a temperature of 300 Kand a density of 1 g cm − .∆H vap P D τ NMR (kcal/mol) (kbar) (10 − m s − ) (ps)ANI-1ccx 11.8 ± ± ± ± a b c a. Reference [49] b. Reference [53] c. Reference [54] finds the cage and book structure to have equal ener-gies, rather than the cage being lower. The TIP4P/2005model, like other empirical liquid state models,[48] doesnot correctly give the lowest energy hexamer structure.The properties of liquid water as predicted by theANI-1ccx potential are given in Table II. The enthalpyof vaporization is calculated from ∆H vap = (cid:104) E (cid:105) gas −(cid:104) E (cid:105) liquid /N +RT, where (cid:104) E (cid:105) gas is the energy of an iso-lated gas-phase molecule, (cid:104) E (cid:105) liquid is the average energyof the liquid with N molecules and R is the ideal gasconstant. It is in pretty good agreement with the ex-perimental value of 10.5 kcal/mol,[49] although overesti-mated by about one kcal/mol. The pressure is negative,which, along with the overestimated ∆H vap , indicatesthat the balance between attractive and repulsive inter-actions is shifted towards attractive interactions. AIMDwith simple generalized gradient approximation (GGA)functionals tend to underestimate the density, consis-tent with an overestimated pressure at the experimentaldensity,[50, 51] The addition of dispersion interactionsleads to improvement.[17, 30, 50–52] More recent sim-ulations with meta-GGA functionals, which have betterexchange interactions, give larger densities, above the ex-perimental value.[32, 35] The ANI-1ccx potential givesgood agreement for both dynamical properties, the trans-lational diffusion constant and the rotatational constant, τ NMR , compared to experiment.[53, 54]Given that the potential has no charges and no elec-trostatic moments of any kind, the model does not havea dielectric response. The dielectric constant will be one.By enlarging the training set to include ions, a NN modelcould be trained to respond to the field of an ion. Thismight require extending the range of the potential, al-though to fully capture the dielectric response of wa-ter, methods that explicitly treat long-ranged interac-tions may be necessary.[55, 56]The radial distribution functions for water are shown inFigure 2. The oxygen-oxygen,[57] oxygen-hydrogen,[58]and hydrogen-hydrogen[58] experimental results are alsoshown. Also shown for comparison are the AIMD results,at CCD/aug-cc-pVDZ level, of Liu, et al. , which give verygood agreement with experiment.[34] The agreement forthe ANI-1ccx model, particularly for the oxygen-oxygeng(r), is very good for both peak positions and heights.The model does give a small probability for pairs at lowr values. This region below r < g OO (r) g O H (r) g HH (r) FIG. 2: Water-water radial distribution functions for theANI-1ccx model (red lines), AIMD[34] (blue dotted lines),X-ray[57] and neutron diffraction data[58, 60] (black dashedlines). pears a bit noisy in the g(r) since these close pairs arerelatively uncommon. In these structures, the hydrogenbetween the two oxygen atoms is shifted away from theline connecting the two oxygens, so that there is not alinear hydrogen bond. Other than the oxygen pairs be-ing very close, there is no atypical structures, like threewater rings, formed. The oxygen-hydrogen distance doesnot get short in this structure, so the radial distributionfunctions involving hydrogens do not show any unusualfeatures. The close structures result from an oxygen-oxygen potential that is too soft at small distances (Fig-ure 1(A)) and, like the pressure results, indicates that therepulsive interactions are underestimated. The potentialis only too soft at distances below 2.6 ˚A, so the oxygen-oxygen radial distribution function for most of the firstpeak, at distances larger than 2.6 ˚A, is unaffected. Forthe oxygen-hydrogen and hydrogen-hydrogen radial dis-tribution functions, the peaks are a little too structured,as are some of the results from AIMD.[30, 35, 59] OtherML potentials have given very accurate radial distribu-tion functions for water as well.[18, 19, 21]Nuclear quantum effects (NQE) from the quantum na-ture of the hydrogen atom lead to measurable struc-tural and thermodynamic differences between light andheavy water.[60, 61] In running classical dynamics using ab initio -derived potentials, these effects are neglected.Simulations using path integrals to include quantum ef-fects find that the inclusion of quantum effects decreasesthe heights and broadens with widths of the O-H andH-H RDFs,[21, 62, 63] but has a smaller effect on theO-O RDF.[21, 62, 64] Consistent with these results, neu-tron diffraction results for light and heavy water findthat heavy water is more structured that light water.[65]The negative pressure of the ANI-1ccx potential maybe partially due the classical treatment of hydrogensas well. NQEs have been shown to increase the den-sity of water,[21, 63, 66] and therefore increase the pres-sure. Other studies suggest that quantum effects de-crease the density.[67, 68] Experimentally, heavy waterhas smaller density (1.1044 g cm − [69] at 25 ◦ ) than wouldbe predicted using the density of light water (0.9971 gcm − [70]) and multiplying by the mass ratio (20/18), togive 1.1079 g cm − , implying that NQEs increase thedensity. Path integral simulations could be run with theANI-1ccx potentials, but keeping with more typical, andless costly, simulation methods, purely classical simula-tions were used.Hydrogen bonds are defined using the geometric defi-nition of Luzar and Chandler,[71] in which the oxygen-oxygen distance is less that 3.5 ˚A and the angle, β , is lessthan 30 ◦ , where β is the angle between the vector con-necting to the hydrogen bond donating oxygen and thehydrogen bond accepting oxygen and the vector connect-ing the donating oxygen to the hydrogen involved in thehydrogen bond. The model gives an average number ofhydrogen bonds of 3.65 ± TABLE III: The average number of hydrogen bonds and theaverage tetrahedral order parameter for liquid water at 300K and 1 g cm − . hydrogen bonds QANI-1ccx 3.65 ± ± a b AIMD 3.60 − c − d TIP4P/2005 3.66 0.668a. Reference [72] b. Reference [65] c. Reference [33]d. References [33, 35] diffraction data[72] gives an average hydrogen bond num-ber of 3.58 and AIMD with various DFT methods givehydrogen bond numbers of 3.60 to 3.85,[33] all using thesame hydrogen bond definition as used here. The re-sults for the distribution of β for neighbors within thefirst solvation shell, at a distance less than the minimumof the oxygen-oxygen pair correlation function (3.4 ˚A),are shown in Figure 3(A). Combining measurements of P ( β ) β (degrees)ANI-1ccxexperimentTIP4P/2005(A)00.511.522.5 P ( Q t ) t ANI-1ccxTIP4P/2005(B)
FIG. 3: (A) Distribution of the hydrogen bond angles forliquid water from the ANI-1ccx model (red line), the abinitio /experimental result[73] (blue dashed line) and theTIP4P/2005 model (green dotted line). (B) Distribution ofthe tetrahedral order parameter for liquid water from theANI-1ccx model (red line) and the TIP4P/2005 model (greendashed line). the proton magnetic shielding tensor with DFT calula-tions, Modig, Pfrommer, and Halle constructed distribu-tion functions of the hydrogen bond angle, β .[73] Ourresults are in pretty close agreement with those results,except for a tail of the distribution at angles beyond40 ◦ , indicating nearest-neighbors which are not hydro-gen bonded. For TIP4P/2005 model, this tail is larger,with a second peak, near 50 ◦ . AIMD results[30] also findthis tail, while the NN potential of Behler and co-workersdoes not.[17]The structure of a water molecule’s nearest neighborscan be characterized using the tetrahedral order param-eter, Q, of Errington and Debenedetti[74] Q i = 1 − (cid:88) j =1 4 (cid:88) k = j +1 (cid:18) cos θ ijk + 13 (cid:19) (4)where the sums are over the four nearest neighborsof molecule i and θ ijk is the angle between the threemolecules, with molecule i at the center. This parameteris zero for a random arrangement of neighbors, as in anideal gas, and one for a perfectly tetrahedral structure, asin ice Ih. The distribution of Q is shown in Figure 3(B),with the TIP5P/2005 results for comparison. The aver-age value of Q is 0.66 ± Methane in water.
The water-methane dimer energyis shown in Figure 4. The agreement with the ab initio results[42] are very good. The properties of the methanewater solution are shown in Figure 5. Integrating theoxygen-carbon RDF over this first solvation shell, out to5.4 ˚A, gives 20 nearest neighbors. Simulation studies findcoordinate numbers of 19 to 23,[75–77] with AIMD giving17.[78] Experimental estimates from neutron diffractionresults give 16[79] and 19.[80] This value is sensitive tothe cut-off distance, the AIMD result of Montagna, etal. ,[78] and the experimental results of Koh, et al. ,[79]which give lower coordination numbers, use a 5.0 ˚A limit.The others use a limit of around 5.5 ˚A. Using a 5.0 ˚A limitwith the ANI-1ccx results gives a coordination numberof 17.From Figure 5 it is evident that the structure of waterin the first solvation shell of methane, from about 3.5 ˚Ato 5.4 ˚A, is different from the bulk. The number of hydro-gen bonds increases to over 3.9 and the tetrahedral orderincreases as well. The volume available to a molecule as -2-10123456 d i m e r ene r g y ( kc a l / m o l ) FIG. 4: Minimized energy of the water-methane pair as afunction of oxygen-carbon separation comparing the ANI-1ccx model to ab initio results[42].. given by Voronoi volume increases if the standard tessel-lation is used, but decreases if the S-cell method is used,which takes into account different molecular size, as hasbeen seen previously.[81] A number of simulations usingAIMD,[78, 82] all-atom force fields,[83–87] and coarse-grained models[81] and experimental studies[82, 88–90]have examined the structure of water around hydropho-bic molecules. The results of these studies are mixed,with some finding a more[78, 81–85, 87, 89, 90] and somefinding a less[86, 88] structured first solvation shell. Ourresults, finding an increased number of hydrogen bondsand an increased tetrahedral order, imply more orienta-tional order of water next to methane. A number of sim-ulation studies have found enhanced tetrahedral orderin the first solvation shell.[82, 85, 87] AIMD[78, 82] andall-atom potential model simulations[84, 87] also find anincreased number of hydrogen bonds, while other modelsfind no difference in hydrogen bonds.[85, 86]
IV. CONCLUSIONS
The NI-1ccx potential reproduces the many ofthe structural properties of liquid water and themethane/water solution, as well as dynamical proper-ties of liquid water, at the standard density and 300 K.It should be stressed that the model was not developedspecifically for water, liquid or otherwise, but as partof a broadly applicable potential for molecules contain-ing H,C,N, and O atoms, trained on quantum data formolecules.[9] No experimental or theoretical data fromthe liquid phase was used in the training of the model.For water, the potential reproduces dynamical proper-ties for both translational and rotational motion (Ta- g C O (r) (A) 01020304050 N (r) h y d r ogenbond s (B)2 4 6 8 10r ( ˚A)0.620.640.660.680.7 Q t (r) (C) 2 4 6 8 10r ( ˚A)242628303234 V o r ono i v o l u m e (˚ A ) (D) FIG. 5: (A) Carbon-oxygen radial distribution function (sold line) and running coordination number, N(r), (blue dashed line)(B) hydrogen bond numbers, (C) tetrahedral order parameter, and (D) standard (solid line) and S-cel (blue dashed line) Voronoivolumes as a function of carbon-oxygen distance. ble II). The structure of the liquid as given by the ra-dial distribution functions (Figure 2) are in good agree-ment with experiment, as good or better than many pair-potentials and AIMD simulations.[30, 91] Many-body po-larization effects are also described accurately. The hy-drogen bonded trimer has hydrogen bonds about 12%stronger than those of the dimer, in aggreement of the abinitio results.[45] Non-additive effects are also evident inthe liquid phase energy, as given by ∆H vap , compared tothe dimer energy. The ANI-1ccx potential has the inter-action strengths correct for both the dimer and the liq-uid phase. Non-polarizable potentials, like TIP4P/2005,have to overestimate the pair energy to get an accu-rate ∆H vap (Figure 1 (B)). For the methane/water so-lution, the model gives an accurate number of nearest-neighbor water molecules and an increase in the numberof hydrogen bonds and tetrahedral order for the solvat-ing waters, in agreement with experimental and AIMD results.[78, 79, 82] The one main inaccuracy of the modelis the negative pressure for water at the experimentaldensity, indicating the predicted density under standardconditions would be too high.The quality of the results is dependent on three factors:the mathematical form, the training, and the execution ofthe model. The one major approximation of the form ofthe potential is that it is short-ranged, going to zero at adistance of 5.2 ˚A. It also ignores polarization beyond thethree-body level. The training, using single molecules ordimers with energy minimized structures and structuresfrom normal mode sampling to generate structures awayfrom minima, appears to do well near the energy minimabut not at close distances (Figures 1 and 4). For theexecution of the model, all atoms were treated classically,so NQEs are neglected. The three factors then lead tothree specific aspects of the potential as applied to liquidwater: no interactions beyond 5.2 ˚A, too soft at smalldistances, and no NQEs.The influence of NQEs on the properties of water arewell-known. As discussed above, treating the hydrogenatoms as quantum particles would increase the density(and pressure)[21, 63, 66] and broaden the widths of theoxygen-hydrogen and hydrogen-hydrogen radial distribu-tion functions, without changing the oxygen-oxygen verymuch.[21, 62–64] These would all improve the results ofthe ANI-1ccx potential, but the effects are not all thatlarge for the density change. Any deficiencies of themodel are not likely to be due to a lack of long-rangedinteractions. For one thing, addition of long-ranged in-teractions would only make the pressure more negative,correcting it in the wrong direction. Secondly, the liquidinteractions quickly become small above 3.5 ˚A (Figure1(B)), and accurate NN models for water have been de-veloped that do not have longed-ranged interactions.[17](Even NN models for protonated water clusters were de-veloped which did not have charges.[92, 93]) Long-rangedinteractions could be built into a NN model by addingCoulombic interactions with either assigned[19] or NNgenerated[22, 56] charges, but this does not appear tobe necessary, at least for the systems studied here. Thenegative pressure, the most inaccurate result of the ANI-1ccx model for the water simulations, is likely due to thepotential being not repulsive enough at small distances.This could be corrected by training with dimer structures at these distances. The ANI-1ccx potential was trainedfor dimers only near the energy minimum.[9, 12] Alterna-tively, corrections to the NN potential could be done asthe simulation proceeds, on-the-fly, using query by com-mittee methods[9] or other[7] methods to determine whenthe model may be inaccurate.The success of the ANI-1ccx potential in simulating achallenging system like water demonstrates the promiseof the approach. Using carefully selected data sets ofonly single molecules and molecule pairs, models can betrained which work well in the liquid phase. This en-ables the study of novel materials for which there maynot be experimental data, as well as simulating react-ing systems. The present results indicate that improvedmodels could be trained by using dimer structures thatprobe the repulsive part of the potential along with thepotential minimum.
Acknowledgements
This material is based upon work supported by theU.S. Department of Energy, Office of Science, Basic En-ergy Sciences, under EPSCoR Grant No. DE-SC0012432with additional support from the Louisiana Board of Re-gents. [1] R. Car, M. Parrinello, Unified approach for moleculardynamics and density-functional theory, Phys. Rev. Lett.55 (1985) 2471–2474.[2] E. Paquet, H. L. Viktor, Computational methods for abinitio molecular dynamics, Advances in Chemistry 2018(2018) 9839641.[3] T. B. Blank, S. D. Brown, A. W. Calhoun, D. J. Doren,Neural network models of potential energy surfaces, J.Chem. Phys. 103 (10) (1995) 4129–4137.[4] S. Manzhos, X. Wang, R. Dawes, T. Carrington, A nestedmolecule-independent neural network approach for high-quality potential fits, J. Phys. Chem. A 110 (16) (2006)5295–5304.[5] J. Behler, M. Parrinello, Generalized neural-networkrepresentation of high-dimensional potential-energy sur-faces, Phys. Rev. Lett. 98 (2007) 146401.[6] O. A. von Lilienfeld, R. Ramakrishnan, M. Rupp,A. Knoll, Fourier series of atomic radial distributionfunctions: A molecular fingerprint for machine learningmodels of quantum chemical properties, Int. J. QuantumChem. 115 (16) (2015) 1084–1093.[7] Z. Li, J. R. Kermode, A. De Vita, Molecular dynamicswith on-the-fly machine learning of quantum-mechanicalforces, Phys. Rev. Lett. 114 (2015) 096405.[8] J. Behler, Perspective: Machine learning potentials foratomistic simulations, J. Chem. Phys. 145 (2016) 170901.[9] J. S. Smith, B. Nebgen, N. Lubbers, O. Isayev, A. E.Roitberg, Less is more: Sampling chemical space withactive learning, J. Chem. Phys. 148 (2018) 241733.[10] J. S. Smith, B. T. Nebgen, R. Zubatyuk, N. Lubbers, C. Devereux, K. Barros, S. Tretiak, O. Isayev, A. E.Roitberg, Approaching coupled cluster accuracy with ageneral-purpose neural network potential through trans-fer learning, Nat. Commun. 10 (2019) 2903.[11] A. P. Bart´ok, R. Kondor, G. Cs´anyi, On representingchemical environments, Phys. Rev. B 87 (2013) 184115.[12] J. S. Smith, O. Isayev, A. E. Roitberg, Ani-1: an exten-sible neural network potential with dft accuracy at forcefield computational cost, Chem. Sci. 8 (2017) 3192–3203.[13] S. Chmiela, H. E. Sauceda, K. M¨uller, A. Tkatchenko,Towards exact molecular dynamics simulations withmachine-learned force fields, Nat. Commun. 9 (2018)3887.[14] M. J. Willatt, F. Musil, M. Ceriotti, Atom-density rep-resentations for machine learning, J. Chem. Phys. 150(2019) 154110.[15] S. Liu, J. Li, K. C. Bennett, B. Ganoe, T. Stauch,M. Head-Gordon, A. Hexemer, D. Ushizima, T. Head-Gordon, Multiresolution 3d-densenet for chemical shiftprediction in nmr crystallography, J. Phys. Chem. Lett.10 (16) (2019) 4558–4565.[16] J. Behler, First principles neural network potentials forreactive simulations of large molecular and condensedsystems, Angew. Chem. Int. Ed. Engl. 56 (42) (2017)12828–12840.[17] T. Morawietz, A. Singraber, C. Dellago, J. Behler, Howvan der waals interactions determine the unique proper-ties of water, Proc. Natl. Acad. Sci. (USA) 113 (2016)8368–8373.[18] S. K. Natarajan, J. Behler, Neural network molecular dynamics simulations of solid-liquid interfaces: water atlow-index copper surfaces, Phys. Chem. Chem. Phys. 18(2016) 28704.[19] H. Wang, W. Yang, Force field for water based on neuralnetwork, J. Phys. Chem. Lett. 9 (2018) 3232–3240.[20] A. Singraber, T. Morawietz, J. Behler, C. Dellago, Den-sity anomaly of water at negative pressures from firstprinciples, J. Phys.: Condens. Matter 30 (2018) 254005.[21] B. Cheng, E. A. Engel, J. Behler, C. Dellago, M. Ceriotti,Ab initio thermodynamics of liquid and solid water, Proc.Natl. Acad. Sci. (USA) 116 (2019) 1110–1115.[22] T. Morawietz, V. Sharma, J. Behler, A neural networkpotential-energy surface for the water dimer based onenvironment-dependent atomic energies and charges, J.Chem. Phys. 136 (2012) 064103.[23] Y. Li, H. Li, F. C. Packard IV, B. Narayanan, F. Sen,M. K. Y. Chan, S. Sankaranarayanan, B. R. Brooks,B. Roux, Machine learning force field parameters from abinitio data, J. Chem. Theory Comput. 13 (2017) 4492–4503.[24] T. Bereau, R. A. DiStasio, A. Tkatchenko, O. A. vonLilienfeld, Non-covalent interactions across organic andbiological subsets of chemical space: Physics-based po-tentials parametrized from machine learning, J. Chem.Phys. 148 (24) (2018) 241706.[25] T. T. Nguyen, E. Sz´ekely, G. Imbalzano, J. Behler,G. Cs´anyi, M. Ceriotti, A. W. G¨otz, F. Paesani, Com-parison of permutationally invariant polynomials, neuralnetworks, and gaussian approximation potentials in rep-resenting water interactions through many-body expan-sions, J. Chem. Phys. 148 (24) (2018) 241725.[26] M. E. Taylor, P. Stone, Transfer learning for reinforce-ment learning domains: a survey, J. Mac. Learn. Res. 10(2009) 1633–1685.[27] C. Riplinger, P. Pinski, U. Becker, E. F. Valeev, F. Neese,Sparse maps-a systematic infrastructure for reduced-scaling electronic structure methods. ii. liner scaling do-main based pair natural orbital couple cluster theory, J.Chem. Phys. 144 (2016) 024109.[28] A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E.Castelli, R. Christensen, M. Du(cid:32)lak, J. Friis, M. N.Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C.Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L.Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B.Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Pe-terson, C. Rostgaard, J. Schiøtz, O. Sch¨utt, M. Strange,K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter,Z. Zeng, K. W. Jacobsen, The atomic simulation environ-ment—a python library for working with atoms, Journalof Physics: Condensed Matter 29 (27) (2017) 273002.[29] H. J. C. Berendsen, J. R. Grigera, T. P. Straatsma, Themissing term in effective pair potentials, J. Phys. Chem.91 (1987) 6269.[30] Z. Ma, Y. Zhang, M. E. Tuckerman, Ab initio moleculardynamics study of water at constant pressure using con-verged basis sets and empirical dispersion corrections, J.Chem. Phys. 137 (2012) 044506.[31] S. Y. Willow, M. A. Salim, K. S. Kim, S. Hirata, Ab ini-tio molecular dynamics of liquid water using embedded-fragment second-order many-body perturbation theorytowards its accurate property prediction, Sci. Rep. 5(2015) 14358.[32] L. R. Pestana, N. Mardirossian, M. Head-Gordon,T. Head-Gordon, Ab initio molecular dynamics simula- tions of liquid water using high quality meta-gga func-tionals, Chem. Sci. 8 (2017) 3554.[33] L. Zheng, M. Chen, Z. Sun, H.-Y. Ko, B. Santra, P. Dhu-vad, X. Wu, Structural, electronic, and dynamical prop-erties of liquid water by ab initio molecular dynamicsbased on scan functional within the canonical ensemble,J. Chem. Phys. 148 (16) (2018) 164505.[34] J. Liu, X. He, J. Z. H. Zhang, L. W. Qi, Hydrogen-bondstructure dynamics in bulk water: insights from ab ini-tio simulations with coupled cluster theory, Chem. Sci. 9(2018) 2065–2073.[35] M. Chen, H.-Y. Ko, R. C. Remsing, M. F. Calegari An-drade, B. Santra, Z. Sun, A. Selloni, R. Car, M. L. Klein,J. P. Perdew, X. Wu, Ab initio theory and modelingof water, Proc. Natl. Acad. Sci. (USA) 114 (41) (2017)10846–10851.[36] I. C. Yeh, G. Hummer, System-size dependence of diffu-sion coefficients and viscosities from molecular dynamicssimulations with periodic boundary conditions, J. Phys.Chem. B 108 (2004) 15873–15879.[37] K. R. Harris, L. A. Woolf, J. Chem. Eng. Data 49 (2004)1064.[38] R. W. Impey, P. A. Madden, I. R. McDonald, Spectro-scopic and transport properties of water. model calcu-lations and the interpretation of experimental results.,Mol. Phys. 46 (1982) 513–539.[39] S. V. Anishchik, N. N. Medvedev, Three-dimensionalapollonian packing as a model for dense granular sys-tems, Phys. Rev. Lett. 75 (1995) 4314–4317.[40] V. P. Voloshin, N. N. Medvedev, M. N. Andrews, R. R.Burri, R. Winter, A. Geiger, Volumetric properties ofhydrated peptides: Voronoi–delaunay analysis of molec-ular simulation runs, J. Phys. Chem. B 115 (48) (2011)14217–14228.[41] F. Franks, in: F. Franks (Ed.), Water-A ComprehensiveTreatise, Plenum, New York, 1972, pp. 1–20.[42] M. P. Metz, K. Szalewicz, J. Sarka, R. T´obi´as, A. G.Cs´asz´ar, E. M´atyus, Molecular dimers of methaneclathrates: ab initio potential energy surfaces and vari-ational vibrational states, Phys. Chem. Chem. Phys. 21(2019) 13504.[43] E. B. Guidez, M. S. Gordon, Dispersion correction de-rived from first principles for density functional theoryand hartree-fock theory, J. Phys. Chem. A 119 (2015)2161–2168.[44] J. L. F. Abascal, C. Vega, A general purpose model forthe condensed phases of water: TIP4p/2005, J. Chem.Phys. 123 (2005) 234505.[45] F. N. Keutsch, J. D. Cruzan, R. J. Saykally, The watertrimer, Chem. Rev. 103 (2003) 2533–2577.[46] E. Miliordos, E. Apr´a, S. S. Xantheas, Optimal geome-tries and harmonic vibrational frequencies of the globalminima of water clusters (h2o)n, n = 2–6, and severalhexamer local minima at the ccsd(t) level of theory, J.Chem. Phys. 139 (11) (2013) 114302.[47] B. S. Gonz´alez, E. G. Noya, C. Vega, L. M. Ses´e, Nuclearquantum effects in water clusters: The role of molecularflexibility, J. Phys. Chem. B 114 (2010) 2484–2492.[48] T. James, D. J. Wales, J. Hern´andez-Rojas, Global min-ima of water clusters, [h2o]n, n <
21, described by a five-site empirical potential, Chem. Phys. Lett. 415 (2005)302–307.[49] N. E. Dorsey, Properties of Ordinary Water-Substancein all its Phases: Water Vapor, Water, and all the Ices, Reinhold Publishing, New York, 1940.[50] J. Schmidt, J. VandeVondele, I.-F. W. Kuo, D. Sebas-tiani, J. I. Siepmann, J. Hutter, C. J. Mundy, Isobaric-isothermal molecular dynamics simulations utilizing den-sity functional theory: An assessment of the structureand density of water at near-ambient conditions, J. Phys.Chem. B 113 (35) (2009) 11959–11964.[51] A. P. Gaiduk, F. Gygi, G. Galli, Density and compress-ibility of liquid water and ice from first-principles sim-ulations with hybrid functionals, J. Phys. Chem. Lett.6 (15) (2015) 2902–2908.[52] S. Y. Willow, X. C. Zeng, S. S. Xantheas, K. S. Kim,S. Hirata, Why is mp2-water “cooler” and “denser” thandft-water?, J. Phys. Chem. Lett. 7 (4) (2016) 680–684.[53] K. Krynicki, C. D. Green, D. W. Sawyer, Pressure andtemperature dependence of self-diffusion in water, Dis-cuss. Faraday Soc. 66 (1978) 199–208.[54] J. Jonas, T. DeFries, D. J. Wilbur, Molecular motionsin compressed liquid water, J. Chem. Phys. 65 (1976)582–588.[55] A. Grisafi, D. M. Wilkins, G. Cs´anyi, M. Ceriotti,Symmetry-adapted machine learning for tensorial prop-erties of atomistic systems, Phys. Rev. Lett. 120 (2018)036002.[56] A. Grisafi, M. Ceriotti, Incorporating long-range physicsin atomic-scale machine learning, J. Chem. Phys. 151(2019) 204105.[57] L. B. Skinner, C. J. Benmore, J. C. Neuefeind, J. B.Parise, The structure of water around the compressibilityminumum, J. Chem. Phys. 141 (2014) 214507.[58] A. K. Soper, The radial distribution functions of waterand ice from 220 to 673 k and at pressures up to 400mpa, Chem. Phys. 258 (2000) 121–137.[59] K. Forster-Tonigold, A. Groß, Dispersion corrected rpbestudies of liquid water, J. Chem. Phys. 141 (6) (2014)064501.[60] A. K. Soper, The radial distribution functions of wateras derived from radiation total scattering experiments: Isthere anything we can say for sure?, ISRN Phys. Chem.2013 (2013) 279463.[61] S. Herrig, M. Thol, A. H. Harvey, E. W. Lemmon, A ref-erence equation of state for heavy water, J. Phys. Chem.Ref. Data 47 (4) (2018) 043102.[62] H. A. Stern, B. J. Berne, Quantum effects in liquid water:Path-integral simulations of a flexible and polarizable abinitio model, J. Chem. Phys. 115 (16) (2001) 7622–7628.[63] T. Spura, C. John, S. Habershon, T. D. K¨uhne, Nu-clear quantum effects in liquid water from path-integralsimulations using an ab initio force-matching approach,Molecular Physics 113 (8) (2015) 808–822.[64] O. Marsalek, T. E. Markland, Quantum dynamics andspectroscopy of ab initio liquid water: The interplay ofnuclear and electronic quantum effects, J. Phys. Chem.Lett. 8 (2017) 1545–1551.[65] A. K. Soper, C. J. Benmore, Quantum differences be-tween heavy and light water, Phys. Rev. Lett. 101 (2008)065502.[66] G. S. Fanourgakis, S. S. Xantheas, Development oftransferable interaction potentials for water. v. extensionof the flexible, polarizable, thole-type model potential(ttm3-f, v. 3.0) to describe the vibrational spectra of wa-ter clusters and liquid water, J. Chem. Phys. 128 (2008)074506.[67] G. S. Fanourgakis, S. S. Xantheas, The flexible, polariz- able, thole-type interaction potential for water (ttm2-f)revisited, J. Phys. Chem. A 110 (11) (2006) 4100–4106.[68] F. Paesani, S. Iuchi, G. A. Voth, Quantum effects in liq-uid water from an ab initio-based polarizable force field,J. Chem. Phys. 127 (2007) 074506.[69] M. Nakamura, K. Tamura, S. Murakami, Isotope effectson thermodynamic properties: mixtures of x(d2o or h2o)+ (1 - x)ch3cn at 298.15 k, Thermochimica Acta 253(1995) 127 – 136.[70] G. S. Kell, Density, thermal expansivity, and compress-ibility of liquid water from 0 ◦ to 150 ◦ c: Correlations andtable for atmospheric pressure and saturation reviewedand expressed on 1968 temperature scale, J. Chem. Eng.Data 20 (1975) 97–105.[71] A. Luzar, D. Chandler, Hydrogen-bond kinetic in liquidwater, Nature 379 (1996) 55–57.[72] A. K. Soper, F. Bruni, M. A. Ricci, Site-site pair cor-relation functions of water from 25 to 400 ◦ c: Revisedanalysis of new and old diffraction data, J. Chem. Phys.106 (1) (1997) 247–254.[73] K. Modig, B. G. Pfrommer, B. Halle, Temperature-dependent hydrogen-bond geometry in liquid water,Phys. Rev. Lett. 90 (2003) 075502.[74] J. R. Errington, P. G. Debenedetti, Relationship betweenstructural order and the anomalies of liquid water, Na-ture 409 (2001) 318–321.[75] L. Lue, D. Blankschtein, Liquid-state theory ofhydrocarbon-water systems: application to methane,ethane, and propane, J. Chem. Phys. 96 (1992) 8582.[76] C. H. Bridgeman, A. D. Buckingham, N. T. Skipper,Temperature dependence of solvent structure around ahydrophobic solute: a monte carlo study of methane inwater, Chem. Phys. Lett. 253 (1996) 209–215.[77] E. C. Meng, P. A. Kollman, Molecular dynamics studiesof the properties of water around simple organic solutes,J. Phys. Chem. 100 (1996) 11460.[78] M. Montagna, F. Sterpone, L. Guidoni, Structural andspectroscopic properties of water around small hydropho-bic solutes, J. Phys. Chem. B 116 (2012) 11695–11700.[79] C. A. Koh, R. P. Wisbey, X. Wu, R. E. Westacott, A. K.Soper, Water ordering around methane during hydrateformation, J. Chem. Phys. 113 (2000) 6390.[80] P. H. K. de Jong, J. E. Wilson, G. W. Neilson, A. D.Buckingham, Hydrophobic hydration of methane, Mol.Phys. 91 (1) (1997) 99–104.[81] N. Islam, M. Flint, S. W. Rick, Water hydrogen degreesof freedom and the hydrophobic effect, J. Chem. Phys.150 (2019) 014502.[82] J. Grdadolnik, F. Merzel, F. Avbelj, Origin of hydropho-bicity and enhanced water hydrogen bond strength nearpurely hydrophobic solutes, Proc. Natl. Acad. Sci. (USA)114 (2) (2017) 322–327.[83] T. M. Raschke, M. Levitt, Nonpolar solutes enhance wa-ter structure within hydration shells while reducing in-teractions between them, Proc. Natl. Acad. Sci. (USA)102 (19) (2005) 6777–6782.[84] A. Godec, J. C. Smith, F. Merzel, Increase of both orderand disorder in the first hydration shell with increasingsolute polarity, Phys. Rev. Lett. 107 (2011) 267801.[85] N. Galamba, Water’s structure around hydrophobic so-lutes and the iceberg model, J. Phys. Chem. B 117 (7)(2013) 2153–2159.[86] J. Kim, Y. Tian, J. Wu, Thermodynamic and struc-tural evidence for reduced hydrogen bonding among wa-1