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

2019-10-25T20:17:28Z (GMT) by Braden Kelly William Smith
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.