Quantum mechanical free energy calculations using path integral molecular dynamics

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

Hybrid Path Integral Monte Carlo/Molecular Dynamics Approach for the Simulation of Rigid Rotating and Translating Molecules

The path integral formulation of quantum statistical mechanics [1] is widely used for formulating the partition function and diverse thermodynamic properties. The basic objective of this thesis is to present the benchmark calculations and develop new methods to study the quantum rotating and translating system. Prior to the development of a new method, we perform benchmark calculations on newly implemented Path Integral Monte Carlo code in the Molecular Modelling Toolkit software package [2] against exact diagonalization calculation on Hydrogen Fluoride in an electric field. An energy convergence study proved the validity of the new integrated Path Integral Monte Carlo code within statistical error. Then, we benchmark our PIMC code against the direct Classical Monte Carlo (dCMC) calculations to study the interacting quantum rotors at finite temperature. Then, the application of the direct Classical Monte Carlo code for a 1D water chain system at different temperature and various lattice spacing is performed. Two phase transitions and three specific phase regions have been observed from the calculations. Then, we combine our well-established Path Integral Molecular Dynamics code [3, 4] with the Path Integral Monte Carlo code to develop a novel hybrid Path Integral Molecular Dynamics/Path Integral Monte Carlo program for quantum system simulation including rotation and translation. This proposed methodology is successfully benchmarked for CO 2 He system at a large tau region. The second Path Integral Monte Carlo methodology introduces a cluster-update algorithm applying on quantum rotational system. This methodology has been previously applied to lattice spin system [5, 6] and hard sphere continuum system [7]. This is the first time to discuss and develop a cluster-update method on quantum rotors. A primary cluster-update Path Integral Monte Carlo algorithm and related test is presented. The preliminary results for the test of Hydrogen Fluoride molecules with dipole-dipole interaction shows an existing problem within our code, and the problem is fully discussed in this thesis for future study. Overall, this thesis focuses on the development of new methodologies to study quantum molecules. Noted that all methodologies require more tests and benchmark calculations to prove its validity. In spite of that, this discussion of a novel methodology and application will lead to an exciting further study