Breaking the Coupled Cluster Barrier for Machine-Learned Potentials of Large Molecules: The Case of 15-Atom Acetylacetone

Machine-learned potential energy surfaces (PESs) for molecules with more than 10 atoms are typically forced to use lower-level electronic structure methods such as density functional theory (DFT) and second-order Møller–Plesset perturbation theory (MP2). While these are efficient and realistic, they fall short of the accuracy of the “gold standard” coupled-cluster method, especially with respect to reaction and isomerization barriers. We report a major step forward in applying a Δ-machine learning method to the challenging case of acetylacetone, whose MP2 barrier height for H-atom transfer is low by roughly 1.1 kcal/mol relative to the benchmark CCSD(T) barrier of 3.2 kcal/mol. From a database of 2151 local CCSD(T) energies and training with as few as 430 energies, we obtain a new PES with a barrier of 3.5 kcal/mol in agreement with the LCCSD(T) barrier of 3.5 kcal/mol and close to the benchmark value. Tunneling splittings due to H-atom transfer are calculated using this new PES, providing improved estimates over previous ones obtained using an MP2-based PES.

larger systems. Our group has made use of this approach for PIP PESs of N-methylacetamide, 9,10 glycine, 11 and tropolone. 12 We also recently reported PIP and fragmented PIP PESs (see below for some details) for 15-atom acetylacetone (AcAc), 13 using MP2 energies and gradients from Meuwly and co-workers, 14 supplemented by us with roughly 500 additional configurations. This PES has a barrier for symmetric H-atom transfer of 2.13 kcal/mol (745 cm −1 ), in close agreement with the direct MP2 value of 2.18 kcal/mol (763 cm −1 ). However, that value of the barrier is more than 1 kcal/ mol below the reported CCSD(T)/aug-cc-pVTZ barrier of 3.2 kcal/mol. 15 It is expected that this error in the MP2-based PES leads to a large overestimate of the tunneling splitting for the ground vibrational state H-atom transfer. Nevertheless, the splitting was obtained with the MP2-based PES, using fulldimensional diffusion Monte Carlo calculations. The splitting is 160 cm −1 with an uncertainty of 15 cm −1 . Using a simple 1d model, a splitting of 113 cm −1 was obtained for a barrier of 2.2 kcal/mol and 74 cm −1 for a scaled barrier of 3.2 kcal/mol. This simple 1d estimate for the larger barrier is not expected to be quantitative; however, it does confirm that a large decrease in the splitting upon increasing the barrier by 1 kcal/mol can be expected.
This magnitude of the error in chemical barriers is typical for MP2 and DFT accuracy compared to benchmark CCSD(T) results. In general, MP2 and DFT geometries and harmonic frequencies are relatively more accurate than barrier heights. Thus, there is a strong motivation to improve a PES based on a lower-level method such as DFT and MP2 when the focus is on barriers.
Recent approaches to do this, using ML, aim to bring a PES based on a low-level of electronic structure theory to a higher level. There are two approaches currently being investigated to accomplish this goal. One is transfer learning (TL), which has been developed extensively in the context of artificial neural networks, 16 and much of the work in that field has been brought into chemistry. 14,17−20 The basic idea of TL is that a fit obtained from one source of data (perhaps a large one) can be corrected for a related problem by using limited data and by making small training alterations to the parameters obtained in the first fit. Therefore, in the present context of PES fitting, an ML-PES fit to low-level electronic energies/gradients can be reused as the starting point of the model for an ML-PES with the accuracy of a high-level electronic structure theory. As noted, this is typically done with artificial neural networks, where it is hoped that weights and biases trained on lower-level data require minor changes in response to additional training using high-level data.
The other approach is Δ-machine learning. In this approach, a correction is made to a property obtained using an efficient, low-level ab initio theory. 18−21 The focus of most work on TL or Δ-learning has been on developing transferable force fields, with applications mainly in the thermochemistry and molecular dynamics simulations.
Meuwly and co-workers applied TL using thousands of local CCSD(T) energies to improve their MP2-based neural network PESs for malonaldehyde, acetoacetaldehyde, and acetylacetone (AcAc). 14 We recently proposed and tested a Δ-learning approach that uses a small number of CCSD(T) energies to correct a PES based on DFT electronic energies and gradients. 22 The method was validated for PESs of small molecules, CH 4 and H 3 O + , and for 12-atom N-methylacetamide. In all cases, the coupled cluster energies were obtained over the same large span of configurations used to get the lower-level PES. For N-methylacetamide, these included the cis and trans isomers and the saddle points separating them.
Here we apply this Δ-learning approach to 15-atom AcAc, C 5 H 8 O 2 . The approach is to construct a high-level, coupledcluster-level PES starting from a lower level MP2 one by using a correction PES. This correction PES is a fit to a small number of high-level ab initio energies that span the same range of configurations used to obtain the lower-level PES. Explicitly, the corrected high-level PES, denoted V LL→CC , is given by where V LL is the lower-level PES and ΔV CC−LL is the correction PES. In the present application to AcAc, we calculated 2151 LCCSD(T)-F12/cc-pVTZ-F12 energies 23 and performed training on subsets of these ranging in size from 430 to 1935. By contrast, V LL was fit using a data size of 250 884 MP2 energies and gradients.
In the PIP approach to fitting, 24 the potential V is represented in compact notation by where c i are linear coefficients, p i are PIPs, n p is the total number of polynomials for a given maximum polynomial order, and x represents Morse variables. For example, x αβ is given by exp(− r αβ /λ), where r αβ is the internuclear distance between atoms α and β. The range (hyper) parameter, λ, was chosen to be 2 bohr. The linear coefficients are obtained using standard least-squares methods for large data sets of electronic energies (and for large molecular gradients as well) at scattered geometries.
For molecules of more than 10 atoms, the size of the PIP basis can become a computational bottleneck. This size depends in a complicated and nonlinear way on the maximum polynomial order, the number of Morse variables, and the order of the symmetric group. 24 While we have been able to use a full PIP basis even for a 15-atom tropolone 12 and AcAc, 13 we have shown that the fragmented PIP, which can be applied to larger molecules, performs very well and runs faster than the full PIP basis.
The fragmented PIP basis is obtained by fragmenting a molecule into groups of atoms. A PIP basis for each group can be calculated rapidly and then combined with those of other groups to provide a compact and still precise representation of the PES. 9 Indeed this has been verified for N-methylacetamide, 9,10 tropolone, 12 and AcAc. 13 Note that these PESs use the most recent software that includes gradients in the fit and produces gradients on output. 25−27 The histogram of the MP2 database from our AcAc PES 13 is shown in Figure 1(a); it includes 5454 geometries. The MP2 energies are relative to the MP2 global minimum. The two lowest bars of the histogram are primarily the 454 points added by our group via an AIMD trajectory or random grids starting from or centered on the global minimum (GM) and H-atom transfer transition-state saddle point geometries but are not specifically chosen to be along the reaction coordinate. Hereafter, unless indicated otherwise, the H-atom transfer saddle point is denoted as "SP". The AIMD trajectory was run at 400 cm −1 for 200 steps starting from the GM, and 50 geometries/energies were used. The grids were made from taking the Cartesian coordinates of the GM and SP and varying, for each geometry, all of the 45 coordinates by a random number between the coordinate ± del . For each of these starting points, approximately 75 geometries were generated with del = 0.001 Å, and approximately 150 geometries were generated with del = 0.004 Å. In all, our group calculated 454 points for the MP2 database. A majority of the geometries and energies (5000 points) were kindly provided to us by the Meuwly group on the basis of their AIMD trajectories.
The database we used to calculate the difference potential, ΔV CC−LL , is described by the energy histogram in Figure 1(b). The LCCSD(T) energies are relative to the LCCSD(T) energy at the MP2 global minimum structure. The database is composed of the 500 geometries that are clustered near the GM, as determined by choosing the RMS bond difference for the 105 bond lengths to be less than 0.3 Å when comparing the geometry to that of the GM, and another 1651 taken by (a) choosing a random integer that indicated their position on a list of the remaining MP2 points, (b) discarding choices whose MP2 energies were more than 30 000 cm −1 , and (c) stopping the selection when the requisite number of choices had been The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter made. Again, the aim was to produce a difference PES that spans the same range of configurations as the original lowerlevel PES, namely, one that includes points near the GM and SP but not specifically chosen to be along the reaction coordinate. The difference database shows a peak near a zero energy difference with a distribution ranging from approximately −1500 to 800 cm −1 . Figure 1(c) shows the correlation between the MP2 energies and those of the difference potential (LCCSD(T)-MP2). Note that the MP2 energies at which LCCSD(T) energies were calculated span the range from 0 to about 30 000 cm −1 . It is clear that there is no systematic correlation between the MP2 energy and the difference energy. The largest energy differences are for geometries whose MP2 energies are between 10 000 and 20 000 cm −1 , whereas small differences are seen for MP2 energies of as large as 30 000 cm −1 . The scatter of energies makes it clear that the difference potential is global in nature. Furthermore, as can be seen, the density of points for the difference potential mirrors the density of MP2 energies. This is already an indicator that the difference potential spans a very similar range of configurations as for the original MP2 energies. Further support comes from looking at the distribution of geometries discussed next. Figure 2 shows that the distribution of geometries where LCCSD(T) calculations were performed overlaps the distribution where MP2 calculations were performed, as it must in order to construct a difference potential. The plot axes show the two OH distances, where H is the transferring H atom. The numbering of atoms is shown in the Supporting Information (SI), and because the distribution is permutationally symmetric, we have shown only the upper half by taking O2−H1 to be larger than O3−H1. Here, H1 is the hydrogen atom transferred and O2 and O3 are the symmetric oxygens to which it can be bonded. Note that other distances are changing as well, and this is indicated schematically by the structures shown in the figure. As seen, the LCCSD(T) set is sparse and dispersed. Note also that the clustering observed in both panels is due to the selection procedures adopted by the Meuwly group in the construction of their data set.
In addition to the LCCSD(T)-F12 energies calculated for the PES database, we performed two benchmark calculations at the global minimum and at the saddle point for H-atom transfer. These two calculations found the optimum geometries and energies and also determined the harmonic vibrational frequencies and normal coordinates. While the LCCSD(T)-F12 calculations for just the energy at a single geometry took approximately 30 min using 12 cores of the 2.4 GHz Intel Xeon processors, the full optimization and frequency calculations took on the order of 73 days using the same number of processors. This computational cost certainly underscores the infeasibility of doing even LCCSD(T) calculations for an AcAc PES.
The fit to the difference potential has to take account of the small data set (i.e., a maximum of 2151 energies). Therefore, the number of terms in the PIP basis has to be significantly less than this number to avoid overfitting. With this is mind, we   Figure S1 of the SI, along with an explanation of the meaning of the symmetry designation. This basis contains just 85 PIPs and thus 86 linear coefficients to be determined by standard least-squares regression. This is the smallest PIP we have ever used, and the bonus is that we can examine small training data sets without concerns about overfitting. Furthermore, the evaluation of such a small PIP basis is fast. We note that in order to test the performance of the PES on configurations outside the training set, the largest training set consists of only 90% of the total 2151 LCCSD(T) data points (N = 1935 shown in Table 1); the remaining 10% (every 10th point in the full data set of 2151 points) is reserved for testing. Three smaller data sets are also used to fit the PES: for N = 1075, 717, and 430, every kth point among the 2151 ones is picked as a training point, where k = 2, 3, and 4 for N = 1075, 717, and 430, respectively, while the remaining points not included in the training set are for testing.
Numerous metrics of the performance of the Δ-ML approach are given in Table 1. Beginning with training and testing errors of the ΔV CC−LL PESs using training sets of different sizes, we note again that the testing set consists of points not used for training, and the testing error is for the differences between LCCSD(T) and MP2 electronic energies. It can be seen that the testing error increases monotonically as the number of training points decreases due to smaller coverage by the training data, as expected. However, this increase in the testing error is relatively small.
Next we consider the equilibrium geometries and normalmode frequencies of both the global minimum and the saddle point to H transfer as well as the barrier height. For geometries, we computed the root-mean-square (RMS) difference between the 105 internuclear distances from the PES and direct LCCSD(T)-F12-optimized geometries. For harmonic frequencies, we calculated the mean absolute error (MAE) by comparing frequencies from PESs with direct LCCSD(T)-F12 frequencies. All of these are listed in Table 1 for four PESs with different training sets as well as for the lowlevel MP2 PES. We get excellent agreement for the geometries: in all four Δ-ML PESs, the RMS differences of the 105 internuclear distances between PES and direct LCCSD(T) geometries are around 0.008 and 0.004 Å for GM and SP, respectively, and the RMS differences between V LL and LCCSD(T) geometries are 0.0115 and 0.0026 Å for GM and SP, respectively. Therefore, the geometry of GM is slightly improved using the Δ-ML approach, and the SP geometry is still in good agreement with the LCCSD(T)-F12 geometry despite a slightly increased RMS difference (only 0.0014 Å). A plot of V LL→CC internuclear distances vs direct ab initio distances is shown in Figure S2 of the SI. It is perhaps worth noting that the largest distances (nearly 7 Å) are between the H atoms in the two methyl rotors.
The barrier height of the H-transfer motion on all four corrected PESs, each based on a different training set, is around 1218 cm −1 (3.49 kcal/mol), in excellent agreement with the direct LCCSD(T)-F12 value of 1234 cm −1 (3.54 kcal/mol). The best estimate of this barrier height is 1148 cm −1 (3.29 kcal/mol), based on CCSD(T)-F12/aug-cc-pVTZ single-point calculations at the LCCSD(T)-optimized geometries. Therefore, the Δ-ML PES slightly overestimates the barrier height; nevertheless, it is a significant improvement over the MP2-based PES, 13 which has a barrier height of 745 cm −1 (2.13 kcal/mol).
The energies of seven low-lying stationary points (including the GM and SP) are shown in Table 2. TS(T)-I/II/III represents three transition state saddle points with respect to the torsion of the two methyl rotors, and TS(HT)-I/II represents two higher-order saddle points with imaginary frequencies in both H-transfer motion and methyl torsion. In nearly all cases, the energies of the stationary points are better captured by the V LL→CC PESs than by the V LL PES.
The harmonic frequencies of the global minimum and Htransfer saddle point from the MP2 PES (V LL ), the corrected PES (V LL→CC ) using 1935 training points, and direct LCCSD(T)-F12 calculations are listed in Table S1 of the SI. For most of the modes, the differences between V LL and V LL→CC frequencies are small, but for mode 32 of GM (OH stretch) and the imaginary-frequency mode of the H-transfer SP, the improvement in the Δ-ML PES is significant. Again, the four Δ-ML PESs based on different training sets achieved similar MAE in frequencies (around 13 cm −1 for GM and 25 cm −1 for the H-transfer SP; see Table 1), and that is a significant improvement over the low-level PES, which has MAEs of 17.3 and 35.6 cm −1 for GM and the H-transfer SP, respectively.   These results show that the Δ-ML approach indeed improves the PES and brings it closer to coupled-cluster level of accuracy. This approach significantly improves the barrier height of H transfer, moderately improves the harmonic frequencies of GM and H-transfer SP, and slightly improves the optimized geometries of GM. More importantly, even with a training set as small as 430 points, the corresponding V LL→CC PES is almost as good as the one fitted to 1935 points. Nevertheless, the fit using 1935 points is still our best one in terms of coverage of configurations and testing error, so the results for zero-point energy (ZPE) and H-transfer tunneling splitting shown next are based on this fit.
Diffusion Monte Carlo (DMC) calculations were employed to compute the ground-state tunneling splitting of AcAc. Specifically, the simple unbiased algorithm 28,29 was used to calculate the ground-state energy, while DMC with a fixednode approximation 30 was used to calculate the energy of the excited state with respect to the H-transfer motion.
In the simple unbiased algorithm we use, an ensemble of random walkers is used to represent the nuclear wave function of the molecule. At each step, a random displacement in each degree of freedom is assigned to each walker, and this walker may remain alive (and may give birth to new walkers) or be killed by comparing its potential energy, E i , with a reference energy, E r . For the ground state, the probability of birth or death is given as where Δτ is the step size in imaginary time. In a fixed-node approximation for excited states, in addition to the process described above, any walker that crosses a node is instantly killed. In most cases, the node is unknown in Cartesian coordinates, but for certain modes such as H-transfer in a symmetric double well, a very reasonable approximation can be made for the node as described in detail below. After removing all dead walkers, the reference energy is updated using the equation where τ is the imaginary time; ⟨V(τ)⟩ is the average potential over all the walkers that are alive; N(τ) is the number of live walkers at time τ; and α is a parameter that can control the fluctuations in the number of walkers and the reference energy. Finally, the average of the reference energy over the imaginary time gives an estimate of ZPE (or the energy of the excited state in a fixed-node calculation).
DMC calculations were performed in Cartesian coordinates in full dimensionality. For fixed-node calculations, we assume that the node is r H1O2 = r H1O3 (using the numbering scheme shown in Figure S1 of the SI). Initially, the H1 atom is closer to one O, say O3, so if r H1O3 of a walker becomes larger than r H1O2 , then that walker crosses the node and would be instantly removed. An additional correction was made for the excited state by taking recrossing into account. 30 Ten DMC simulations were performed for each state, and in each simulation, 30 000 walkers were equilibrated for 5000 steps and then were propagated for 50 000 steps to compute the energy, with a step size of 5.0 au. Thus in these simulations, ∼10 10 potential energy evaluations are required; clearly these cannot be done without an efficient PES.
The zero-point energy of the corrected PES using 10 DMC calculations is 26 741 ± 7 cm −1 , while the energy of the excited state for the H-transfer motion from 10 fixed-node DMC calculations is 26 773 ± 10 cm −1 . Therefore, the tunneling splitting is 32 cm −1 with an uncertainty of roughly ±10 cm −1 . As a comparison, the splitting using the MP2-based PES (i.e., V LL ) is 160 cm −1 . Such a significant decrease in the tunneling splitting is expected because the barrier height of V LL→CC , 1218 cm −1 , is significantly higher than the barrier height of V LL , 745 cm −1 .
The ground-state wave function of hydrogen from DMC calculations using the new PES is shown in Figure 3. It can be seen that the largest magnitude of the transferring H atom is closer to one of the O atoms, which is a character of the GM. However, the orientations of the two methyl rotors are the same as for the H-transfer SP. On the basis of the DMC wave function, we cannot draw a definitive conclusion as to whether the ground state of AcAc has a C s or C 2v structure. In fact, the ground-state structure of AcAc is under debate. Meuwly et al. finds a slight preference for the C s structure by the inclusion of a MP2 ZPE difference in the CCSD(T) electronic energies of   Caminati and Grabow support a C 2v structure. 31 We also performed DMC calculations for the singly deuterated isotopologue of AcAc, but the energies of the ground state and excited state are so close that the splitting is smaller than the uncertainty in the DMC calculations. Therefore, we could not obtain a reliable estimate of the tunneling splitting for the deuterated AcAc using DMC.
We also applied an approximate 1d approach to obtain the tunneling splittings. This method has been described in detail previously. 32 It was used in our work on AcAc based on the MP2 PES (V LL ). 13 Briefly, a 1d potential, denoted V(Q im ), which is the minimum-energy path as a function of the imaginary-frequency mode (Q im ) of the H-transfer saddle point, was obtained by optimizing all of the other coordinates at fixed Q im values using the V LL→CC PES except for the methyl rotors, which cannot be described using rectilinear normal coordinates. These are held fixed at the saddle point values all the way along the path. Because of the fixed methyl orientation and fitting error, the barrier height of this 1d Q im path is 1055 cm −1 and it is 179 cm −1 lower than the LCCSD(T) value. Therefore, we "morphed" this 1d potential using the same strategy as described previously 13 so that it gives the correct barrier height (1234 cm −1 ). The two mass-scaled 1d potentials employed (one for H and the other for D transfer) are shown in Figure 4.
The splittings are obtained simply using 1d-DVR calculations 33 of the energies of the ground and first excited states on the morphed V(Q im ) paths, and thus the splitting. We note that this same approach was applied to obtain the tunneling splittings of H-atom and D-atom transfer in malonaldehyde. 32 The results were within roughly 10% of the rigorous diffusion Monte Carlo splittings.
Using this 1d approach, the ground-state tunneling splittings are 37.9 and 8.2 cm −1 for H and D, respectively. The Hsplitting is in reasonably good agreement with the 32 cm −1 obtained from the DMC calculation, and the D-splitting is also consistent with the fact that it is smaller than the uncertainty of the DMC calculations.
It is clearly seen that the Δ-ML approach can indeed bring the PES closer to the coupled-cluster level of accuracy, especially as applied to the barrier height of the H-atom transfer. The improvements in geometries and harmonic frequencies are relatively smaller, as the MP2 results are already quite close to the new LCCSD(T) ones presented here.
Of course these improvements are not achieved without extra cost. First, there is an extra cost to compute the LCCSD(T)-F12 energies for 2151 configurations. A singlepoint calculation of the energy takes about 30 min using 12 cores of the 2.4 GHz Intel Xeon processor, and all 2151 points can be completed within 1 week using 7 nodes for these computations. As we have shown above, far fewer points can be used to obtain a high-quality correction PES, so this cost is affordable and minor.
Perhaps of more importance is the extra cost in evaluating the corrected PES. In the Δ-ML approach, the extra cost is the calculation of the energy correction, ΔV CC−LL . For AcAc, the ΔV CC−LL PES uses a maximum polynomial order of 2 and costs about 10% of the V LL PES, but this additional cost is also a small price to pay for bringing about the accuracy of the PES, especially the barrier height for the H-atom transfer, to near the coupled-cluster level.
The success in bringing an MP2-based full-dimensional PES for 15-atom acetyacetone to coupled-cluster quality is very encouraging. Given that a small number (around 500 to 2000) of coupled-cluster energies were needed for the correction makes it clear that the approach should be readily applicable to molecules with more than 10 atoms. The coupled-cluster approach used was the relatively efficient local method LCCSD(T) available in the 2015 version of Molpro that we use. 34 Other efficient CCSD(T) methods are also available in both Molpro and other software packages.
Overall, the success of the current application of the Δ-ML method to 15-atom AcAc and the applications reported earlier on smaller molecules and application to cis and trans-Nmethylacetamide suggests that the method can have wide applicability and ease of use. Since the correction PES is not localized around a minimum or a reaction path, it can be used in general anharmonic vibrational analyses of polyatomic molecules. The present application to full-dimensional DMC calculations is already evidence of this.
Finally, it is also worth commenting on the current Δ-ML approach and the recent application of transfer learning (TL) by Meuwly and co-workers 14 to MP2-based PESs for AcAc. The MP2-based PES we reported used a slightly extended database of MP2 energies and gradients from that group. Thus, the PES we reported is not the same as the earlier one. However, they are similar (e.g., the H-atom transfer barrier heights are 2.13 and 2.17 kcal/mol for the PIP PES and the neural-network (NN) PES, respectively). These are in very good agreement with the direct MP2 result of 2.18 kcal/mol. Several TL NN models, based on random training data sets (PNO-LCCSD(T) energies) of different sizes, were considered by Meuwly and co-workers. 14 In Table 3 of that paper, results both from a single TL NN model and from an average of several TL NN models were given at the optimized geometries of each model. A TL-NN model gave a barrier height of 0.92 kcal/mol using 100 high-level energies, 2.4 kcal/ mol from a single TL-NN model, 1.80 kcal/mol averaged from several TL-NN models with 1000 energies, and 2.66 kcal/mol from a single TL-NN model using 5000 energies. Barrier heights of 3.31 and 3.32 kcal/mol were obtained with a single and multiple TL-NN models using 15 000 energies; these are in excellent agreement with the benchmark barrier height of 3.25 kcal/mol. From these results, two apparent conclusions emerge. The first is that the TL-NN may produce a worse result (i.e., a lower barrier height) than the original NN PES, and the second is that roughly 15 000 high-level energies are needed to obtain a TL-NN PES with an accurate barrier height. (These authors did note some improvement in the TL-NN model using 1000 energies by the addition of 100 energies along the minimum-energy path (i.e., the barrier height increased from 1.8 to 2.72 kcal/mol).
Thus, on the basis of the above, it appears that the Δ-ML approach performs well and is a reasonable alternative to TL for this example. However, unlike the TL-NN approach, the Δ-ML approach we applied does not produce worse results than the low-level PES. We attribute this to the fact that the difference potential is small relative to the low-level and highlevel potentials. If this is not true (and this can be checked, of course), then the current approach would probably require a larger database of high-level energies to achieve a satisfactory result. These are early conclusions, and more work on both approaches is clearly warranted.
Atom-numbering scheme; table of benchmark harmonic frequencies; and fitting precision plots (PDF)