Path Integral Molecular Dynamics Calculations Research Articles

Diffusion Monte Carlo and PIMD calculations of radial distribution functions using an updated CCSD(T) potential for CH5 +

We report a new permutationally invariant polynomial (PIP) fit to CCSD(T)/aug-cc-pVTZ energies for CH 5 + that provides fast analytical forces, unlike the original PIP fit of 2006. The new fit, which is slightly more precise than the previous one, is used in classical molecular dynamics, diffusion Monte Carlo and path integral molecular dynamics calculations of CH and HH radial distribution functions (RDFs) at temperatures of 10 and 50 K . Corresponding RDFs are also reported for CD 3 H 2 + . The well-known quantum localisation of the H atoms to the diatomic site from previous and present DMC calculations is maintained at 10 K but shows additional delocalisation to the ‘ CH 3 + ’ site at 50 K .

Modeling Liquid Water by Climbing up Jacob's Ladder in Density Functional Theory Facilitated by Using Deep Neural Network Potentials.

Within the framework of Kohn-Sham density functional theory (DFT), the ability to provide good predictions of water properties by employing a strongly constrained and appropriately normed (SCAN) functional has been extensively demonstrated in recent years. Here, we further advance the modeling of water by building a more accurate model on the fourth rung of Jacob's ladder with the hybrid functional, SCAN0. In particular, we carry out both classical and Feynman path-integral molecular dynamics calculations of water with the SCAN0 functional and the isobaric-isothermal ensemble. To generate the equilibrated structure of water, a deep neural network potential is trained from the atomic potential energy surface based on ab initio data obtained from SCAN0 DFT calculations. For the electronic properties of water, a separate deep neural network potential is trained by using the Deep Wannier method based on the maximally localized Wannier functions of the equilibrated trajectory at the SCAN0 level. The structural, dynamic, and electric properties of water were analyzed. The hydrogen-bond structures, density, infrared spectra, diffusion coefficients, and dielectric constants of water, in the electronic ground state, are computed by using a large simulation box and long simulation time. For the properties involving electronic excitations, we apply the GW approximation within many-body perturbation theory to calculate the quasiparticle density of states and bandgap of water. Compared to the SCAN functional, mixing exact exchange mitigates the self-interaction error in the meta-generalized-gradient approximation and further softens liquid water toward the experimental direction. For most of the water properties, the SCAN0 functional shows a systematic improvement over the SCAN functional. However, some important discrepancies remain. The H-bond network predicted by the SCAN0 functional is still slightly overstructured compared to the experimental results.

Thermodynamics of the insulator-metal transition in dense liquid deuterium

Recent dynamic compression experiments [M. D. Knudson et al., Science 348, 1455 (2015); P. M. Celliers et al., Science 361, 677 (2018)] have observed the insulator-metal transition in dense liquid deuterium, but with an approximately 95-GPa difference in the quoted pressures for the transition at comparable estimated temperatures. It was claimed in the latter of these two papers that a very large latent heat effect on the temperature was overlooked in the first, requiring correction of those temperatures downward by a factor of 2, thereby putting both experiments on the same theoretical phase boundary and reconciling the pressure discrepancy. We have performed extensive path-integral molecular dynamics calculations with density functional theory to directly calculate the isentropic temperature drop due to latent heat in the insulator-metal transition for dense liquid deuterium and show that this large temperature drop is not consistent with the underlying thermodynamics.

Path integral molecular dynamic simulation of flexible molecular systems in their ground state: Application to the water dimer.

