{'fraction_non_alphanumeric': 0.03007191109174112, 'fraction_numerical': 0.010852037480932665, 'mean_word_length': 4.622641509433962, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Minimal Experimental Bias on the Hydrogen Bond Greatly Improves Ab Initio Molecular Dynamics Simulations of Water.

Experiment Directed Simulations (EDS) is a method within a class of techniques seeking to improve molecular simulations by minimally biasing the system Hamiltonian to reproduce certain experimental observables. In a previous application of EDS to ab initio molecular dynamics (AIMD) simulation-based on electronic density functional theory (DFT), the AIMD simulations of water were biased to reproduce its experimentally derived solvation structure. In particular, by solely biasing the O-O pair correlation functions, other structural and dynamical properties that were not biased were improved. In this work, the hypothesis is tested that directly biasing the OH pair correlation, will provide an even better improvement of DFT-based water properties in AIMD simulations. The logic behind this hypothesis is that for most electronic DFT descriptions of water the hydrogen bonding is known to be deficient due to anomalous charge transfer and over polarization in the DFT. Using recent advances to the EDS learning algorithm, we thus train a minimal bias on AIMD water that reproduces the O-H radial distribution function derived from the highly accurate MB-pol model of water. It is then confirmed that biasing the O-H pair correlation alone can lead to improved AIMD water properties, with structural and dynamical properties in even closer to experiment than the previous EDS-AIMD model.


Introduction
Ab initio molecular dynamics (AIMD) simulation 1 has become a popular tool to understand water (see, e.g., refs [2][3][4][5] ) and aqueous solutions/environments (see, e.g., refs 6-7 ). By using AIMD, as opposed to classical empirical water models, electronic structure calculations -primarily within the Kohn-Sham density functional theory (DFT) framework -can be used to solve for the electronic ground state. Electron densities can then respond to the surrounding electric field and thus account for polarization of the electron cloud, and forces can be calculated on-the-fly to propagate the dynamics of the nuclei. Furthermore, AIMD becomes particularly valuable for simulating systems while accounting for chemical reactivity, as it does not require specifying a defined bonding topology, and as such allows one to study dynamics of aqueous systems containing excess protons or hydroxide ions (see, e.g., refs [8][9][10][11][12][13] ).
Accurately modeling the hydrogen bond between water molecules has proven challenging for AIMD water simulations that use generalized gradient approximations (GGA) such as PBE 14 or BLYP [15][16] as the exchange-correlation functional in the DFT. 17 GGA functionals result in overpolarization due to having a small energy gap between Kohn-Shan virtual orbitals and occupied orbitals, [18][19] and they also exhibit partial covalency (or charge-transfer) within the intermolecular interactions of water molecules. 17,[20][21] These inaccuracies in modeling the hydrogen bond result in a water model that is overstructured at room temperature and diffusion coefficients that can be an order or orders of magnitude slower than found in experiment. 9,[22][23] It is common practice in the AIMD community to mitigate the deficiencies of GGA functionals by increasing the simulation temperature or by going beyond GGA functionals. 2-3, 5, 24 It has previously been shown for the PBE functional that simulation temperatures of 400 K are necessary to mimic room temperature structure and dynamics of water. 23,[25][26] It has been the case that increasing the temperature is not only an ad hoc remedy for the glassy behavior of GGA functionals, but additionally used to mimic nuclear quantum effects in GGA AIMD in an ad hoc fashion. 27 By contrast, going beyond simple GGA functions and using hybrid or meta-GGA functionals has been shown to improve the hydrogen bond properties and hence the water properties, 2, 24 but doing so significantly increases the computational cost. It should also be noted that it is still commonplace to combine non-GGA functionals and increased simulation temperatures in AIMD as an ad hoc fix for nuclear quantum effects in room temperature water. [2][3] One way to improve the accuracy of MD simulations (including AIMD) is to add a biasing potential to an observable of the system in order to improve agreement between that observable and one in target system. Pitera and Chodera showed using the method of Lagrange multipliers that there exists a linear bias on observables that minimizes the relative entropy of an ideal probability distribution for an observable to a target experimental one. 28 Using this idea, White and Voth developed the Experiment Directed Simulation (EDS) method as a means to parameterize these unknown linear terms in the Hamiltonian on-the-fly within a single MD simulation using a stochastic gradient descent algorithm. 29 Since then, there have been additional advances to EDS and similar methods, 30-31 and these have been recently reviewed by Amirkulova and White. 32 In related work, the Coarse-Grained Directed Simulation (CGDS) method was developed where coarse-grained observables in a molecular sub-system are biased based on data from coarse-grained simulations of a larger supramolecular complex. 33 In this work, the stochastic gradient descent algorithm of the original EDS method was unable to find coupling constants fast enough due to longer timescale and limited configurations of the CG observables that are otherwise available in isotropic systems such as water. In the CGDS paper, however, several variants of the learning algorithm were developed, where the stochastic gradient descent algorithm was replaced with a gradient descent using the covariance of all deviations of collective variables from target observables, and by a Levenberg-Marquadt minimization. In the studies presented here, we will take advantage of these more efficient methods which have better convergence properties.
In previous work, 34 we developed the idea that EDS could be applied to adjust the solvation structure of a relatively low cost DFT AIMD approach, rather than employing the previously mentioned alternatives of increasing the simulation temperature or going beyond GGA functionals.
In that work 34 the solvation structure of BLYP and dispersion corrected BLYP water was biased to reproduce the experimental O-O radial distribution function (RDF). This EDS bias, henceforth known as "EDS-AIMD(OO)", was able to improve the targeted solvation structure of water without needing an increase of simulation temperature or computational cost. Importantly, other properties that were not biased also improved as a consequence of having a better O-O structure.
For example, when adding an excess proton to the system, the EDS-AIMD(OO) approach improved the ratio of the hydrated excess proton to water diffusion coefficient and did not disrupt the other properties of the system. To confirm that employing EDS to a more accurate model would not adversely affect our results, the method was also tested on a dispersion corrected functional. It was found that, as expected, the linear bias learned from the gradient descent algorithm for the BLYP-D3 system was smaller than the EDS bias of normal BLYP due to the already better agreement with experiment for BYTP-D3, and all results for the biased EDS-BLYP-D3 model were equivalent or better.
The present work builds on this previous study, but this time by biasing the hydrogen bond (intermolecular O-H coordination properties) as opposed to the oxygen coordination structure of water. We chose to pursue this approach due to the intuitive feeling that the having a correct description of hydrogen bonding in DFT water is the underlying deficiency leading to other poor properties, due to the over-polarization and anomalous charge transfer of GGA functionals.
Correcting the hydrogen bond structure with EDS is more difficult than the O-O coordination, due to the anisotropic nature of hydrogen bonding, and the need to distinguish an O-H pair as either covalently bonded or hydrogen bonded. By biasing the hydrogen bond, however, we can more directly target the charge transfer and over-polarization found in GGA water simulations which, among other things, gives rise to its glassy behavior and slow diffusion. In this work, we thus bias BLYP and BLYP-D3 simulations to reproduce the O-H RDF of the highly and demonstrably accurate MB-pol water model, [35][36][37][38] which includes many-body interactions parameterized from high level CCSD(T) electronic structure calculations and many-body polarization effects. MB-pol accurately reproduces many properties of water and is gaining recognition as one of the most -if not the most -accurate water models available. By using the classical MB-pol O-H RDF as the refence in EDS-AIMD, we are therefore able to effectively include higher order correlations into the BLYP-level AIMD water simulations.
The remainder of this paper is organized as follows: In the Methods section, we first include a subsection that reviews the theory regarding the EDS method and how we apply it to bias the hydrogen bond. In the next subsection of Methods, we describe the details of the EDS algorithm as used, while in the Simulations Procedure subsection, we describe the simulation setup and details for the EDS learning algorithm and the production runs. In the Results and Discussion section, we then present the results and give further discussion. We lastly present conclusions and the outlook for the use of experimental-bias methods in AIMD.

Experiment Directed Simulation Method
EDS modifies the system Hamiltonian by adding a biasing potential that is parameterized to reproduce a target observable. 29 In the present case, the hydrogen bond is biased to reproduce the O-H radial distribution function (RDF) of the MB-pol water model at 298 K. As in the previous work, the RDF is corrected by biasing its statistical moments such as the coordination number (the zeroth moment). The bias on each moment is a function of three variables: a coupling constant, α ! , a function characterizing the pairwise distances between all hydrogen-oxygen pairs, " ( ⃗ #$ ), and its target average value, " ' . The EDS potential for each hydrogen atom is therefore: where M represents the number of moments being biased, , is the number of oxygen atoms in the system, and we take " ( ⃗ #$ ) to be a product of the O-H pairwise distance and mollified step function. The step function for biasing the hydrogen bond was chosen as: where ) was set to 2.125 Å, was set to 0.7, and 0 was set to 1.
2, and 34 . %& / is the target O-H radial distribution function. Note that by k th moment we are specifically referring to the power in Eq. 3, and the zeroth moment corresponds to the coordination number.

Determining EDS Coupling Constants
The determination of EDS coupling constants to bias the hydrogen bond characterized by these AIMD potentials was nontrivial and required improvements to the EDS learning algorithm to improve the rate of convergence of EDS parameters and collective variables (CVs). 33 EDS learns the coupling constants during the progression of a simulation where the coupling constants at time segment + 1 are updated similar to a stochastic gradient descent algorithm 8 where 5 is the learning rate, 5 is the step size, and Δ ! " = 〈 ! ( ) " 〉 − ! , , and bold font represent vectors. Each time segment is set for a predefined period of simulation timesteps by the user, and should be about double the autocorrelation time of its corresponding collective variable.
Our current implementation of EDS uses the learning rate of White and Voth, and is calculated as where Ak is the maximum value a coupling constant can change. The simplest step size uses a derivative term that is proportional to the covariance of the collective variable error at time . With In the other case, the Levenberg-Marquardt (LM) algorithm was used for the step-size, where the derivative term in Eq. 4 is replaced with: where diag( ̿ 7 ̿ ) is a purely diagonal matrix with elements of ̿ 7 ̿ , and λ 6 is known as a mixing parameter which can tune the style of the step size. In this work, the LM algorithm is used to determine the coupling constants for EDS-BLYP and EDS-BLYP-D3 AIMD simulations. with a timestep of 0.5 fs. Coupling constants were determined in the constant NVT ensemble using three Nosé-Hoover chains at 298 K and with a time constant of 3000 cm -1 . EDS simulations were run until the CVs were determined to be converged, see Table 1

Results and Discussion
By examining Figure 1a, we can see that biasing the hydrogen bond with EDS can greatly improve the O-H RDF of BLYP. In fact, we see the first peak has nearly perfect agreement with the MBpol reference, and the structural properties past the first peak come into closer agreement with MBpol. Additionally, by biasing the hydrogen bond in BLYP water, we see greater improvements in the O-O RDF (Fig. 1b) and the H-H RDF (Fig. 1c). As shown in the O-H RDF, we see that the first peak in these other RDFs come into nearly perfect agreement with MB-pol, even though these other properties were not biased but they are still improved due to their strong correlations with the hydrogen bond. We note that the hydrogen bond EDS bias is targeting the close-range interactions between O-H pairs within a hydrogen bond. It should also be noted, however, that EDS greatly improves the structural properties past the first peak in all three RDFs (Fig. 1a-c), although not perfectly.
An improvement is also seen in the structural properties of dispersion corrected BLYP (BLYP-D3) over uncorrected BLYP when we bias the hydrogen bond with EDS. Similar to biasing the hydrogen bond in BLYP water simulations, we see that biasing the hydrogen bond in BLYP-D3 water simulations leads to RDFs whose first peaks nearly perfectly reproduces the first peaks in MB-pol (Fig 1d-f). Again, we emphasize that only the hydrogen bond is biased in these EDS-  Å, which is the radial distance where the coordination number is ~ 4. We see that in both the BLYP (Fig. 2a) and BLYP-D3 (Fig. 2b) cases, the unbiased simulations have a large probability above 100°, which is typically attributed to a strong tetrahedral ordering. By using EDS to bias the simulations, we see a decrease in both the biased BLYP and BLYP-D3 plots that is very close to MB-pol, with the most prominent improvement arising from biasing the hydrogen bond (OH) instead of the solvation structure (OO). Also, the peak around 55° in Figure 2a is captured with    38 . EDS-BLYP(OO) diffusion coefficients were taken from Ref. 34 . 9 is the standard deviation for diffusion constant. BLYP time constants show a monoexponential decay even when fitting with a bi-exponential. + and . . are time constants obtained via bi-exponential fits of the hydrogen bond ACD (Fig. 3).  Another goal in biasing the hydrogen bond in BLYP and BLYP-D3 AIMD water simulations was to see an improvement in the dynamics. The self-diffusion coefficient of BLYP water has been reported to be an order (or orders) of magnitude slower 9, 22-23 than that found in experiment (0.23 Å 2 /ps). 45 In the AIMD simulations reported here, the BLYP water oxygen self-diffusion coefficient was found to be 0.01 ± 0.01Å 2 /ps, and with an improvement to 0.06 ± 0.05 Å 2 /ps for the dispersion into even better agreement with experiment.) One can also confirm that the increased diffusion is not a result of a spurious increase in simulation temperature (see Table 2).  Table 2).  green line shows the distribution from using EDS to bias the solvation structure, and red line shows the distribution function from using EDS to bias the hydrogen bond. The hydrogen bond ACF comes from averaging over the first 40 ps of the simulation. The hydrogen bond ACF for EDS-BLYP-D3(OO) was not included due to the more limited statistical sampling in that prior work. Time constants from a biexponential fit of these curves can be found in Table 2. In Figure 4a we show the EDS potential (solid) and force (dashed) for EDS-BLYP(OH) and EDS-BLYP-D3(OH) simulations. The EDS bias in both cases has a repulsive "ramp" at short distances, which is followed by a potential energy well at larger distances. This indicates that using EDS to bias the hydrogen bond is not as simple as adding a repulsive bias between the pairs as was found when biasing the O-O solvation shell. 34 The EDS-BLYP potential has a sharper repulsive wall and deeper attractive well than the dispersion corrected EDS-BLYP-D3 case, which results in stronger repulsive and attractive force magnitudes. We also separated the EDS potentials into the 0 th and 2 nd moment contributions in Figure 4b. It is seen in both moments that the EDS-BLYP potential has larger magnitudes than EDS-BLYP-D3, which adds additional explanation for the differences in the EDS-BLYP and EDS-BLYP-D3 potentials in Fig. 4a.  In this work, the Experiment Directed Simulation (EDS) method has been implemented to bias the O-H hydrogen bond in BLYP and BLYP-D3 AIMD water to reproduce a target O-H RDF of the highly accurate MB-pol model (taken as the "exact" result). The rationale for this EDS bias is that most GGA-level DFT functionals result in hydrogen bonds that are too strong and hence AIMD water at 298 K that is over structured and much too slowly diffusing. By biasing only the hydrogen bond correlations, we were able to improve other structural properties and provide a very significant improvement in the water self-diffusion coefficient. Biasing the hydrogen bond opposed to the O-O solvation structure directly targets the over-polarization and charge transfer present in GGA DFT to obtain a better AIMD water model, and produces superior dynamical properties as shown in the prior text and Supporting Information. We note that once the EDS bias is determined, this results in the addition of only a classical correction term and hence negligible change to the computational cost of the AIMD.

