Bridging semiempirical and ab initio QM/MM potentials by Gaussian process regression and its sparse variants for free energy simulation

Abstract

Free energy simulations that employ combined quantum mechanical and molecular mechanical (QM/MM) potentials at ab initio QM (AI) levels are computationally highly demanding. Here, we present a machine-learning-facilitated approach for obtaining AI/MM-quality free energy profiles at the cost of efficient semiempirical QM/MM (SE/MM) methods. Specifically, we use Gaussian process regression (GPR) to learn the potential energy corrections needed for an SE/MM level to match an AI/MM target along the minimum free energy path (MFEP). Force modification using gradients of the GPR potential allows us to improve configurational sampling and update the MFEP. To adaptively train our model, we further employ the sparse variational GP (SVGP) and streaming sparse GPR (SSGPR) methods, which efficiently incorporate previous sample information without significantly increasing the training data size. We applied the QM-(SS)GPR/MM method to the solution-phase SN2 Menshutkin reaction, NH3+CH3Cl→CH3NH3++Cl-, using AM1/MM and B3LYP/6-31+G(d,p)/MM as the base and target levels, respectively. For 4000 configurations sampled along the MFEP, the iteratively optimized AM1-SSGPR-4/MM model reduces the energy error in AM1/MM from 18.2 to 4.4 kcal/mol. Although not explicitly fitting forces, our method also reduces the key internal force errors from 25.5 to 11.1 kcal/mol/Å and from 30.2 to 10.3 kcal/mol/Å for the N-C and C-Cl bonds, respectively. Compared to the uncorrected simulations, the AM1-SSGPR-4/MM method lowers the predicted free energy barrier from 28.7 to 11.7 kcal/mol and decreases the reaction free energy from -12.4 to -41.9 kcal/mol, bringing these results into closer agreement with their AI/MM and experimental benchmarks.

FIG. 1.

Menshutkin reaction between ammonia and methyl chloride.

FIG. 2.

AM1/MM and indirect B3LYP/MM free energy profiles for the Menshutkin reaction. The AM1/MM free energy profile was obtained along the string MFEP determined at the same level (with α = 0 for reactant and 1 for product). The indirect free energy profile was obtained by free energy perturbation (FEP) from the AM1/MM to the B3LYP/6-31+G(d,p)/MM level based on the configurations sampled along the AM1/MM MFEP. Even with AI/MM-based free energy corrections, the indirect FEP results still do not adequately reproduce the experimental results. Details on the calculations of the error bars can be found in the supplementary material, Sec. S4.

FIG. 3.

The GPR-predicted energy corrections (ΔU^pred) plotted against the reference energy differences (ΔU^ref) between the AM1/MM and B3LYP/6-31+G(d,p)/MM levels for configurations sampled along the AM/MM MFEP. The color scale shows the number of configurations found in each region.

FIG. 4.

(a) The AM1/MM and AM1-GPR/MM MFEPs. (b) The AM1/MM and AM1-GPR/MM free energy profiles.

FIG. 5.

(a) The AM1/MM and AM1-SSGPR/MM MFEPs. (b) The AM1/MM and AM1-SSGPR/MM free energy profiles.

FIG. 6.

Internal force differences between the B3LYP/6-31+G(d,p)/MM and AM1-SSGPR-4/MM levels; the force differences on the C–Cl bond (a) and N–C bond (b) were evaluated based on configurations sampled along the AM1-SSGPR-4/MM MFEP.