Assessing the Accuracy of Machine Learning Thermodynamic Perturbation Theory: Density Functional Theory and Beyond

Machine learning thermodynamic perturbation theory (MLPT) is a promising approach to compute finite temperature properties when the goal is to compare several different levels of ab initio theory and/or to apply highly expensive computational methods. Indeed, starting from a production molecular dynamics trajectory, this method can estimate properties at one or more target levels of theory from only a small number of additional fixed-geometry calculations, which are used to train a machine learning model. However, as MLPT is based on thermodynamic perturbation theory (TPT), inaccuracies might arise when the starting point trajectory samples a configurational space which has a small overlap with that of the target approximations of interest. By considering case studies of molecules adsorbed in zeolites and several different density functional theory approximations, in this work we assess the accuracy of MLPT for ensemble total energies and enthalpies of adsorption. The problematic cases that were found are analyzed and it is shown that, even without knowing exact reference results, pathological cases for MLPT can be detected by considering a coefficient that measures the statistical imbalance induced by the TPT reweighting. For the most pathological examples we recover target level results within chemical accuracy by applying a machine learning-based Monte Carlo (MLMC) resampling. Finally, based on the ideas developed in this work, we assess and confirm the accuracy of recently published MLPT-based enthalpies of adsorption at the random phase approximation level, whose high computational cost would completely hinder a direct molecular dynamics simulation.

The use of machine learning (ML) techniques could be highly beneficial for MD simulations and, since the seminal work of Behler and Parrinello (BP) in 2007 22 , has seen an increasing popularity in this field [23][24][25][26] . While keeping a level of accuracy comparable to ab initio calculations, ML approaches can be used to replace most of the expensive quantum mechanical calculations with numerically cheap predictions. This allows for an increase in the system size and timescale that would be normally accessible by traditional ab initio MD simulations. However, ML models typically employed in MD simulations require a significant amount of data to be trained and this represent an issue for the most expensive quantum mechanical approximations.
In a recent work, some of us applied a scheme that couples machine learning techniques with thermodynamic perturbation theory (MLPT) to compute enthalpies of adsorption of molecules in the zeolite chabazite at the RPA level of theory 27 (the same methodology can be used in the context of free energies of activation 28,29 ). Achieving these results with a brute force molecular dynamics simulation would be completely out of reach. Indeed, by considering only the cost of the RPA energy calculations for the zeolite, completing a 200000 steps MD simulation would require 65M CPU hours and more then 110 years on 64 cores at 2.6 GHz (this estimate does not include the significant additional cost of computing RPA forces 30 ). With our MLPT technique as few as 10 single point RPA energy calculations were sufficient to train a machine learning model that was then used to predict the RPA energy for several other configurations 27 .
In this previous work the van der Waals (vdW) corrected PBE+D2 functional 31,32 was used to generate a production MD trajectory (in general, a numerically inexpensive approximation should be chosen as production method). Based on the ∆-ML method 33 , a model was thereafter trained to predict the difference between RPA and PBE+D2 energies. Following the thermodynamic perturbation approach 34-36 the energy differences were subsequently used to reweight the statistical weight of the PBE+D2 configurations in order to obtain the RPA canonical distribution and compute the RPA ensemble energies and enthalpies (in this context the RPA is considered the "target" level of theory). Thanks to the inexpensive predictions of the ML model, this technique involves a computational cost that is several orders of magnitude smaller than a full RPA MD simulation.
While the error involved in the ML procedure is rather well controlled 27,28 , the application of thermodynamic perturbation theory (TPT) might be at the origin of a bias in the results estimated at the target level of theory. Indeed, the configurational space sampled by the production MD might have suboptimal superposition with the target configurational space 34,35 . In certain cases this might lead to a strong loss of statistical significance when the contributions of the production configurations are reweighted to obtain the target distribution (namely only few configurations could contribute to the whole target level statistics).
In this work, we access the accuracy of the MLPT approach using the adsorption of molecules in zeolites as a test case. Indeed, many separation and adsorption processes re-quire the improvement of the dedicated materials, and DFT calculations often provide suitable suggestions of optimized formulations in this regard [37][38][39][40] . However, additional efforts have to be done to find a better compromise between accuracy of prediction and calculation cost. By considering five different DFT functionals, including generalized gradient approximation (GGA), meta-GGA, and non-local vdW corrected functionals, we perform full MD simulations to provide "exact" references values for the ensemble energies and enthalpies of adsorption. We then use each one of these five functionals as a production method and apply MLPT to obtain the energy/enthalpy estimates for the remaining four functionals, which are considered as target approximations.
By comparing MLPT estimates with exact reference results the level of accuracy of this approach is established. Some pathological cases are found, for example when GGA functionals are used as starting point for meta-GGA functionals and vice versa (in certain cases the MLPT estimate of the ensemble energy can deviate by more than 8 kcal/mol). By using machine learning dimensionality reduction algorithms, we qualitatively visualize the relative distributions of the configurations from the different MD simulations. This analysis confirms that the most problematic cases can be ascribed to a poor superposition of the configurational spaces. In order to detect possible failure of MLPT we propose a diagnostic test, the I w index, which is evaluated from ML energy predictions and does not rely on the knowledge of any exact reference results.
We then propose a scheme to significantly improve results even when the I w coefficient is close to 0 (lowest production-target superposition). This approach, denoted as MLMC, is based on a Monte Carlo resampling of the target configurational space that reuse the machine learning model already trained for MLPT. Without any additional target level calculations this scheme decreases even the largest deviations within the threshold of chemical accuracy.
Finally, we apply the tools developed in this work to analyze our previous results based on the RPA 27 , where reference values cannot be produced. The relatively large values of the I w index hints that the PBE+D2 functional provides a reliable starting point for RPA target properties. The results are stable even if a full MLMC resampling of the configurational space is performed. This shows that, with a proper choice or the production approximation, MLPT allows for a quick and accurate estimate of target level properties. The MLPT approach opens the way to a more systematic application of accurate but expensive DFT and quantum chemical approximations in finite temperature simulations of materials.
The layout of the paper is organized as follows. In the following section, Sec. 2, we discuss the methodological part of this work and provide computational details. In Sec. 3, we present numerical results and analyze them. Sec. 4 contains our conclusions.

