Doubly Polarized QM/MM with Machine Learning Chaperone Polarizability

Abstract

A major shortcoming of semiempirical (SE) molecular orbital methods is their severe underestimation of molecular polarizability compared with experimental and ab initio (AI) benchmark data. In a combined quantum mechanical and molecular mechanical (QM/MM) treatment of solution-phase reactions, solute described by SE methods therefore tends to generate inadequate electronic polarization response to solvent electric fields, which often leads to large errors in free energy profiles. To address this problem, here we present a hybrid framework that improves the response property of SE/MM methods through high-level molecular-polarizability fitting. Specifically, we place on QM atoms a set of corrective polarizabilities (referred to as chaperone polarizabilities), whose magnitudes are determined from machine learning (ML) to reproduce the condensed-phase AI molecular polarizability along the minimum free energy path. These chaperone polarizabilities are then used in a machinery similar to a polarizable force field calculation to compensate for the missing polarization energy in the conventional SE/MM simulations. Because QM atoms in this treatment host SE wave functions as well as classical polarizabilities, both polarized by MM electric fields, we name this method doubly polarized QM/MM (dp-QM/MM). We demonstrate the new method on the free energy simulations of the Menshutkin reaction in water. Using AM1/MM as a base method, we show that ML chaperones greatly reduce the error in the solute molecular polarizability from 6.78 to 0.03 ų with respect to the density functional theory benchmark. The chaperone correction leads to ∼10 kcal/mol of additional polarization energy in the product region, bringing the simulated free energy profiles to closer agreement with the experimental results. Furthermore, the solute-solvent radial distribution functions show that the chaperone polarizabilities modify the free energy profiles through enhanced solvation corrections when the system evolves from the charge-neutral reactant state to the charge-separated transition and product states. These results suggest that the dp-QM/MM method, enabled by ML chaperone polarizabilities, provides a very physical remedy for the underpolarization problem in SE/MM-based free energy simulations.

Figure 1.

Solution-phase molecular polarizability as a function of the reaction coordinate for the Menshutkin reaction: AM1 (squares), B3LYP/aug-cc-pVTZ (circles), and their difference (triangles). The means (solid curves) and standard deviations (vertical bars) are computed based on the samples within each string image (see SI.5 for the tabulated statistical distributions).

Figure 2.

Regressions of the molecular polarizabilities from AM1 (square with a dashed line) and from the chaperone-corrected AM1 (AM1+Δα^C; circles with a solid line) against those from B3LYP/aug-cc-pVTZ; the corresponding root-mean-square errors (RMSEs) relative to the B3LYP reference values are also shown.

Figure 3.

Free energy profiles as a function of the reaction coordinate for the Menshutkin reaction: AM1/MM (dashed line) and dp-AM1/MM with polarizabilities corrected to B3LYP/aug-cc-pVTZ (solid line). The error bars relative to the free energy in the reactant state (α = 0) along the string MFEP are estimated using a procedure developed by Zhu and Hummer (Ref. 70), slightly modified for nonuniform collective-variable grids (see SI.1 for details)

Figure 4.

The minimum free energy path (MFEPs) as a function of the collective variables, i.e., the C-Cl and N-C bonds: AM1/MM (dashed line) and the dp-AM1/MM with chaperone polarizabilities corrected to the B3LYP/aug-cc-pVTZ level (solid line). The locations of free energy transition states are also marked: AM1/MM (open square) and dp-AM1/MM (open circle).

Figure 5.

Atomic chaperone polarizabilities as a function of the reaction coordinate for each solute atom in the Menshutkin reaction.

Figure 6.

Convergence of atomic polarizability on the chlorine atom with respect to basis sets and AI methods.

Figure 7.

Radial distribution functions between the chlorine atom and water oxygens (O_w) in the reactant (R), transition-state (TS), and product (P) regions of the Menshutkin reaction: AM1/MM (dashed line) and dp-AM1/MM with polarizabilities corrected to the B3LYP/aug-cc-pVTZ level (solid line).

Figure 8.

Decomposition of polarization energy to each water molecule around the chlorine atom: AM1/MM (squares), chaperone polarization energy (triangles), and B3LYP/aug-cc-pVTZ/MM (circles).

Scheme 1.

Schematic representation of the Menshutkin reaction from the charge-neutral reactant state to charge-separated product state (NH₃ + CH₃Cl → NH₃CH₃⁺ + Cl⁻).

Scheme 2.

Topology of artificial neural network containing two CV inputs, ten hidden neurons, and nine atomic chaperone polarizabilities in output layer