Conclusions
An overall advantage of AIMD is it can also account for chemical reactions. This becomes particularly useful for simulating the hydrated excess proton or hydroxide, which diffuses through the Grotthuss mechanism. [46][47][48] In previous work, 34 we also implemented an EDS-AIMD O-O solvation structure bias to simulate AIMD water having an excess proton added to it, making the approximation that such a bias is still appropriate for the water solvent even once an excess proton is introduced into the system. It was found that such an approximation is quite accurate (and it is expected to be even more accurate for any system in which a solute or solutes are chemically more distinct from the water solvent). On the other hand, the same approximation cannot be simply implemented in the O-H EDS bias approach taken in this paper. This is because the EDS-AIMD hydrogen bond bias applies forces to all pairs of hydrogens and oxygens within the specified cutoff distances, but one would not want to apply this bias to any special hydrogens considered as part of a hydronium or hydroxide ion complex. Thus, in future work we plan to rigorously combine the EDS-AIMD hydrogen bond bias into a continuous potential for the hydrated excess proton (by combining EDS with a bond order analysis) in such a way that forces are applied to only select classes of hydrogens and oxygens. Such a result can also likely be generalized to other forms of chemical reactions in EDS-AIMD(OH) water in which the water participates in the chemistry. This research is currently in progress.

