Quantum mechanical free energy calculations using path integral molecular dynamics

Abstract

Free energy calculations are one of the most powerful tools within modern theoretical chemistry and are often used to make comparisons with experimental results. Existing free energy calculations are typically performed for classical molecular dynamics simulations but there are certain systems where nuclear quantum effects play an integral role. Specifically, systems with light atoms or low temperatures are the most influenced by such nuclear quantum effects and the development of Feynman path integrals [1] has been effective in accurately describing the quantum nature of these nuclei [2–8]. The primary objective of this thesis is the development of a pair of methodologies to calculate free energies utilizing path integral molecular dynamics to account for nuclear quantum effects.

Prior to the development of these free energy methodologies, this thesis presents a communication interface between the OpenMM and MMTK software packages that has been previously published [9]. This interface allows for users of MMTK to take advantage of the performance of OpenMM without major modifications to existing simulation scripts. Notably, the serial OpenMM integrator is shown to provide a 3x performance gain in comparison to a standard MMTK simulation while the GPU implementations of OpenMM provide over a 400x performance gain for larger systems with periodic boundary conditions.

The first path integral free energy methodology of this thesis combines the existing um- brella sampling technique [10,11] with path integral molecular dynamics. This methodology has been previously published and proposes that the umbrella sampling biasing potential only needs to be applied to a single path integral bead [12]. Furthermore, this proposed methodology is successfully benchmarked for a pair of Lennard-Jones dimer systems before being applied to the more difficult water dimer. The free energy profiles obtained from simulation are then used to calculate a free energy difference of -12.90 ± 0.05 kJ/mol for the MB-Pol potential in comparison to the experimental dissociation energy of -13.2 ± 0.12 kJ/mol [13].

The second path integral free energy methodology introduces a constraint within the path integral molecular dynamics simulations as opposed to an umbrella sampling restraint. Specifically, this methodology applies a constraint to an individual path integral bead in a manner that is similar to the concept of thermodynamic integration for classical simulations [14]. Formal estimators for the derivative of the free energy have been developed by Iouchtchenko et al. [15] and the results presented in this thesis analyze the effectiveness of these estimators for molecular dynamics simulations of Lennard-Jones and water dimers. Additionally, a new estimator is developed and the resulting free energy profiles are used to evaluate a free energy difference for the water dimer of -13.03 ± 0.14 kJ/mol, which is within the errors of the experimental dissociation energy [13].

Overall, this thesis provides a theoretical framework to study the free energy of weakly bound systems over a broad range of temperatures. It is important to note that these methodologies were insufficient below 25 K and it remains more practical to use reaction coordinates that are not distances at such temperatures. Nevertheless, the extension and application of these methodologies to more complicated systems remains an area of exciting development.