Alchemical Hydration Free Energy Calculations Using Molecular Dynamics with Explicit Polarization and Induced Polarity Decoupling: An OTFP Approach

We present an algorithm to calculate hydration free energies in explicit solvent that incorporates polarization of the solute molecule in conjunction with the use of a classical fixed--charge force field. The goal is to improve the accuracy over the alternative approach of developing a polarizable force field with adjustable parameters. We incorporate polarization by implementing on--the--fly periodic updating of the solute's partial charges during a standard molecular dynamics (MD) alchemical change simulation by the use of mixed QM/MM calculations. We decouple the polarizing solvent's electric field along with the normal MD solute Coulomb decoupling to calculate the free energy difference between an unpolarized solute in vacuum and a fully polarized solute in solution. This approach is in contrast to the common approach of GAFF, which calculates the difference between a solute in vacuum that is over--polarized by the use of fixed charges calculated using HF/6-31G*, and correspondingly under--polarized by the same partial charge set in the solution phase. We apply our methodology to a test set of 31 molecules, ranging from small polar to large drug--like molecules. We find that results using our method with Minimum Basis Iterative Stockholder (MBIS) charges and using RESP charges with B3LYP/cc-pVTZ are superior to results calculated using the current ``gold standard" AM1--BCC method. We show results using MBIS partial charges using B3LYP/cc-pVTZ and MP2/cc-pVTZ, RESP partial charges using B3LYP/cc-pVTZ and HF/6-31G*, and AM1-BCC partial charges. Our method using MBIS in conjunction with MP2/cc-pVTZ yields an AAD that is 2.91 kJ$\cdot$mol$^{-1}$ (0.70 kcal$\cdot$mol$^{-1}$) lower than that of AM1--BCC for our test set. AM1-BCC was within experimental uncertainty on 13 \% of the data, while our method using MP2 was within experimental uncertainty on 43 \% of the data. We conjecture that results can be further improved by using Lennard--Jones and torsional parameters that are fitted to the MBIS charge method and that using RESP with our method can be improved by using a higher level of theory than B3LYP, for instance MP2 or $\omega$B97X-D.