Supporting Information
Additional

S2
Section S1. Discussion on the Choice of Biased Moments We used EDS to bias the 0 th and 2 nd moment of the O-H RDF in AIMD simulations. In this discussion, we show that it is not necessary to bias the 1 st moment in these simulations to reproduce the O-H RDF of MB-pol In Figure S1 we show the 0 th (biased), 1 st (unbiased), and 2 nd (biased) moments in the constant NVE ensemble of the various simulation methods used in this study. We first note that both the 0 th moment and 2 nd moment in the EDS simulations are much closer to the target value than the unbiased simulations Additionally, we find the 1 st moment is much closer to the target value without a direct bias due to the strong correlations between the moments. This seems to suggest that directly biasing the 1 st moment is not necessary to improve its distribution.
We also found that including the 1 st moment into the EDS learning algorithm had a negligible effect on the resulting EDS potential energy. In Figure S2, we show the EDS-BLYP (OH) potential energy that results from biasing the 0 th and 2 nd moment, and the resulting EDS-BLYP (OH) potential after a 30 ps learning simulation where the 0 th , 1 st , and 2 nd moments were biased. We see excellent agreement between the two potential energy functions except at short distances where

Section S2. Calculation of the EDS Potential and Force
The EDS potential for a given pairwise distribution is defined in the main text as the product of the coupling constant ( ! / ! ) and the H-O pairwise function. (See Eq. 1). In practice, the PLUMED2 EDS module calculates the mean pairwise distribution over all O-H pairs and subtracts the target value from this mean value. The EDS bias is then the product of this difference and the coupling constant. Although this changes the potential, the forces remain unchanged due to the target values being a constant value. In the main text (Fig 4) and the supporting information (Fig.  S4), we chose to present the EDS potential between a single O-H pair using Eq. S1, where we have not subtracted the target value from the pairwise distribution, and have removed any dependence on the number of hydrogen bonds in the system. We show in Fig. S3 the corresponding EDS potential and force plots where the target value is subtracted from the pairwise function. EDS potential plots are the average over all three EDS simulations and the EDS force plots are calculated using finite difference of the potential. In this section, we present other important quantities that improved due to the EDS H-bond (O-H) biasing procedure. These are the: potential of mean force (PMF) (S4.1), and hydrogen bond 2D PMF (S4.2). We note that in all cases, biasing the hydrogen bonds using EDS was superior to biasing the oxygens (O-O). In parenthesis of each label signifies where the bias was on the solvation structure (OO) or the hydrogen bond (OH).  Table I and Equation S1.

S6
Section S3.2 Hydrogen-Bond 2D PMF. The free-energy is calculated via Eq. S3, where C is a constant chosen to set the minimum of the free-energy to 0 kcal/mol. By using EDS to bias the hydrogen bond, we see a broadening in the hydrogen-bond 2D PMF compared to the unbiased simulations.