QM/MM free energy simulations: recent progress and challenges

Abstract

Due to the higher computational cost relative to pure molecular mechanical (MM) simulations, hybrid quantum mechanical/molecular mechanical (QM/MM) free energy simulations particularly require a careful consideration of balancing computational cost and accuracy. Here we review several recent developments in free energy methods most relevant to QM/MM simulations and discuss several topics motivated by these developments using simple but informative examples that involve processes in water. For chemical reactions, we highlight the value of invoking enhanced sampling technique (e.g., replica-exchange) in umbrella sampling calculations and the value of including collective environmental variables (e.g., hydration level) in metadynamics simulations; we also illustrate the sensitivity of string calculations, especially free energy along the path, to various parameters in the computation. Alchemical free energy simulations with a specific thermodynamic cycle are used to probe the effect of including the first solvation shell into the QM region when computing solvation free energies. For cases where high-level QM/MM potential functions are needed, we analyze two different approaches: the QM/MM-MFEP method of Yang and co-workers and perturbative correction to low-level QM/MM free energy results. For the examples analyzed here, both approaches seem productive although care needs to be exercised when analyzing the perturbative corrections.

Figure 1

Thermodynamic cycles used to compute the solvation free energy of a QM solute in (a) a pure MM environment and (b) when nearby solvent molecules are also described at the QM level. In both cases, a pure MM solvation free energy calculation (bottom horizontal process) is used as a reference; to obtain the desired QM solvation free energies, various vertical processes need to be studied in which the solute (colored red) or a micro droplet (colored yellow) is converted between a MM and a QM model either in the gas phase or in an MM solution environment (in blue).

Figure 2

Two different routes to compute the free energy change as the potential function is switched between a low level and a high level. (a–b) illustrate the routes for alchemical free energy calaculations and (c–d) illustrate the corresponding calculations along a path (e.g., an optimized string). In the NBB route, a low-level window (e.g., ${\lambda }_1^L$) is connected with a neighboring high-level (e.g., ${\lambda }_0^H$) by re-weighting the ${\lambda }_0^L$ data. In the FEP(+BAR) route, free energy change due to a switch in the potential function is calculated only for the same (e.g., λ₀) window; different λ/s windows (e.g., λ₀/λ₁ or s₀/s₁) are connected with the low-level potential function using BAR. See text for the corresponding equations for the free energy changes (compare Eqs. 21 and 23).

Figure 3

Results for the proton transfer in a model channel (see Fig. 4) from 1D REUS simulation along the ζ coordinate. (a) Convergence of the PMF as a function of simulation time per replica; (b) Illustration of the exchange behavior of replicas during the REUS simulation, using replica 2 as an example.

Figure 4

2D metadynamics simulations for proton transfer in a model channel to illustrate the coupling between hydration level change and proton transfer. (a)–(b) Snapshots for the reactant and transition state regions, respectively. The “bulk” MM water molecules are shown in blue, and the QM water molecules are colored by atom type. The dipoles are shown in green. (c)–(d) Two dimensional PMF from well-tempered metadyanmics simulations using DFTB3/3OB (barrier 9.8 kcal/mol) and DFTB3/3OBw (barrier 11.5 kcal/mol), respectively. ζ describes the proton transfer and s describes the level of channel hydration (see text).

Figure 5

Study of a proton transfer process inspired by the D-channel in cytochrome c oxidase. (a) Representative structures along the optimized string with key distances (in Å) labeled (water molecules are excluded for clarity); (b) PMF for the first step of proton transfer (approximately corresponding to 1→3 in panel a) studied by umbrella sampling; (c) PMFs for the first step of proton transfer studied by the string method with different collective variables and restraining force constants (indicated in the legend in kcal/mol/Ų).

Figure 6