We extend the Langevin equation Path Integral Ground State (LePIGS), a ground state quantum molecular dynamics method, to simulate flexible molecular systems and calculate both energetic and structural properties. We test the approach with the H2O and D2O monomers and dimers. We systematically optimize all simulation parameters and use a unity trial wavefunction. We report ground state energies, dissociation energies, and structural properties using three different water models, two of which are empirically based, q-TIP4P/F and q-SPC/Fw, and one which is ab initio, MB-pol. We demonstrate that our energies calculated from LePIGS can be merged seamlessly with low temperature path integral molecular dynamics calculations and note the similarities between the two methods. We also benchmark our energies against previous diffusion Monte Carlo calculations using the same potentials and compare to experimental results. We further demonstrate that accurate vibrational energies of the H2O and D2O monomer can be calculated from imaginary time correlation functions generated from the LePIGS simulations using solely the unity trial wavefunction.

Exploring Systematic Discrepancies in DFT Calculations of Chlorine Nuclear Quadrupole Couplings.

Previous studies have revealed significant discrepancies between density functional theory (DFT)-calculated and experimental nuclear quadrupolar coupling constants (CQ) for chlorine atoms, particularly in ionic solids. Various aspects of the computations are systematically investigated here, including the choice of the DFT functional, basis set convergence, and geometry optimization protocol. The effects of fast (fs) time-scale dynamics are probed using molecular dynamics (MD) and nuclear quantum effects (NQEs) are considered using path-integral MD calculations. It is shown that the functional choice is the most important factor related to improving the accuracy of the quadrupolar coupling calculations, and that functionals beyond the generalized gradient approximation (GGA) level, such as hybrid and meta-GGA functionals, are required for good correlations with experiment. The influence of molecular dynamics and NQEs is less important than the functional choice in the studied systems. A method which involves scaling the calculated quadrupolar coupling constant is proposed here; its application leads to good agreement with experimental data.

Equilibrium fractionation of H and O isotopes in water from path integral molecular dynamics

The equilibrium fractionation factor between two phases is of importance for the understanding of many planetary and environmental processes. Although thermodynamic equilibrium can be achieved between minerals at high temperature, many natural processes involve reactions between liquids or aqueous solutions and solids. For crystals, the fractionation factor α can be theoretically determined using a statistical thermodynamic approach based on the vibrational properties of the phases. These calculations are mostly performed in the harmonic approximation, using empirical or ab-initio force fields. In the case of aperiodic and dynamic systems such as liquids or solutions, similar calculations can be done using finite-size molecular clusters or snapshots obtained from molecular dynamics (MD) runs. It is however difficult to assess the effect of these approximate models on the isotopic fractionation properties. In this work we present a systematic study of the calculation of the D/H and 18O/16O equilibrium fractionation factors in water for the liquid/vapour and ice/vapour phases using several levels of theory within the simulations. Namely, we use a thermodynamic integration approach based on Path Integral MD calculations (PIMD) and an empirical potential model of water. Compared with standard MD, PIMD takes into account quantum effects in the thermodynamic modeling of systems and the exact fractionation factor for a given potential can be obtained. We compare these exact results with those of modeling strategies usually used, which involve the mapping of the quantum system on its harmonic counterpart. The results show the importance of including configurational disorder for the estimation of isotope fractionation in liquid phases. In addition, the convergence of the fractionation factor as a function of parameters such as the size of the simulated system and multiple isotope substitution is analyzed, showing that isotope fractionation is essentially a local effect in the investigated system.

Strong Isotope Effect in Phase II of Dense Solid Hydrogen and Deuterium

Quantum nuclear zero-point motions in solid H(2) and D(2) under pressure are investigated at 80 K up to 160 GPa by first-principles path-integral molecular dynamics calculations. Molecular orientations are well defined in phase II of D(2), while solid H(2) exhibits large and very asymmetric angular quantum fluctuations in this phase, with possible rotation in the (bc) plane, making it difficult to associate a well-identified single classical structure. The mechanism for the transition to phase III is also described. Existing structural data support this microscopic interpretation.

Improving the convergence of closed and open path integral molecular dynamics via higher order Trotter factorization schemes

