QM/CG-MM: Systematic Embedding of Quantum Mechanical Systems in a Coarse-Grained Environment with Accurate Electrostatics

Abstract

Quantum Mechanics/Molecular Mechanics (QM/MM) can describe chemical reactions in molecular dynamics (MD) simulations at a much lower cost than ab initio MD. Still, it is prohibitively expensive for many systems of interest because such systems usually require long simulations for sufficient statistical sampling. Additional MM degrees of freedom are often slow and numerous but secondary in interest. Coarse-graining (CG) is well-known to be able to speed up sampling through both reduction in simulation cost and the ability to accelerate the dynamics. Therefore, embedding a QM system in a CG environment can be a promising way of expediting sampling without compromising the information about the QM subsystem. Sinitskiy and Voth first proposed the theory of Quantum Mechanics/Coarse-grained Molecular Mechanics (QM/CG-MM) with a bottom-up CG mapping. Mironenko and Voth subsequently introduced the DFT-QM/CG-MM formalism to couple a Density Functional Theory (DFT) treated QM system and to an apolar environment. Here, we present a more complete theory that addresses MM environments with significant polarity by explicitly accounting for the electrostatic coupling. We demonstrate our QM/CG-MM method with a chloride-methyl chloride S_N2 reaction system in acetone, which is sensitive to solvent polarity. The method accurately recapitulates the potential of mean force for the substitution reaction, and the reaction barrier from the best model agrees with the atomistic simulations within sampling error. These models also have generalizability. In two other reactive systems that they have not been trained on, the QM/CG-MM model still achieves the same level of agreement with the atomistic QM/MM models. Finally, we show that in these examples the speed-up in the sampling is proportional to the acceleration of the rotational dynamics of the solvent in the CG system.

Figure 1:. QM/CG-MM umbrella sampling of the chloride-methyl chloride system.

(A) The reaction PMF or single point energy of the reactive system. (B) The schematic plot of the chloride-methyl chloride system and the definition of the reaction coordinate (also referred to as CV). Error bars are from a block analysis of 5 blocks, shown as 95% CI.

Figure 2:. One-site acetone model.

Left: CG potential between beads. Right: RDF between beads

Figure 3:. Mapping of the acetone molecule.

(A) Five 2D tomography slices (snapshots) of the electrostatic potential in the QM region, calculated from one sample of the 100 snapshots in atomistic and two different CG resolutions. Blue means negative value, and red represents positive. (B) Coupling parameters from fitting. (C) The two-site CG mapping scheme. Everything else but the oxygen atom was mapped into a single bead. (D) The three-site CG mapping scheme. The CH₃─C─CH₃ part was split into two equal beads, breaking the carbonyl carbon atom into two beads.

Figure 4:

CG potentials and radial distribution functions, respectively. The dashed lines showed atomistic simulation RDF (reference), and the solid line is from CG simulation. (A–B) Two-site model. (C–D) Three-site model.

Figure 5:

Reaction PMF for the bromide-methyl chloride system

Figure 6:. Kemp decarboxylation reaction in acetone.

(A) The 2-D reaction PMF as a function of the two CVs. (B) The 1-D reaction PMF integrated from the 2-D PMF, following the minimum free energy path. (C) The reaction scheme and the definition of the two CVs.

Figure 7:. The convergence of QM/CG-MM.

(A) The PMF from truncated trajectories. (B) The dipole auto-correlation function for different acetone models in a pure solvent simulation. (C) Dielectric relaxation time (in ps) from a linear fit of ln $C(t)$ over $t$. The shown uncertainty is 95 % CI.