Computed PMFs (red lines) for the dissociation of pNPP²⁻ in solution along the order parameter d_O^lgP based on (a) DFTB3/MM and (b) B3LYP/MM simulations; pNPP²⁻ is treated as QM and water with TIP3P. For DFTB3, the 3OB/OPhyd parameterization is used; for B3LYP, the 6–31+G(d,p) basis set is used. Also shown in black are the results from adiabatic PMF mapping along d_O^lgP using DFTB3/MM and B3LYP/MM simulations. See text for additional details.

Figure 7

Structural properties of pNPP²⁻ from (a–b) B3LYP/MM and (c–d) DFTB3/MM simulations. The dihedral angle is the P-O-C-C angle indicated in panel (e), which illustrates a “parallel” type of configuration; in (f), a “perpendicular” configuration is illustrated. The histograms in the left (right) column are collected from the d_O^lgP =1.9 (3.1) Å window; note the different numbers of configurations from B3LYP/MM and DFTB3/MM windows. The vertical black line and the numerical value in each histogram indicate the result of adiabatic PMF mapping calculation for the corresponding d_O^lgP window.

Figure 8

Properties from B3LYP/MM simulations for pNPP²⁻ dissociation in solution. The left (right) columns are for d_O^lgP =1.9 (3.1) Å, and each window includes 400 MM configurations. (a-b) illustrate distributions of the QM internal energy (i.e., QM/MM electrostatic energy from self-consistent electrostatic embedding calculation minus the classical QM-MM electrostatic interaction energy), (c-d) illustrate distributions of ESP charge on P. The vertical black line and numerical value indicate results from mean-field calculations (e.g., see Eq. 14 for the internal QM energy).

Figure 9

Potential of mean force (PMF) computed via standard umbrella sampling (US), and perturbation of the PMF following the NBB and FEP routes (see Fig. 2) for (a) intra-molecular proton transfer in malondiadehyde in solution (b) p-nitrophenyl phosphate (pNPP²⁻) dissociation in solution. For the NBB and FEP calculations, the “low-level” method is taken to be DFTB3/MM for both reactions, while the “high-level” method is AM1/MM for (a) and B3LYP/MM for (b). For the NBB/FEP perturbations at the endpoints, only two reactant windows and two transition state windows are used, as shown in Fig. 2, leading to “corrected” barrier heights indicated by dashed lines (also see text). For (a), s₀ = −1.0 Å, s₁ = −0.9 Å, s_n−1 = −0.1 Å, s_n=0.0 Å; for (b), s₀ = 1.6 Å, s₁ = 1.7 Å, s_n−1 = 3.1 Å, s_n=3.2 Å.

Figure 10

Corrected free energy barrier height via NBB/FEP using different numbers of data points per window. Convergence plots are shown for (a) intra-molecular proton transfer in malondiadehyde in solution (b) p-nitrophenyl phosphate (pNPP²⁻) dissociation in solution. NBB/FEP at end points (i.e., reactant and transition state windows) only requires adequate overlaps between s₀-s₁ and s_n−1-s_n windows, thus the results are less dependent on the number of data points.

Figure 11

Distribution of energy difference between the low- and high-levels of theory, ΔU^LH, using trajectories from the two levels of theory; also plotted is the integrand of the FEP expression using the low-level trajectory (i.e., exp(–βΔU^LH)P_L(ΔU^LH)). The ΔU^LH data are shifted with a different amount in (a–d) for plotting purposes; see Table 3 for absolute free energy perturbation results. The plots are for: (a) reactant in malondiadehyde proton transfer (b) transition state in malondiadehyde proton transfer (c) reactant in pNPP²⁻ dissociation (d) transition state in pNPP²⁻ dissociation. In (a)-(b), P_AM1(ΔU) has a similar shape from P_DFTB3(ΔU)exp(-βΔU)), thus FEP and NBB calculations yield reliable barrier corrections. In (c)-(d), the ΔU^LH distributions are remarkably broader and data at different levels have almost vanishing overlap, thus it seems challenging for FEP and NBB to obtain reliable barrier corrections, although the Gaussian nature of the ΔU^LH distributions (fit shown explicitly in green) suggests that a linear response model works well (see Table 3).