Higher order factorization schemes are developed for path integral molecular dynamics in order to improve the convergence of estimators for physical observables as a function of the Trotter number. The methods are based on the Takahashi-Imada and Susuki decompositions of the Boltzmann operator. The methods introduced improve the averages of the estimators by using the classical forces needed to carry out the dynamics to construct a posteriori weighting factors for standard path integral molecular dynamics. The new approaches are straightforward to implement in existing path integral codes and carry no significant overhead. The Suzuki higher order factorization was also used to improve the end-to-end distance estimator in open path integral molecular dynamics. The new schemes are tested in various model systems, including an ab initio path integral molecular dynamics calculation on the hydrogen molecule and a quantum water model. The proposed algorithms have potential utility for reducing the cost of path integral molecular dynamics calculations of bulk systems.

Path Integral Computation of Quantum Free Energy Differences Due to Alchemical Transformations Involving Mass and Potential.

Thermodynamic integration, perturbation theory, and λ-dynamics methods were applied to path integral molecular dynamics calculations to investigate free energy differences due to "alchemical" transformations. Several estimators were formulated to compute free energy differences in solvable model systems undergoing changes in mass and/or potential. Linear and nonlinear alchemical interpolations were used for the thermodynamic integration. We find improved convergence for the virial estimators, as well as for the thermodynamic integration over nonlinear interpolation paths. Numerical results for the perturbative treatment of changes in mass and electric field strength in model systems are presented. We used thermodynamic integration in ab initio path integral molecular dynamics to compute the quantum free energy difference of the isotope transformation in the Zundel cation. The performance of different free energy methods is discussed.

Proton transfer dynamics in the propionic acid dimer from path integral molecular dynamics calculations

The double proton transfer process in the cyclic dimer of propionic acid in the gas phase was studied using a path integral molecular dynamics method. Structures, energies and proton trajectories were determined. Very large amplitude motions of the skeleton of a propionic acid molecule were observed during the simulations, and almost free rotation of the C(2)H(5) group around the C(α)-C bond. A double-well symmetric potential with a very small energy barrier was determined from the free energy profile for the proton motions. Infrared spectra for different isotopomers were calculated, and comparative vibrational analysis was performed. The vibrational results from CPMD appear to be in qualitative agreement with the experimental ones.

Path Integral Molecular Dynamics Study of Small H2 Clusters in the Large Cage of Structure II Clathrate Hydrate: Temperature Dependence of Quantum Spatial Distributions

We report a path integral molecular dynamics (PIMD) study of the temperature dependence of the spatial distribution of two and four H2 molecules inside the large cage of the structure II clathrate hydrate. The PIMD calculations were performed at five temperatures ranging from 25 to 200 K. Their results were combined with those from an earlier diffusion Monte Carlo (DMC) study of this system at T = 0 K [Sebastianelli, F.; Xu, M.; Bačić, Z. J. Chem. Phys. 2008, 129, 244706]. The spatial distribution of the confined H2 molecules at each of the temperatures considered was characterized with the help of several one-dimensional (1D) and three-dimensional (3D) distribution functions of suitably chosen intermolecular coordinates, generated by the PIMD and DMC calculations. The 1D distribution that proved to be the most strongly temperature dependent, and also the most revealing about the structural properties of the system as a function of temperature, involves the H2−cage center−H2 angle. In the case of four caged H2 molecules, this angular distribution provides clear evidence that between 50 and ∼100 K the system undergoes a qualitative change. At 50 K and below, the system is fully localized in the global minimum of the intermolecular potential, corresponding to a tetrahedral configuration of H2 molecules with a unique orientation relative to the cage frame. At temperatures of 75−100 K and higher, nearly degenerate local minima ∼200 cm−1 above the global minimum become accessible and are increasingly sampled by the system. The 3D spatial distributions also show this growing delocalization above 75−100 K. Our findings are in accord with the localization−delocalization transition observed experimentally to occur at 50 K for four D2 molecules in the large cage [Lokshin, K. A.; Zhao, Y.; He, D.; Mao, W. L.; Mao, H. K.; Hemley, R. J.; Lobanov, M. V.; Greenblatt, M. Phys. Rev. Lett. 2004, 93, 125503].