Methodology
The enthalpies of adsorption(∆ ads H) of the CH 4 molecule in protonated chabazite (HCHAB) and CO 2 molecule in siliceous chabazite (SiCHAB) were investigated by employing ab initio molecular dynamics (AIMD) and machine learning thermodynamic perturbation theory [27][28][29] (MLPT). The models of the adsorbed molecules are shown in Fig. (1). The enthalpy of adsorption is defined by where E is the internal energy computed as ensemble average of potential energy via the AIMD simulations or MLPT (M denotes the molecules, Z the zeolites, and M@Z the adsorbed system), k B is the Boltzmann constant, and T is the system temperature (equal to 300K in all our simulations). The goal of this work is to benchmark the accuracy of MLPT by comparing its estimates with "exact" reference values obtained from full MD simulations. We chose five DFT functionals in order to span a range of different characteristics but keeping also into account their numerical cost, as long and stable MD trajectories have to be produced for reference. non-local vdW corrections, SCAN+rVV10 46 . We will also discuss results based on the random phase approximation (RPA), as implemented in the VASP code [47][48][49] ; due to the high numerical cost a direct MD simulation at this level of theory is completely out of reach and the MLPT method becomes instrumental to estimate finite temperature properties.
All AIMD simulations were performed within the NVT ensemble, whereby temperature was maintained by Andersen 50 thermostat with collision probability of 0.05. An integration time step of 0.5 fs was considered and the total simulation time of each MD run was 100 ps, i.e. ∼200·10 3 configurations were generated with each method for each system investigated in this work. The initial part of 10 ps of each trajectory was considered as an equilibration period and the corresponding data were discarded. The remaining part of the trajectory was The results produced by AIMD are used as a reference for our MLPT method. Within this approach an MD simulation is first performed at a certain production level of theory and then "corrected" to obtain results at one or multiple target levels of theory. One of the advantages of MLPT is that an estimate of a finite-temperature property for the target method can be obtained with only a few tens or few hundreds of single point calculations 27 .
Out of about 1900·10 3 configurations generated in production runs, a reduced sample of about 19·10 3 evenly separated configurations (X i ) and corresponding potential energies obtained with the production method (E i ) are selected to be used for the predictions of the energies at the new target functional level. We note that this configuration reduction is not strictly necessary but, since MD configurations close in time are correlated, this procedure speeds up the predictions of the ML model without compromising the accuracy of prediction.
In order to compute target level ensemble energies E , MLPT applies TPT to reweight the production statistical distribution 35 where E i denotes the target level energies of the ith configuration, to consider multiple target level theories or highly expensive approximations (e.g., the RPA).
As discussed in the next sections MLPT does not provide satisfactory results in certain cases. This happens for the estimate of SCAN and SCAN+rVV10 adsorption enthalpies from PBE production calculations. This issue is related to the unsatisfactory superposition of the target and production configurational spaces. To improve the results in these cases we have reused the ML model from the MLPT application to resample the target level configurational space using the Monte Carlo (MC) algorithm [61][62][63][64] . We will denote this procedure as MLMC.
The MC approach used to sample the canonical ensemble is equivalent to MD but is simpler to implement and does not require the calculation of forces (for certain methods forces are not available in most solid state physics software implementations). By using our in-house software, the MLMC approach recursively sample new configurations following these steps: (1) A new geometric configuration is randomly generated; (2) A production level ab initio calculation is performed to obtain E i for this new configuration; (3)  we will discuss the accuracy of MLPT for the case of computationally affordable DFT functionals, for which an "exact" reference can be easily obtained. At the end of this section we will consider the RPA enthalpies of adsorption.
In relative errors with respect to experiment higher than 35%. However, the main purpose of this work is not to compare the accuracy of different functionals compared to experiment but to assess the performance and limitations of MLPT to reproduce exact AIMD results obtained using different quantum mechanical approximations.   As the enthalpies of adsorption are computed from differences of internal energies of interacting and non-interacting systems, fortuitous error cancellations might arise that could hide the real accuracy of MLPT predictions. To shed some light onto this problem, we report in Table (3) MLPT total ensemble energies for the individual components involved in the adsorption process (adsorbed system, zeolite alone, molecule alone). These results confirm the conclusions previously drawn from Table (2), with the predictions made for the cases where the production and target methods are both from the same group of functionals (i.e., either PBE/PBE+D2/vdW-DF2 or SCAN/SCAN+rVV10) being well within the chemical accuracy.
MLPT becomes unreliable when predictions mix these two groups, with deviations that often reach several kcal/mol (these observations do not hold for the standalone CH 4 and CO 2 molecules, whose deviations are always below 0.1 kcal/mol). The worst performance is confirmed to correspond to the prediction of SCAN and SCAN+rVV10 from the PBE production trajectory of CH 4 @HCHAB. For these cases, we notice that the largest deviations in the enthalpies of adsorption (>8 kcal/mol) are not only due to the low accuracy of the MLPT ensemble energy estimates for CH 4 @HCHAB and HCHAB, but also to the opposite sign in deviations of predictions made for these two systems, which leads to an error accumulation when evaluating Eq. (1).
While MLPT allows, in principle, for a quick evaluation of finite-temperature properties at one or more target levels of theory from a single MD production run (typically based on the most computationally inexpensive approximation), the results in Tables (2-3) show that such a strategy does not always allow one to achieve the required level of accuracy.
Accordingly, a special care should be taken in choosing a production method suitable for the target approximation(s) of interest. The failure of MLPT in certain cases can be explained by the limitations of the thermodynamic perturbation theory (TPT) itself, rather than by an inaccurate ML model. Indeed, if the high-probability configurational spaces of the production and target approximations do not overlap sufficiently, the TPT has low predictive power 35 . In the most pathological cases it can happen that this overlap is so poor that in the reweighting procedure (Eq. (2)) only one or few individual configurations effectively contribute to the ensemble average. In order to analyze qualitatively this behavior we consider all the production MD simulations based on the 5 different functionals. By using the t-distributed stochastic neighbor embedding (t-SNE) algorithm 67,68 we present in Fig. (2) two-dimensional visualizations of the high-dimensional configurational space spanned by these trajectories (to improve the readability, figures were created using only 500 uncorrelated structures selected from each trajectory). To be consistent with the ML learning algorithm used in this study for the regression, the t-SNE approach was applied using a definition of distance D between two configurations χ A and χ B based on the normalized SOAP kernel K: this definition also inherits some of the properties of the SOAP kernel, such as the rotational and translational invariance required in materials and molecular modelling.
It can be noticed in Fig. (2) that for CH 4 @HCHAB, CO 2 @SiCHAB, HCHAB, and SiCHAB the trajectories generated by the two groups of functionals (PBE/PBE+D2/vdW-DF2 and SCAN/SCAN+rVV10) form two clusters with very limited overlap. This shows that the configurational spaces spanned by the two groups of approximations are to a large extent different, which is at the origin of the poor performance of MLPT in the cases discussed above. For the gas phase molecules CH 4 and CO 2 , the configurational space has a much simpler structure and the five different approximations produce trajectories which always largely overlap; this is also consistent with the results in Table (3), where MLPT is highly accurate for all molecular energies. The t-SNE visualization provides useful insights on the potential failures of TPT and, consequently, MLPT. However, a diagnostic approach that requires the full MD trajectories is highly impractical in reality, since the goal of MLPT is indeed to completely avoid the MD simulation at the target level of theory. In order to overcome this difficulty and introduce a tool to analyze and predict the failures of MLPT, we follow the approach introduced in Ref. 28 and introduce the I w index. By considering the weight factor (w i ) in Eq. (2) where M is the total number of configurations sampled in the production run and N is the lowest integer which satisfies the condition where the weights w i have been sorted in ascending order. In practice, the I w index reflects how many configurations contribute after the reweighting in Eq. (2) to half of the total statistical weight.
In the case of an extremely small overlap between configurational spaces sampled by the production and target approximations, we would have N ≈ M (since most weights w i would be close to 0, N has to approach M to satisfy the condition in Eq. (5)) and I w ≈ 0.
This corresponds to a significant loss in the statistical significance since the large number of configurations sampled in the production run is reduced to a few individual structures. In the opposite case of excellent overlap we would have N ≈ M/2 and I w ≈ 0.5; if for example the production and target approximations provides energies that differ only by a constant we would obtain I w = 0.5.
The I w values for all systems investigated in this work are presented in Table ( 4); in order to simplify the analysis, we are also reporting in parenthesis the deviations of the MLPT ensemble energies with respect to the exact MD references (same as in Table (3)).
All the the I w values rounded to the value 0.00 actually correspond to values that are always smaller than 0.003. As expected the I w index is equal to 0.5 when the production and target methods are the same. The values of I w allow for an interpretation of the deviations previously reported in Table (3). All the largest errors of MLPT (often significantly above chemical accuracy) correspond to I w = 0.00, namely to a small correlation/overlap of production and target methods. In few cases, small deviations are found for systems with I w = 0.00 (for example the MLPT evaluation of the SCAN ensemble energy from PBE+D2 for HCHAB) but these results must be considered as coincidental, whereby the reweighting of a very small number of configurations led, by chance, to a satisfactory reconstruction of the target statistical distribution.
As discussed above, the MLPT estimates become reliable for the production and target method combinations from within the two groups of functionals PBE/PBE+D2/vdW-DF2 and SCAN/SCAN+rVV10. As shown in Table (      SCAN+rVV10 (MLMC(SCAN+rVV10)) functionals using PBE as production method.
We have shown insofar that the MLMC can significantly improve the sampling of the geometries at the target level. We now show in Once again, the largest improvements have been observed for the adsorbate with substrate systems, followed by the clean zeolites, while the effect of MLPT or MLMC is negligible for the isolated molecular systems for which most of the observed deviations are close to zero and probably are related to numerical errors. In Tables (5-6), all the predictions of SCAN/SCAN+rVV10 energies from PBE were significantly improved by the MLMC method.
The general improvement reduced the relative internal energy deviations below 1.0 kcal/mol.

Assessing the accuracy of MLPT for the random phase approximation
In a previous work (Ref. 27 ) some of us applied MLPT to compute enthalpies of adsorption at the RPA level of theory. Beyond the two systems studied in this work, also CH 4 @CHAB and CO 2 @HCHAB were considered in this previous paper. The RPA enthalpies of adsorption showed an excellent agreement with experiment and for the molecules adsorbed in HCHAB zeolite a significant improvement with respect to the production method (PBE+D2) was found (the results for ensemble energies and enthalpies of adsorption will not be repeated here, where we will rather discuss their reliability). Because of the high computational cost involved in the RPA calculations in this case, it is not possible to produce reference trajectories and use them to repeat some of the previous analysis (indeed MLPT is meant to make possible finite-temperature calculations involving computationally expensive approximations).
However, the I w index can be straightforwardly obtained and we present its values in

Conclusions
We assessed the accuracy of the machine learning thermodynamic perturbation theory approach by computing enthalpies of adsorption of molecules in zeolites. By considering a set of five DFT functionals with different characteristics, MLPT estimates are compared to reference results produced from full MD simulations. The largest deviations are found when the configurational space accessible to the target level has a small overlap with that sampled by the production method. Even when an exact reference result is not available, the most problematic cases can be detected by using a diagnostic test, the I w index. In these cases the results can be significantly improved by coupling the Monte Carlo approach to the MLPT model to resample the configurational space at the target level of theory which is, however, performed at the cost of the production method.
Finally, the analysis is extended to some recently published MLPT results that used the RPA as target approximation 27 ; in this case the high computational cost of the RPA makes the generation of exact reference ensemble energies and enthalpies completely unpractical.
The application of the I w diagnostic test and of the MC resampling highlights the full reliability of those previous results, which were based on PBE+D2 production MD.
The main advantage of the MLPT approach stands in the possibility of evaluating finite temperature properties at a certain target level of theory by performing only a minimal number of single-point calculations. Additionally, MLPT is also suitable for approximations where the computation of gradients is not available or adds a significant computational cost 30 .
This approach opens the possibility of systematically applying high-accuracy/high-cost ab initio methodologies to achieve a new level of predictive power in materials simulations.