Multireference Generalization of the Weighted Thermodynamic Perturbation Method

Abstract

We describe the generalized weighted thermodynamic perturbation (gwTP) method for estimating the free energy surface of an expensive "high-level" potential energy function from the umbrella sampling performed with multiple inexpensive "low-level" reference potentials. The gwTP method is a generalization of the weighted thermodynamic perturbation (wTP) method developed by Li and co-workers [J. Chem. Theory Comput. 2018, 14, 5583-5596] that uses a single "low-level" reference potential. The gwTP method offers new possibilities in model design whereby the sampling generated from several low-level potentials may be combined (e.g., specific reaction parameter models that might have variable accuracy at different stages of a multistep reaction). The gwTP method is especially well suited for use with machine learning potentials (MLPs) that are trained against computationally expensive ab initio quantum mechanical/molecular mechanical (QM/MM) energies and forces using active learning procedures that naturally produce multiple distinct neural network potentials. Simulations can be performed with greater sampling using the fast MLPs and then corrected to the ab initio level using gwTP. The capabilities of the gwTP method are demonstrated by creating reference potentials based on the MNDO/d and DFTB2/MIO semiempirical models supplemented with the "range-corrected deep potential" (DPRc). The DPRc parameters are trained to ab initio QM/MM data, and the potentials are used to calculate the free energy surface of stepwise mechanisms for nonenzymatic RNA 2'-O-transesterification model reactions. The extended sampling made possible by the reference potentials allows one to identify unequilibrated portions of the simulations that are not always evident from the short time scale commonly used with ab initio QM/MM potentials. We show that the reference potential approach can yield more accurate ab initio free energy predictions than the wTP method or what can be reasonably afforded from explicit ab initio QM/MM sampling.

Figure 1:

Non-enzymatic RNA 2′-O-transesterification model reactions explored in this work (atomic numbering of 2′ and 5′ positions in the model systems reflect their analogous positions in RNA). (a) The ethylene phosphate model reaction with a methoxide leaving group. (b) The native model reaction used in ref. with an ethoxide leaving group. The reaction coordinate used in this work is ξ_PT = R_P-O5’ − R_P-O2’.

Figure 2:

Free energy surfaces (ΔA) of the ethylene phosphate model reaction (Figure 1a) calculated with MBAR and wTP analysis. Reweighting entropies (RE) are shown in the subplots below each FES. Parts (a), (e), (i), and (m) are the 4 MNDO/d QM/MM+DPRc potentials (ML0, ML1, ML2 and ML3) evaluated from MBAR analysis. Parts (b), (f), (j), and (n) are their corresponding wTP estimates of the PBE0/6–31G* FES. Parts (c), (d), (g), (h), (k), (l), (o), and (p) are the reweighting entropy of the corresponding MBAR and wTP analysis. The error bars are 95% confidence intervals.

Figure 3:

Free energy surfaces (ΔA) of the ethylene phosphate model reaction (Figure 1) calculated with gwTP analysis. Reweighting entropies (RE) are shown in the subplots below each FES. (a) The gwTP FES estimated from the combined umbrella sampling produced from the 4 MNDO/d QM/MM+DPRc potentials (ML0, ML1, ML2 and ML3). (b) The PBE0/6–31G* surface estimated from gwTP analysis of all available PBE0/6–31G* and MNDO/d QM/MM+DPRc umbrella sampling. Parts (c) and (d) are the reweighting entropy of the corresponding MBAR and gwTP analysis. The error bars are 95% confidence intervals.

Figure 4:

Free energy surfaces (ΔA) of the native model reaction (Figure 1). (a) Comparison of the PBE0/6–31G* FES to the average FES of the 4 DFTB2/MIO QM/MM+DPRc potentials (ML0, ML1, ML2 and ML3). The listed times are the total amount of sampling per umbrella window used in the FES calculation. (b) Comparison of the PBE0/6–31G* FES to those predicted from gwTP analysis. The ΔA(PBE0;ML*) gwTP surfaces are evaluated using the sampling from all 4 DPRc potentials. The ΔA(PBE0;All) surface is calculated from the 1.2 ns/window of DPRc sampling and the 100 ps/window PBE0/6–31G* sampling. The error bars are 95% confidence intervals. The RE values shown in parts (c) and (d) are the reweighting entropies.

Figure 5:

Free energy surfaces (ΔA) of the native model reaction (Figure 1). (a) Comparison of the PBE0/6–31G* FES to the average FES of the 4 DFTB2/MIO QM/MM+DPRc potentials (ML0, ML1, ML2 and ML3). (b) Comparison of the PBE0/6–31G* FES to those predicted from gwTP analysis. The surfaces labeled “last 600 ps” refer to aggregate sampling taken from the last 50 ps of the 12 DFTB2/MIO QM/MM+DPRc simulations. The error bars are 95% confidence intervals. The RE values shown in parts (c) and (d) are the reweighting entropies.