Machine learning potentials from transfer learning of periodic correlated electronic structure methods: Application to liquid water with AFQMC, CCSD, and CCSD(T)

Obtaining the atomistic structure and dynamics of disordered condensed phase systems from first principles remains one of the forefront challenges of chemical theory. Here we exploit recent advances in periodic electronic structure to show that, by leveraging transfer learning starting from lower tier electronic structure methods, one can obtain machine learned potential energy surfaces for liquid water from the higher tier AFQMC, CCSD, and CCSD(T) approaches using $\le$200 energies. By performing both classical and path integral molecular dynamics simulations on these machine learned potential energy surfaces we uncover the interplay of dynamical electron correlation and nuclear quantum effects across the entire liquid range of water while providing a general strategy for efficiently utilizing periodic correlated electronic structure methods to explore disordered condensed phase systems.


I. INTRODUCTION
Ab initio molecular dynamics (AIMD) simulations, where the forces and energies are generated at each time-step by performing an electronic structure calculation, provide an appealing route to simulate reactive chemical dynamics. However, for disordered condensed phase systems an accurate description typically requires many 100's of atoms to obtain a chemically reasonable description of bulk systems (e.g., water) and this grows into the 1000's for more heterogeneous systems (e.g., those with interfaces). Since AIMD simulations require an electronic structure calculation to be performed at each time-step, to statistically converge even simple thermodynamic properties necessitates many tens of thousands of ab initio calculations (10's of picosecond timescale at a ∼1 fs time step) and for slower converging properties many millions are needed (nanosecond timescale). The computational expense of these simulations is further compounded if one wants to incorporate nuclear quantum effects (NQE) via ab initio path integral molecular dynamics simulations (PIMD). Due to its reasonable compromise between accuracy and efficiency, density functional theory (DFT) is currently the most frequently employed electronic structure method in condensed phase AIMD studies. However, the results depend-sometimes sensitively-on the choice of the exchange-correlation functional and the inclusion of dispersion corrections. 1,2 This issue motivates the use of beyond-DFT electronic structure theories, such as those based on the many-electron wavefunction. For example, work performed almost a decade ago demonstrated the AIMD simulations of liquid water using second-order Möller-Plesset perturbation a) Electronic mail: t.berkelbach@columbia.edu b) Electronic mail: drr2103@columbia.edu c) Electronic mail: tmarkland@stanford.edu theory (MP2). 3 However, the high cost of more accurate methods precludes their direct use in AIMD.
Machine learned potentials (MLPs) have emerged as an extremely promising approach to accurately model ab initio potential energy surfaces of condensed phase systems while being orders of magnitude more computationally efficient to evaluate. For liquid water, MLPs have been successfully developed at various levels of electronic structure ranging from different levels of DFT 4-10 to more recently using the random phase approximation (RPA) 11 and MP2. 12,13 The modeling of liquid water and other molecular systems with more accurate electronic structure methods, such as coupled-cluster theory or quantum Monte Carlo, has been so far limited to training on finite clusters of molecules. [14][15][16][17][18][19][20][21][22][23][24][25] Recent advances have opened the door to efficiently calculating the properties of periodic condensed phase systems using higher-level methods like coupled cluster singles and doubles without and with perturbative triples (CCSD and CCSD(T)) 26,27 and phaseless auxiliary-field quantum Monte Carlo (AFQMC). 28 However, while these advances allow the energies of many 100's of condensed phase configurations to be evaluated, this is still considerably less than what would be typically required to train an accurate MLP.
Here, we demonstrate that by using an approach based on transfer learning starting from a variety of lower tier electronic structure methods one can generate a highly data-efficient approach to training MLPs using high level electronic structure methods. With this approach, we show that only 200 highquality energies obtained from small periodic boxes containing 16 water molecules provide sufficient data for training our MLPs. Specifically, we train MLPs to periodic electronic structure calculations performed with AFQMC, CCSD, and CCSD(T). The MLPs are then used to perform AIMD and PIMD simulations of larger water boxes for the long times necessary to statistically converge static and dynamic properties using both classical and quantum mechanical treatment of nuclei. This allows us to achieve a careful comparison of the quality of the underlying electronic structure theories in describing water across its entire liquid temperature range and uncover the changes in the properties of liquid water as the dynamical electron correlation captured is increased. In addition, we provide a set of MLPs and a curated training set that can form the basis of future studies of water and aqueous systems. Our data-efficient approach provides a route to accurately obtain the properties of other disordered condensedphase systems by combining high-level electronic structure theory and machine learning.

A. Machine learning
To develop a highly data efficient strategy that only requires a minimal number of energies from periodic high level electronic structure methods, such as CCSD, CCSD(T), and AFQMC, and produce accurate MLPs, we exploited a combination of approaches including an active learning procedure for curating a small but comprehensive training set, the training of MLPs to energies for small periodic boxes (16 waters) and showing that they reproduce properties when used for larger simulations (64 waters) and a transfer learning approach that leverages the transferability of physics from lower tier electronic structure methods.

Developing an efficient training set by leveraging a committee of machine learning potentials
We first applied an active learning procedure to curate a data efficient set of configurations to train our MLPs. We invoked the commonly employed Query-by-Committee (QbC) 8,29,30 approach to iteratively add configurations to the training set. At every iteration, the current dataset of configuration energies was used to train a committee of 8 MLPs, each a Behler-Parrinello neural network potential 4,31 with a different random initialization of weights and random trainvalidation (90-10) splits of the dataset. To mitigate overfitting, we applied early stopping to each of these MLPs in the committee, monitoring the energy prediction error over a validation set. The mean potential energy surface obtained for this committee MLP was used to run a short MD simulation (SI Sec. III) that was terminated either when the system becomes unphysical (i.e., reaching a temperature greater than 400K) or when 50 ps of simulation trajectory was generated. We then selected the 10 configurations where the committee MLP had the largest standard deviation in its potential energy prediction to recalculate at the target level of electronic structure theory and add to the training set for the next iteration of this procedure. To prevent selected configurations from being overly correlated with one another, no two selected configurations were closer than 100 fs apart. The initial dataset of 50 configurations used to initialize this procedure was also selected via an iterative QbC procedure, as applied to a single 100 ps SCC-DFTB 32 (SI Sec. III) trajectory, that started with 10 randomly selected configurations and added 10 additional configurations for each of 4 subsequent iterations.
We generated 5 different 200 configuration datasets by running 5 instances of our QbC active learning scheme where the target level of electronic structure theory was DFT using the revPBE0-D3 33-35 functional (SI Sec. I A) which was chosen due to being computationally efficient compared to the high level electronic structure methods and since it has previously been shown to produce the properties of liquid water accurately when combined with path integral simulations. 36 To select the generated dataset on which the higher tier electronic structure methods would be used we then used each of the 5 datasets to train a committee MLP using revPBE0-D3. We then selected the one that gave the lowest force prediction error (RMSE) when evaluated on a test set of 1000 water configurations (64 water molecules) drawn from previously published AIMD simulations. 36 To check the transferability of the selected 200 configuration dataset, we recalculated the energies and forces for the generated datasets and the test set using the BLYP 37,38 functional (SI Sec. I A). After training and evaluating a new set of BLYP trained MLPs, the same dataset that resulted in the lowest error revPBE0-D3 trained MLP also gave the lowest error BLYP trained MLP (SI Fig. 3). Given the transferability of the relative utility of this dataset, we used this same 200 configuration dataset to train our CCSD, CCSD(T), and AFQMC MLPs.

Training machine learned potentials on configurations of small systems
Given the significant scaling of CCSD, CCSD(T), and AFQMC's computational cost with system size, we sought to reduce the system size of the configurations used in our training set. We investigated the feasibility of training our committee MLP to the energies of a set configurations for small periodic boxes to accurately predict for properties associated with a larger simulation box. For liquid water running periodic molecular dynamics simulations of small water boxes (32 water molecules or fewer) leads to significant artifacts even in simple properties such as radial distribution functions (RDFs) when compared to larger system sizes. 39 However, as demonstrated in SI Sec. VI for the revPBE0-D3 functional, if one trains a MLP on energies of small (16 molecule) periodic water configurations and then uses the resulting model to perform dynamics of a larger system (64 molecule), then the results obtained are in excellent agreement with those obtained from performing a AIMD simulation at the larger system size. Based on this demonstration we therefore trained the MLPs for the higher level methods on periodic configurations of 16 water molecules using the transfer learning approach described in the following section and then report the properties obtained in Sec. III by performing MD and PIMD simulations using 64 water molecules.
FIG. 1. Oxygen-oxygen RDFs for liquid water when running classical and PIMD simulations using NNPs fitted to either CCSD, CCSD(T), and AFQMC energies at both 300K and 370K. The PIMD oxygen-oxygen RDFs when using the CCSD(T) and AFQMC correspond closely with the experimental results shaded in grey at both 295K and 366K.

Transfer learning approach to train the MLPs
With only hundreds of energies from the higher tier periodic electronic structure methods available, to make efficient use of this data we employed a transfer learning approach. 40 To achieve this we first trained a committee MLP at a lower level of electronic structure theory, DFT or Hartree Fock (HF) on 531 configurations using both energies and gradients to improve the fitting. The parameters obtained for these MLPs were then used as the starting point for the fits to the higher level methods. This strategy exploits the idea that while the lower level methods may not reach the levels of chemical accuracy required for some applications they do contain fundamental physical information (e.g., about the fact that O and H when in close proximity form a high frequency chemical bond) that can be used to structure the neural network that underlies the MLP. Hence while the high level data is used to tune the accuracy of the MLP, it is leveraging the physics learnt by the initial training to the lower level method.
While using transfer learning to make efficient use of very small amounts of high level electronic structure data has distinct advantages over starting from a randomly initialized MLP, one must be mindful of the risk of hysteresis in the final MLP, i.e., that by biasing the weights by taking them from a model trained on a low level electronic structure method, that the final MLP will incorrectly contain remnants of the failures of the low level method. Hence to assess our transfer learning approach's ability to produce a final model that accurately reproduces the properties of the target high level electronic structure method, we initialized the procedure starting from a range of different low level methods that each give very different properties of liquid water. By comparing the structural and dynamical properties of liquid water obtained by performing MD and PIMD simulations on the final models obtained from these different starting points, one can thus evaluate which properties are obtained universally across the models initialized from different low level methods and thus accurately reflect the high level electronic structure approach.
In practice, following the convergence of the lower level fits to both energy and gradient data (see SI Sec. II), we retrained the committee MLP (8 separate MLPs) to the energies of the higher level method employing the extended Kalman filter optimizer as implemented in the n2p2 41 package and using a 90-10 train-validation split in order to monitor the prediction error over the validation set for early-stopping each individual fit to prevent overfitting. Before applying this transfer learning approach to the high level CCSD, CCSD(T) and AFQMC methods for which AIMD is not possible for the timescales required, in SI Sec. VI we tested it by training an MLP with revPBE0-D3 as the higher level method and BLYP and HF as the lower level methods from which the transfer learning was performed. We chose DFT with the revPBE0-D3 exchangecorrelation functional as the higher level in this benchmark due to its accurate description of the properties of liquid water and since we can compare the results of the transfer learned MLPs to AIMD and AI-PIMD trajectories that we have previously obtained. 36 We chose BLYP and HF as the lower level methods since the former gives an incredibly overstructured and dynamically sluggish description of water while the latter gives the opposite.
As shown in SI Sec. VI, MD simulations using the final transfer learned MLP models of revPBE0-D3 using only 200 energies starting from either BLYP or HF correctly capture the target RDFs and VDOS for liquid water at 300K obtained from AIMD simulations, which are both markedly different from those given by low level methods themselves (BLYP water has a more structured oxygen-oxygen RDF and higher wavenumber hydrogen VDOS bend and stretch peak positions than revPBE0-D3 water, and vice-versa for HF water). The agreement with the reference revPBE0-D3 results for both models is particularly strong for the RDFs. For the VDOS, the BLYP-initialized model outperforms the HF-initialized model in capturing the high frequency O-H stretch peak. Both models also accurately reproduce the RDFs obtained from AI-PIMD simulations of revPBE0-D3 for liquid water 36 at 300K (SI Fig. 5) with the only discrepancy again being in the VDOS (SI Fig. 10) where the BLYP-initialized MLP captures the full spectrum whereas HF shows an overstructured and blue shifted OH stretch region.
Additionally, we compared our transfer learning approach to two common alternative machine learning approaches: directly training a committee MLP on the high level energies starting from randomly initialized weights and training a committee delta-learning model that corrects from the lower to higher level method. To allow for a fair comparison, all three types of models were trained to the same 200 configuration training set of energies (revPBE0-D3) for liquid water using the same MLP architecture and optimization settings. The randomly initialized model resulted in a potential that was unstable for the purposes of running MD for a periodic simulation box of 64 water molecules, with the instantaneous temperature drifting severely after the first simulation step. For the delta-learning model, we trained a committee MLP to capture the energy difference between revPBE0-D3 and another committee MLP trained to BLYP. The delta-learning model was similarly as unstable as the randomly initialized model. Hence for liquid water with 200 energies at the higher level we found that the transfer learning approach we have detailed above provides the most accurate results.
Given the demonstrated efficacy of our transfer learning procedure, in Sec. III we applied this approach to train committee MLPs to CCSD, CCSD(T), and AFQMC energies for liquid water using three different sets of MLP initializations: HF, BLYP, and revPBE0-D3.

B. Correlated electronic structure methods
We perform correlated electronic structure calculations of liquid water with periodic boundary conditions at the Γpoint using AFQMC, CCSD, and CCSD(T). Periodic CCSD and CCSD(T) calculations were performed using PySCF 42,43 where electron-repulsion integrals are handled using the range-separated density fitting method, 44,45 and AFQMC calculations were performed using QMCPACK 46,47 and ipie. 48 These calculations all began with a periodic spin-restricted HF calculation also perfomed using PySCF. We provide brief details here, and further information can be found in SI Secs. I B and I C.
AFQMC is a projector Monte Carlo method where the ground state of a given Hamiltonian is obtained via imaginary time evolution. Without any approximations, it scales exponentially as the system size grows due to the fermionic sign problem. We use the phaseless approximation 28 to obtain an algorithm that scales with the system size N as for each sample at the expense of introducing errors in the final ground state energy estimate. The phaseless approximation sets a boundary condition in imaginary time evolution using a priori chosen wavefunction called the trial wavefunction. The bias from this approximation becomes smaller as one improves the quality of trial wavefunctions. In this work we employ the simplest trial wavefunction, the spin-restricted Hartree-Fock determinant. A recent benchmark study examined the accuracy of AFQMC with Hartree-Fock trial wavefunctions over 1004 data points, 49 and based on these results at this level of approximation we expect the accuracy of AFQMC to be between CCSD and CCSD(T) for the problems considered here.
Coupled cluster (CC) parameterizes the electronic wavefunction using an exponential ansatz, |Ψ CC = eT |Φ HF , wherê T creates all possible particle-hole excitations from the HF reference and is determined by iteratively solving a set of coupled non-linear equations. 50 CCSD approximates the full CC ansatz by truncating the T -operator at single and doubleexcitation levels,T =T 1 +T 2 . 26 CCSD(T) improves upon CCSD by further including the contribution from triple excitations in a perturbative (non-iterative) manner 27 and is often referred to as the "gold standard" in the quantum chemistry of simple molecules. The computational cost of CCSD and CCSD(T) scales more steeply than that of AFQMC, as O(N 6 ) and O(N 7 ), respectively. In this work, we avoid the high computational cost of full CCSD and full CCSD(T) by systematically compressing the virtual space using the frozen natural orbital (FNO) approximation, 51,52 which we confirmed to introduce a negligible error (SI Sec. I B).

III. RESULTS AND DISCUSSION
Having established the applicability of our data efficient approach to training MLPs that accurately reproduce the target potential energy surfaces on DFT, we now apply it to obtain the static and dynamic properties of liquid water using three correlated methods: CCSD, CCSD(T), and AFQMC. We trained a committee MLP, which as described in Sec. II A 1 is formed of the mean of 8 independently trained MLPs, for each correlated electronic structure method. By performing MD and PIMD simulations in the NVT ensemble at 300 K and 370 K using the committee MLP for each electronic structure method, here we evaluate how different treatments of dynamical electron correlation affect the selected properties when the nuclei are treated classically or quantum mechanically.
To provide an assessment of the accuracy of our transfer learned committee MLPs on different structural and dynamical properties as discussed in Sec. II A 3 we can compare the consistency of the results obtained from models that have been initially trained to different lower level electronic structure methods. Hence, for each of the correlated methods we trained three transfer learned committee MLPs starting from HF, BLYP, and revPBE0-D3 as our lower level methods which are known to understructure, overstructure, and accurately reproduce the properties of water respectively. Performing classical MD at 300K from MLPs starting from these three different methods SI Fig. 13 shows that regardless of the lower level that the transfer learning was initialized with, the properties obtained at the high level, AFQMC in this case, coincide closely although the O-H stretch for the HF-initialized model is blue-shifted slightly. This demonstrates that our training set of 200 energies of periodic boxes consisting of 16 water molecules is large and diverse enough to train an accurate model of liquid water under these conditions using our transfer learning protocol. When nuclear quantum effects are included by peforming PIMD simulations using the committee MLPs, the oxygen-oxygen RDF and VDOS for the revPBE0- D3 and BLYP initialized transfer learning models remain consistent with one another (SI Fig. 14).
In this section we show the results from the committee MLP generated using transfer learning from revPBE0-D3 for each of the correlated methods since, as shown in SI Fig. 2, revPBE0-D3 shows the strongest correlation with the training set energies of CCSD, CCSD(T), and AFQMC out of the three initialization methods we used.
For our AFQMC results, it is important to note that the AFQMC energies contain stochastic error, with each of our N=200 training set AFQMC energies having a corresponding estimate for the standard error that ranges from 1-2 mH depending on the specific training set configuration. To evaluate how this level of error might affect our reported properties, we employed a test where we sample new training sets where the same N=200 configuration are used but a random value is added to each AFQMC energy to reflect the uncertainty of our AFQMC energies (see SI Sec. VII for details). In total, 12 training sets were sampled and a transfer learned committee MLP was fit to each. SI Figs. 11 and 12 serve to quantify the variations in the oxygen-oxygen RDF and hydrogen VDOS obtained from the 12 separate training sets due to the stochastic error in the AFQMC energies, with the grey shading representing the standard deviations. From this test we found that the stochastic error in our AFQMC calculations does introduce noticeable variations in the OH stretch peak of the hydrogen VDOS, particularly around the top of the peak, but the RDFs and VDOS are otherwise consistent for the different training sets.

A. Static properties of water from correlated electronic structure methods
We first compare the static equilibrium properties for liquid water at 300 K and 370 K obtained via classical MD and PIMD simulations using the committee MLP for each correlated electronic structure method. For these properties PIMD exactly includes the NQEs for distinguishable particles, which is a highly reliable assumption for nuclei at this temperature. Figure 1 shows the oxygen-oxygen RDFs for each of CCSD, CCSD(T), and AFQMC as compared to the experi-mental results at 295 K 53 and 366 K. 54 At 300 K, classical MD CCSD(T) and AFQMC both give a first peak that is slightly higher than observed experimentally suggesting the liquid is overstructured. However, once NQEs are accounted for in the PIMD simulations both RDFs become slightly less structured and show better agreement with experiment, with that of CCSD(T) coinciding quantitatively. CCSD gives good agreement with experiment when used in classical MD simulations but is understructured when NQEs are included which is consistent with this electronic structure approach starting from a HF reference which gives a severely understructured liquid with the additional dynamical correlation added through the tiers of CC theory progressively structuring the liquid. At 370 K when used in PIMD simulations all methods give good agreement with the experiment with AFQMC again exhibiting a first peak that is higher than CCSD, CCSD(T) and experiment. SI Fig. 15 and Fig. 16 show the hydrogen-hydrogen and oxygen-hydrogen RDFs at 300 K sampled via classical MD and PIMD, respectively, where all three electronic structure methods give similar results but AFQMC again exhibits a more structured hydrogen bond network with the first intermolecular OH peak at ∼1.85 Å, which corresponds to hydrogen bonds, being slightly higher than the other methods.
The tetrahedral order parameter provides a measure of higher order structural correlations within water's hydrogen bond network beyond the purely radial information encoded in the RDFs. The tetrahedral order parameter q is defined for a given water molecule as, 55 where θ jk is the angle that a given water molecule's oxygen atom makes with two neighboring oxygen atoms j and k.
The tetrahedral order parameter thus ranges from 0 to 1 with higher values indicating that the hydrogen bond network possesses angles closer to that of a perfect tetrahedral arrangement of the four nearest neighbour oxygen atoms around a central water molecule. As shown in Fig. 2 at 300 K for both classical MD and PIMD this property further highlights the understructured hydrogen bond network of CCSD compared to the more accurate correlated methods, CCSD(T) and AFQMC, that are in close agreement. At 370 K the distribution of the tetrahedral order parameter for all three methods shifts to lower values. The comparison of these static equilibrium properties suggests that the inclusion of higher order electron correlation contributions in methods like CCSD(T) and AFQMC, as compared to CCSD or HF, results in a greater degree of structuring in liquid water at 300 K and 370 K. Given the directional nature of this additional structuring, as seen from the greater probability density at higher q in Figure 2, this arises from slightly stronger hydrogen bonds being formed when using the two higher level methods. Our comparisons of the oxygenoxygen RDFs obtained from classical and PIMD simulations at 300 K also show that for these correlated methods the inclusion of NQEs works to slightly destructure liquid water. The relatively subtle overall effect of NQEs on liquid water around 300 K is known to arise from the close balance of competing quantum effects in this system. 56

B. Dynamical properties of water from correlated electronic structure methods
We now turn our attention to how the different correlated electronic structure methods behave when used to compute dynamical properties of liquid water, namely the self diffusion coefficient and VDOS. For these properties we compare the results obtained from classical MD and TRPMD. For these properties, since real time quantum dynamics is intractable for such a large atomistic condensed phase system for the timescales required to compute these properties, we use TRPMD to approximate the role of NQEs. TRPMD has previously been shown to be an accurate way for treating NQEs in the dynamics of condensed phase systems, however it is known to spuriously broaden high-frequency vibrational modes. 36,59 The diffusion coefficients obtained for the three correlated methods are shown in Table I at 300 K and 370 K. Unlike the other properties we report, diffusion coefficients exhibit a notable scaling with system size that must be corrected for to make comparisons with experiment. Hence, as described in SI Sec. IV A, the diffusion coefficients were computed using simulations of 64 water molecules and then extrapolated to the infinite system size limit using the previously derived system size scaling relation 60 and the experimental viscosity of water. 61 At 300 K, where the experimentally measured diffusion coefficient is 2.41±0.05×10 −9 m 2 /s, 57 when classical MD is used AFQMC and CCSD(T) yield smaller diffusion coefficients than observed experimentally, 2.09±0.05 and 2.21±0.06 ×10 −9 m 2 /s respectively, while CCSD gives a larger value 2.60±0.08 ×10 −9 m 2 /s. This behavior is consistent with the trends observed for the electronic structure approaches in the structural properties, where CCSD formed a slightly understructured hydrogen bond network compared to the more accurate correlated methods, which here leads to faster dynamics. Upon including NQEs using TRPMD, the diffusion coefficients for all three electronic methods increase, consistent with the disruption of the hydrogen bond network upon including zero-point energy, which brings the CCSD(T) result (2.30±0.08 ×10 −9 m 2 /s) to within the statistical error bar of the experimentally observed value. Even with NQEs included the AFQMC diffusion coefficient is lower (2.16±0.09 ×10 −9 m 2 /s) than experiment, consistent with it forming a more structured liquid. However, it should be noted that the discrepancy in the diffusion coefficient is very small; to change water's diffusion coefficient from that observed at 300 K via classical MD using our CCSD(T) model to the value obtained by performing TRPMD dynamics would require less than a 2 K change in the temperature of the liquid. 58 In addition, for dynamical properties TRPMD only approximately includes NQEs and hence the better agreement of CCSD(T) with the experimental value could be changed if a different approach was used to treat the quantum dynamics of the nuclei. At 370 K, when TRPMD is used CCSD and CCSD(T) give similar results (8.29±0.14 and 8.16±0.14 ×10 −9 m 2 /s, respectively), both of which are close to the experimental value of 8.26±0.02 ×10 −9 m 2 /s. 58 The diffusion of AFQMC water is again considerably slower, which is consistent with the relatively greater degree of structure we saw in the latter's respective oxygen-oxygen RDF and q distribution.
The VDOS in Figure 3 provides more information on the frequency dependence of the dynamics of water for the three electronic structure methods, since the diffusion coefficient is simply proportional to its zero frequency limit. All three methods give qualitatively similar VDOS around the lower frequency librational band (∼500 cm −1 ) and peak associated with the H-O-H bending mode (∼1600 cm −1 ). At 300 K using classical MD, the main qualitative difference lies in the OH stretch peak (3000-4000 cm −1 ), with CCSD(T) giving a peak that is redshifted by ∼80 cm −1 with respect to the CCSD peak, while the AFQMC peak is slightly broader than the others and is centered closer to the CCSD result. At the low-frequency end of the O-H stretch peak, the VDOS for CCSD(T) and AFQMC coincide with one another. Low frequency O-H stretches are typically associated with stronger hydrogen bonds and the frequency of the O-H stretch peak has previously been demonstrated to be inversely correlated with the tetrahedrality parameter. 5 Hence, the consistency between CCSD(T) and AFQMC at the low frequency part of their respective O-H stretch peaks is consistent with the structural evidence in Sec. III A showing that these two methods similarly structure water via slightly stronger hydrogen bonds as compared to CCSD. The TRPMD simulations at 300 K, which include NQEs via TRPMD simulations, lead to a ∼120 cm −1 red-shift and broadening of the stretch peak and a ∼100 cm −1 shift in the bend for all three methods, which is consistent with observations from previous studies. 36,59 Although some of the broadening likely arises from physical effects, TRPMD is known to introduce spurious broadening of high-frequency vibrational modes. 36,59 At 370 K the most significant difference is the shift in the zero frequency intensity, which we expect since this is directly related to the self-diffusion coefficient and otherwise the VDOS is qualitatively similar with respect to the 300 K results. This is expected since for the high frequency modes the zero point energy greatly exceeds the thermal energy in the the mode, i.e. k B T << ω/2 where k B is the Boltzmann constant, T is the temperature and ω is the frequency of the mode.
Overall, the trends we observe in the dynamical properties largely mirror the evidence we presented for the structural properties showing that CCSD results in an understructured description of liquid water at 300 K, as compared to CCSD(T) and AFQMC, while at 300 K and 370 K AFQMC overstruc-tures water. The differences in the diffusion coefficients in Table I reflect this trend, with CCSD overestimating the experimental diffusion coefficient at 300 K with its understructured description of water and AFQMC underpredicting the diffusion coefficient at 370 K. With the VDOS the main differences between the three correlated methods manifest in the O-H stretch peak positions and breadth, with peaks given by MD at 300 K for both CCSD(T) and AFQMC skewed more to lower frequencies that are associated with stronger hydrogen bonds.

IV. CONCLUSION
In summary, we leveraged developments in high-level periodic electronic structure theory and exploited methods to improve the data efficiency of fitting MLPs to investigate the static and dynamical properties of liquid water at the level of CCSD, CCSD(T), and AFQMC. We devised a data efficient protocol for training MLPs that uses small periodic boxes of water (16 molecules) sampled judiciously via an iterative QbC active learning procedure. To make the most out of the few configuration energies we can afford to compute, we also employed a transfer learning approach that leverages the transferability of physics between lower level electronic structure methods (e.g. DFT with the BLYP functional, revPBE0-D3 functional, or HF) and our target higher-level methods, using MLPs fit to the former to initialize a fine-tuning transfer learning fit to the latter. Using this approach we showed that we can train stable MLPs with as few as 50 configuration energies, capture the RDFs with 100 (SI Fig. 4), and with 200 configuration energies obtain both accurate static and dynamical properties such as the diffusion constant and VDOS (SI Fig. 9). In contrast, with these same 200 energies we were unable to train stable models using delta learning or using random initialization of the model. We used our MLPs trained to CCSD, CCSD(T), and AFQMC to examine how different static and dynamical properties of liquid water are affected by the level of dynamic electron correlation accounted for and the inclusion of NQEs. Our results indicate that CCSD tends to understructure liquid water and overpredict the diffusion coefficient, as compared to experiment. On the other hand, both CCSD(T) and AFQMC give oxygen-oxygen RDFs and diffusion coefficients that are more consistent with experimental values at 300 K, suggesting that the more accurate treatment of dynamical correlation present in these methods is sufficient for describing liquid water. The inclusion of NQEs for our 300 K simulations bring the CCSD(T) and AFQMC results in even closer agreement with experiment and seems to generally manifest as a slight destructuring of liquid water for all three electronic structure descriptions. This small destructuring upon including NQEs for these correlated methods is in contrast to some DFT exchange correlation functionals where due to the overprediction of the anharmonicity of the O-H coordinate the inclusion of NQEs works to structure the liquid phase. 36 Ultimately, we envision that the configurations and energies that form the training dataset, the resulting MLPs, and the pro-tocols we employed here will be useful in their own separate respects for future work in modeling potential energy surfaces for condensed phase systems.   [4][5][6] functionals. Our SCC-DFTB [7] calculations, used to perform an initial MD simulation for sampling the first 50 training set configurations, were also conducted using CP2K. All calculations were run using periodic boundary conditions.
We used Goedecker-Tetter-Hutter (GTH) pseudopotentials [8] for our HF and DFT calculations conducted in CP2K. For our BLYP calculations, we employed a mixed Gaussian and plane wave method [9,10] using the DZVP [11] basis set and an auxiliary plane wave basis with a 300 Ry cutoff. These HF and revPBE0 calculations made use of the TZV2P [11] basis set and a 400 Ry cutoff auxiliary plane wave basis. We employed a truncated Coulomb operator [12] with a cutoff of 3.9Å and the auxiliary density matrix method [13] with the cpFIT3 fitting basis set. The short cutoff we employed for the truncated Coulomb operator was due to the small periodic boxes we used for our training set (side lengths of 7.831Å).
In checking the effects of employing such a short cutoff, and in general these small periodic boxes for the training set, we found that the effects seemed minimal with regards to the properties we are examining in this work (see SI Sec. VI for more details). For the revPBE0 calculations, we applied D3 dispersion correction [14].
Our SCC-DFTB calculations were run using parametrizations provided by CP2K for O and H atoms. We also employed an Ewald sum [15], in the form of a smooth particle mesh using beta-Euler splines [16], for electrostatics.

B. CCSD and CCSD(T)
We performed our CCSD and CCSD(T) calculations using the PySCF [17] software package and the cc-pVQZ basis set. To reduce the high computational cost of these calculations, we use the frozen natural orbital (FNO) approximation [18] as detailed in previous work [19] to systematically compress the virtual space. Particularly, we found that using 200 virtual orbitals is enough to reach a maximum absolute error of 0.03 kcal/mol for CCSD and 0.02 kcal/mol for CCSD(T) or a root-mean-square error of 0.02 kcal/mol for CCSD and 0.01 kcal/mol for CCSD(T) in the relative energy (per water molecule) between 10 ran-domly drawn 16-water molecule configurations. We thus performed all CCSD and CCSD(T) calculations in this work with the FNO approximation and 200 virtual orbitals.
Note that the AFQMC energies we ultimately trained our MLPs on used the HF energies

II. MACHINE LEARNING DETAILS
All of our committee machine learning potentials (MLPs), each consisting of 8 separate MLPs, were fit using the n2p2 package [25] v2.1.0. We employed the same set of input features (symmetry functions) and neural network architectures for all individual MLPs in our committee approach as was used in previous work fitting MLPs for liquid water [26].
In short, for our MLPs we used Behler-Parinello high dimensional neural network potentials employing local atom-centered symmetry function descriptors [27] that decompose the potential energy of a system as a sum of predicted atomic contributions.
When fitting our MLPs to both configuration energies and atomic forces, as was the case with our HF and BLYP committee MLPs, we employed the same training and optimizer hyperparameters (e.g. a 90-10 training to validation set split, early-stopping, Nguyen-Widrow weight initializiation [28], the extended Kalman Filter [29][30][31], etc.) as were used and benchmarked in previously published works [26,32,33]. These force-fitted MLPs were fit to 531 configurations of 16 waters and were used to initialize refits targeting revPBE0-D3, which served as preliminary checks on the effectiveness of this transfer learning procedure. For the refits to CCSD, CCSD(T), and AFQMC, in addition to starting from HF and BLYP, we also used the generation 1 revPBE0-D3 fits [33] from Schran et. al. to initialize the transfer learning models. All of the refits we conducted were trained using just the energies of 200 configurations and employed early-stopping, using a 90-10 training to validation set split, to mitigate overfitting.

III. MOLECULAR DYNAMICS SIMULATION DETAILS
All classical MD and PIMD NVT simulations employing our committee of ML potentials were conducted using the i-PI [34,35] and LAMMPS [36] packages in conjunction, the latter of which is interfaced to n2p2 [32]. Our classical MD simulations were run using a 0. to compute reported properties were performed with boxes of 64 water molecules. For our 300 K simulations, we used the experimental density of liquid water at that temperature, 0.996 g/cm 3 [39,40], to obtain a simulation box length of 12.432Å, and at 370 K we used the experimental density of 0.960 g/cm 3 [40] to obtain a box length of 12.590Å.
The MD simulations we ran as part of our active learning scheme for building our training set were run for 50 ps and for a system consisting of 16 water molecules in a simulation box with side lengths of 7.831Å, with all other settings kept the same as the simulations used to compute properties. The initial 100 ps MD simulation using DFTB that we used to initialize our active sampling scheme was also run using the same settings.

A. Self-diffusion coefficients
We calculated the molecular self-diffusion coefficients by performing a linear fit of the long-time slope of the mean square displacement (MSD) of the oxygen atoms vs. time. We calculate our MSDs out to 10 ps and fit the region from 1 ps to 10 ps. The value we report is the mean diffusion coefficient as calculated over separate 40 ps chunks of the sampled trajectory and the error bars correspond to the standard deviation of the mean.
We used the following relation to correct for the systems size scaling of the simulated diffusion coefficient [41]: (1) L is the length of the simulation box, ξ is a geometry specific constant that is 2.837297 for a cubic box, and η is the shear viscosity (we use the experimental values at T=300 K and 370 K depending on which NVT ensemble we are sampling [42]).

B. Hydrogen vibrational density of states (VDOS)
The hydrogen VDOS was computed by Fourier transforming in time the velocity autocorrelation function. We calculated the velocity autocorrelation function out to 2 ps.
V. TRANSFERABILITY OF RELATIVE TRAINING SET PERFORMANCE SI Fig. 3 demonstrates that the relative performance of our training sets sampled via the active learning procedure we used is transferable between different electronic structure methods. Specifically, we see that training set 1 (blue), which contains configurations sampled with revPBE0-D3 as the target level of electronic structure theory, results in the lowest test set force error predictions for when revPBE0-D3 is the target (left) and also when BLYP is the target (right). Conversely, training set 4 (red) is the highest error training set for both revPBE0-D3 and BLYP and in general the relative performance of training sets is consistent for both DFT functionals. Given the transferability of the relative performance of these sets of training configurations across different electronic structure methods, we whose to use training set 1 for our CCSD, CCSD(T), and AFQMC transfer learning fits.

VI. TESTING TRANSFER LEARNING PROCEDURE ON DFT/REVPBE0-D3
The results in this section show that the transfer learning procedure we employed starting with a committee of MLPs trained to either HF or BLYP energies and forces that are then refitted using revPBE0-D3 energies works successfully in reproducing certain target AIMD and PIMD properties. These results also show that it is possible to train MLPs on small VDOS (shaded in green) [43] when trained on as few as N=200 revPBE0-D3 energies. The first row of SI training on N=200 configurations energies, starting with initial weights from fits to HF and DFT using the BLYP functional. The procedure, when using a training set of N=200 energy calculations, roughly reproduces the q distribution of the target method when using initial weights from a fit to DFT/BLYP. showing that the latter approximately reproduces the hydrogen VDOS of the target method.

VII. TEST EVALUATING EFFECT OF STOCHASTIC ERROR IN AFQMC EN-ERGIES
Here we show results from a test we conducted on how the level of stochastic error in the AFQMC energies affects the properties we are benchmarking. Each of our N=200 training set configurations have a calculated AFQMC energy and a corresponding estimate for the standard error that ranges from 1-2 mH depending on the specific configuration. For this test, we sampled new training sets where the same N=200 configuration are used but a random value is added to each AFQMC energy. These random values were drawn from Gaussian distributions with standard deviations set by the estimated standard error for a given configurations AFQMC energy.
We generated 12 separate training sets in this fashion and the panels on the left side of SI  Note that the grey shading around the mean denotes the standard deviations. Note that the grey shading around the mean denotes the standard deviations.

VIII. COMPARISON OF AFQMC TRANSFER LEARNING RESULTS USING HF, BLYP, AND REVPBE0-D3 MLPS TO INITIALIZE WEIGHTS
The results in the main text were obtained using committee MLPs that were initially trained to revPBE0-D3 and then refitted to energies from either AFQMC, CCSD, or CCSD(T). SI Figs. 13 and 14 show that the results obtained from different transfer learned models, whether initialized with weights trained to HF, BLYP, or revPBE0-D3 energies and forces, that were refitted to our N=200 AFQMC training set give similar oxygen-oxygen RDFs and hydrogen VDOS. The consistency in these results despite the varying initial conditions gives us additional confidence, on top of our benchmark results presented in SI Sec. VI, in the robustness of this transfer learning approach.