Enol Tautomers of Watson−Crick Base Pair Models Are Metastable Because of Nuclear Quantum Effects

Intermolecular enol tautomers of Watson-Crick base pairs could emerge spontaneously via interbase double proton transfer. It has been hypothesized that their formation could be facilitated by thermal fluctuations and proton tunneling, and possibly be relevant to DNA damage. Theoretical and computational studies, assuming classical nuclei, have confirmed the dynamic stability of these rare tautomers. However, by accounting for nuclear quantum effects explicitly through Car-Parrinello path integral molecular dynamics calculations, we find the tautomeric enol form to be dynamically metastable, with lifetimes too insignificant to be implicated in DNA damage.

Vibron hopping and bond anharmonicity in hot dense hydrogen

The Raman-active vibron of dense hydrogen has been shown to exhibit unexpected changes as a function of pressure and temperature to above 100 GPa. To understand these results we have performed supercell-based calculations using Van Kranendonk theory taking into account the renormalization of the hopping parameter by the lattice vibrations. We find that the major temperature dependence at this level of theory comes from the differences in populations of rotational states. The theory provides a fair description of the experimental results up to 70 GPa. We examine in detail a number of assumptions made in the application of the Van Kranendonk model to hydrogen as a function of pressure and temperature. We also present results of hybrid path integral molecular dynamics calculations in the fluid state at a low pressure (7 GPa) near the melting temperature. An amorphous-solid model of the fluid predicts that the Raman vibron frequencies change little upon melting, in agreement with experiment. The Van Kranendonk theory with fixed rotational identities of the molecules tends to predict more peaks in the Raman spectrum than are observed experimentally.

Excess Electron in Water at Different Thermodynamic Conditions

A hydrated electron in water at different densities and temperatures is studied via a set of density functional based molecular dynamics simulations, showing that a localization of an excess electron is still present even at very low densities. Space variations of the molecular dipole moments are analyzed, proposing a simple algorithm to identify the region of localization of the wavefunction relative to the solvated electron in terms of orientation of the H2O molecular dipole moments. Finally, the effects of the self-interaction corrections on the optical absorption spectra are analyzed and compared with both available experimental data and path integral molecular dynamics calculations, showing that a weighted subtraction of the self-interaction yields a systematic improvement in the position of the absorption peak.

DEVELOPMENT OF A THREE-DIMENSIONAL AB INITIO POTENTIAL ENERGY SURFACE FOR THE He–Cl2(X) SYSTEM AND ITS APPLICATION TO SOLVATION STRUCTURES IN THE HenCl2 CLUSTERS

High-quality ab initio electronic structure calculations for the van der Waals interaction of He with Cl 2 in the electronic ground state have been carried out to develop a new three-dimensional potential energy surface for this system. The calculations were performed at the single and double excitation coupled-cluster level of theory with non-iterative perturbational treatment of triple excitations [CCSD(T)] with a very large basis set including an additional set of bond functions. The analytical potential surface developed were then used in the path-integral molecular dynamics calculations for the He n Cl 2 cluster, where quantum solvation structures of helium atoms in clusters were investigated. It has been found that the helium solvation structures are quite different between the electronic ground state and the electronically excited B 3Π state.

Quantum path-integral molecular dynamics calculations of the dipole-bound state of the water dimer anion

The equilibrium structure of the negatively charged water dimer (H 2O) 2 − has been studied using the path-integral molecular dynamics simulation. All the atomic motions as well as the excess electron were treated quantum mechanically, employing a semi-empirical model combining a water–water interatomic potential with an electron–water pseudopotential. It is demonstrated that the molecular structure of (H 2O) 2 − is more flexible than that of (H 2O) 2; both the donor switching and donor–acceptor interchange can more effectively occur in (H 2O) 2 − than in (H 2O) 2. We conclude that this floppy character is a result of the breakdown of the adiabatic Born–Oppenheimer picture.