BAR-Based Multi-dimensional Nonequilibrium Pulling for Indirect Construction of QM/MM Free Energy Landscape: From Semi-Empirical to Ab Initio

The indirect method for the construction of Quantum mechanics (QM)/ molecular mechanics (MM) free energy landscapes provides a cheaper alternative of free energy simulations at QM level. The indirect method features a direct calculation of the free energy profile at a relatively cheap low-level Hamiltonian and a low-level to high-level correction. In the thermodynamic cycle, the direct low-level calculation along the physically meaningful reaction coordinate is corrected via the alchemical method, which is often achieved with perturbation-based techniques. Often the indirect method can lead to about an order of magnitude speedup in free energy simulation. In our previous work, a multi-dimensional nonequilibrium pulling framework is proposed for the indirect construction of QM/MM free energy landscapes. The method relies on bidirectional nonequilibrium pulling and bidirectional reweighting with the statistically optimal estimator. In the previous work, we focus on obtaining semi-empirical QM (SQM) results indirectly from direct MM simulations and MM<->SQM corrections. In this work, we apply this method to obtain results under ab initio QM Hamiltonians by combining direct SQM results and SQM<->QM corrections. The indirect nonequilibrium scheme is tested on a dihedral flipping case and a series of SQM and QM Hamiltonians are benchmarked. It is observed that PM6 achieves the best performance among the low-level Hamiltonians, while AM1 and MNDO perform less well. Therefore, we recommend using PM6 as the low-level theory in the indirect free energy simulation. The comparison between the indirect results from different SQM Hamiltonians could also provide some hints on the development of charge models. As AM1 can be corrected with the bond charge correction (BCC) to provide a cheap and accurate charge model, which is able to accurately reproduce the electrostatic potential (ESP) at HF level, PM6 would be able to do the same thing. Considering its higher similarity to the high-level Hamiltonians, the PM6-BCC model could be more accurate than the existing AM1-BCC model. Another central result in the current work is a basic protocol of choosing the strength of restraints and an appropriate time step in nonequilibrium free energy simulation at the stiff spring limit. We provide theoretical derivations to emphasize the importance of using a sufficiently large force constant and choosing an appropriate time step. It is worth noting that a general rule of thumb for choosing the time step, according to our derivation, is that a time step of 1 fs or smaller should be used, as long as the stiff spring approximation is employed, even in simulations with constraints on bonds involving hydrogen atoms.