Tunneling Splittings in Water Clusters from Path Integral Molecular Dynamics.

We present the aspects of the path integral molecular dynamics (PIMD) method relative to (a) its theoretical fundamental principles, (b) its applicability to model quantum systems and (c) its implementation as a simulation tool. The PIMD method is based on the discretized path integral representation of quantum mechanics. In this representation, a quantum particle is isomorphic to a closed polymer chain. The problem of the indistinguishability of quantum particles is tackled with a non-local exchange potential. When the exact density matrix of the quantum particles is used, the exchange potential is exact. We use a high temperature approximation to the density matrix leading to an approximate form of the exchange potential. This quantum molecular dynamics method allows the simulation of collections of quantum particles at finite temperature. Our algorithm can be made to scale linearly with the number of quantum states on which the density matrix is projected. Therefore, it can be optimized to run efficiently on parallel computers. We apply the PIMD method to the electron plasma in 3-dimension. The kinetic and potential energies are calculated and compared with results for similar systems simulated with a variational Monte Carlo method. Both results show good agreements with each other at all the densities studied. The method is then use to model the thermodynamic behavior of a simple alkali metal. In these simulations, ions and valence electrons are treated as classical and quantum particles, respectively. The simple metal undergoes a phase transformation upon heating. Furthermore, to demonstrate the richness of behaviors that can be studied with the PIMD method, we also report on the metal to insulator transition in a hydrogenoid lattice. Finally, in previous studies of alkali metals, electrons interacted with ions via local pseudo-potentials, an extension of the method to modeling electrons in non-local pseudo-potentials is also presented with applications.

A path integral centroid molecular dynamics method for Bose and Fermi statistics

Tunneling splittings observed in molecular rovibrational spectra are significant evidence for tunneling motion of hydrogen nuclei in water clusters. Accurate calculations of the splitting sizes from first principles require a combination of high-quality inter-atomic interactions and rigorous methods to treat the nuclei with quantum mechanics. Many theoretical efforts have been made in recent decades. This Perspective focuses on two path-integral based tunneling splitting methods whose computational cost scales well with the system size, namely, the ring-polymer instanton method and the path-integral molecular dynamics (PIMD) method. From a simple derivation, we show that the former is a semiclassical approximation to the latter, despite that the two methods are derived very differently. Currently, the PIMD method is considered to be an ideal route to rigorously compute the ground-state tunneling splitting, while the instanton method sacrifices some accuracy for a significantly smaller computational cost. An application scenario of such a quantitatively rigorous calculation is to test and calibrate the potential energy surfaces of molecular systems by spectroscopic accuracy. Recent progress in water clusters is reviewed, and the current challenges are discussed.

The impact of quantum nuclear effects on hydrogen (H-) bond strength has been inferred in earlier work from bond lengths obtained from path integral molecular dynamics (PIMD) simulations. To obtain a direct quantitative assessment of such effects, we use constrained-centroid PIMD simulations to calculate the free energy changes upon breaking the H-bonds in dimers of HF and water. Comparing ab initio simulations performed using PIMD and classical nucleus molecular dynamics (MD), we find smaller dissociation free energies with the PIMD method. Specifically, at 50 K, the H-bond in (HF)(2) is about 30% weaker when quantum nuclear effects are included, while that in (H(2)O)(2) is about 15% weaker. In a complementary set of simulations, we compare unconstrained PIMD and classical nucleus MD simulations to assess the influence of quantum nuclei on the structures of these systems. We find increased heavy atom distances, indicating weakening of the H-bond consistent with that observed by direct calculation of the free energies of dissociation.

In a previous joint experimental and theoretical study of the barrierless chemical reaction C(1D) + H2 at low temperatures (300-50 K) [K. M. Hickson, J.-C. Loison, H. Guo and Y. V. Suleimanov, J. Phys. Chem. Lett., 2015, 6, 4194], excellent agreement was found between experimental thermal rate constants and theoretical estimates based on ring polymer molecular dynamics (RPMD) over the two lowest singlet potential energy surfaces (PESs). Here, we extend this work to one of its isotopologues, C(1D) + D2, over the same temperature range. Experimental and RPMD results are in good agreement when contributions from both PESs to this chemical reaction are included in the RPMD simulations. The deviation between experiment and the RPMD calculations does not exceed 25% and both results exhibit a slight negative temperature dependence. The first excited 1A'' PES plays a more important role than the ground 1A' PES as the temperature is decreased, similar to our previous studies of the C(1D) + H2 reaction but with a more pronounced effect. The small differences in temperature dependence between the earlier and present experimental studies of C(1D) + H2/D2 reactions are discussed in terms of the use of non-equilibrium populations of ortho/para-H2/D2. We argue that RPMD provides a very convenient and reliable tool to study low-temperature chemical reactions.

We present the results obtained from spectroscopic investigations and quantum chemical calculations of the interaction of anisole (methoxybenzene) with small water clusters. The experiments have been carried out using resonant two-photon ionization (R2PI) and IR-UV double-resonance vibrational spectroscopy (IR/R2PI) in the region of the OH stretches. Apart from the vibrational spectra of the water moiety in the clusters, their intermolecular vibrations in the electronically excited S1 state are identified by IR/R2PI hole burning spectroscopy and assigned according to the vibrations calculated for the S1 state and compared with the vibrations calculated for the S0 state. The calculations for the S0 state were carried out at the second order Møller-Plesset level of theory using both the 6-31+G* and aug-cc-pVDZ basis sets and for the S1 state at the configuration interaction singles (CIS) level with the 6-31+G* basis set. In the electronic ground state (S0), the interaction of a water monomer to anisole is mediated through its oxygen atom, and that of a water dimer both through the oxygen atom (σ type of interaction) and the arene ring (π type of interaction). Thus in contrast to the interaction of fluorinated benzenes with water clusters, wherein a conformational transition from an in-plane σ to a on-top π bonding emerges starting with a water trimer, this conformational transition appears in case of anisole already with a water dimer. In the excited state (S1) of the investigated systems, there is a pronounced weakening of the interaction of the water cluster with the aromatic chromophore, which is also responsible for the blue shift of the electronic transitions. Consequently, the structures of the complexes of anisole with a water monomer or dimer are very different in both states. The weakening results from a diminished electron density of the oxygen atom and of the π system of anisole in the excited state. The calculated binding energies of the ground-state conformers indicate that these small water clusters are bound more strongly to anisole than to other π systems like benzene, toluene, fluorobenzene, and p-difluorobenzene. The many-body decomposition of the binding energy (S0) reveals that a progressive increase in the size of the water clusters results in a weakening of the σ O⋯H interaction and a concurrent strengthening of the π-H interaction. The complex containing a cyclic water trimer exhibits in the excited state also a π-type H-bonding interaction, but its stability emerges from a delocalization of the electron density from the water trimer to the anisole oxygen. Excepting the water dimer complexes, there is a good agreement between the calculated OH vibrational frequencies and the experimental IR spectra recorded for the ground state. The surprisingly good agreement of the calculated and the experimentally observed intermolecular modes in the excited state of these complexes provides convincing evidence that the experimental spectra emerge from structures similar to those predicted at the CIS/6-31+G* level.

To fast and accurately compute rate coefficients of the H/D + CH4 → H2/HD + CH3reactions, we propose a segmented strategy for fitting suitable potential energy surface (PES), on which ring-polymer molecular dynamics (RPMD) simulations are performed. On the basis of recently developed permutation invariant polynomial neural-network approach [J. Li et al., J. Chem. Phys. 142, 204302 (2015)], PESs in local configuration spaces are constructed. In this strategy, global PES is divided into three parts, including asymptotic, intermediate, and interaction parts, along the reaction coordinate. Since less fitting parameters are involved in the local PESs, the computational efficiency for operating the PES routine is largely enhanced by a factor of ∼20, comparing with that for global PES. On interaction part, the RPMD computational time for the transmission coefficient can be further efficiently reduced by cutting off the redundant part of the child trajectories. For H + CH4, good agreements among the present RPMD rates and those from previous simulations as well as experimental results are found. For D + CH4, on the other hand, qualitative agreement between present RPMD and experimental results is predicted.

The path-integral molecular dynamics and centroid molecular dynamics methods have been applied to investigate the behavior of liquid water at ambient conditions starting from a recently developed simple point charge/flexible (SPC/Fw) model. Several quantum structural, thermodynamic, and dynamical properties have been computed and compared to the corresponding classical values, as well as to the available experimental data. The path-integral molecular dynamics simulations show that the inclusion of quantum effects results in a less structured liquid with a reduced amount of hydrogen bonding in comparison to its classical analog. The nuclear quantization also leads to a smaller dielectric constant and a larger diffusion coefficient relative to the corresponding classical values. Collective and single molecule time correlation functions show a faster decay than their classical counterparts. Good agreement with the experimental measurements in the low-frequency region is obtained for the quantum infrared spectrum, which also shows a higher intensity and a redshift relative to its classical analog. A modification of the original parametrization of the SPC/Fw model is suggested and tested in order to construct an accurate quantum model, called q-SPC/Fw, for liquid water. The quantum results for several thermodynamic and dynamical properties computed with the new model are shown to be in a significantly better agreement with the experimental data. Finally, a force-matching approach was applied to the q-SPC/Fw model to derive an effective quantum force field for liquid water in which the effects due to the nuclear quantization are explicitly distinguished from those due to the underlying molecular interactions. Thermodynamic and dynamical properties computed using standard classical simulations with this effective quantum potential are found in excellent agreement with those obtained from significantly more computationally demanding full centroid molecular dynamics simulations. The present results suggest that the inclusion of nuclear quantum effects into an empirical model for water enhances the ability of such model to faithfully represent experimental data, presumably through an increased ability of the model itself to capture realistic physical effects.

Efficient numerical algorithms are developed for use with two finite temperature semiclassical approximations to quantum dynamics both of which require trajectories generated on potentials of mean force derived from the path integral expression for the density matrix. The numerical algorithms are formed from the combination of a classical adiabatic relation similar to that used in the Car–Parrinello method and an efficient path integral molecular dynamics scheme. Results on model, an anharmonic oscillator and a realistic, fluid para-hydrogen, problem indicate that semiclassical dynamics can be obtained for virtually the same computational cost as structure and thermodynamics.

Ring polymer molecular dynamics (RPMD) calculations for the C(1D)+H2 reaction are performed on the Zhang-Ma-Bian ab initio potential energy surfaces (PESs) recently constructed by our group, which are unique in very good descriptions of the regions around conical intersections and of van der Waals (vdW) interactions. The calculated reaction thermal rate coefficients are in very good agreement with the latest experimental results. The rate coefficients obtained from the ground ã1A′ ZMB-a PES are much larger than those from the previous RKHS PES, which can be attributed to that the vdW saddles on our PESs have very different dynamical effects from the vdW wells on the previous PESs, indicating that the RPMD approach is able to include dynamical effects of the topological structures caused by vdW interactions. The importance of the excited b̃1A″ ZMB-b PES and quantum effects in the title reaction is also underscored.

The maximum entropy analytic continuation (MEAC) and ring polymer molecular dynamics (RPMD) methods provide complementary approaches to the calculation of real time quantum correlation functions. RPMD becomes exact in the high temperature limit, where the thermal time betavariant Planck's over 2pi tends to zero and the ring polymer collapses to a single classical bead. MEAC becomes most reliable at low temperatures, where betavariant Planck's over 2pi exceeds the correlation time of interest and the numerical imaginary time correlation function contains essentially all of the information that is needed to recover the real time dynamics. We show here that this situation can be exploited by combining the two methods to give an improved approximation that is better than either of its parts. In particular, the MEAC method provides an ideal way to impose exact moment (or sum rule) constraints on a prior RPMD spectrum. The resulting scheme is shown to provide a practical solution to the "nonlinear operator problem" of RPMD, and to give good agreement with recent exact results for the short-time velocity autocorrelation function of liquid parahydrogen. Moreover these improvements are obtained with little extra effort, because the imaginary time correlation function that is used in the MEAC procedure can be computed at the same time as the RPMD approximation to the real time correlation function. However, there are still some problems involving long-time dynamics for which the RPMD+MEAC combination is inadequate, as we illustrate with an example application to the collective density fluctuations in liquid orthodeuterium.

The rate coefficients of the H + H2O2 → H2 + HO2 reaction are calculated using the ring polymer molecular dynamics (RPMD), quasi-classical trajectory (QCT), and canonical variational transition state theory (CVT) with small curvature tunneling (SCT) correction, in conjunction with the recently constructed fundamental invariant-neural network (FI-NN) potential energy surface (PES) [X. Lu et al., Phys. Chem. Chem. Phys. 20, 23095 (2018)]. In RPMD calculations, 32, 16, and 8 beads are used for computing the rate coefficients at 200 K ≤ T ≤ 400 K, 500 K ≤ T ≤ 700 K, and 700 K < T ≤ 1000 K, respectively. Given that the previous experimental rate coefficients vary widely, in particular, at low temperatures, the present RPMD rate coefficients agree well with most of the experimental results. In addition, comparing with some experimental values, the present QCT and CVT/SCT calculations on the FI-NN PES also predict accurate results at some temperatures. These results strongly support the accuracy of the present dynamics calculations as well as the full-dimensional FI-NN PES.

We examine the short-time accuracy of a class of approximate quantum dynamical techniques that includes the centroid molecular dynamics (CMD) and ring polymer molecular dynamics (RPMD) methods. Both of these methods are based on the path integral molecular dynamics (PIMD) technique for calculating the exact static equilibrium properties of quantum mechanical systems. For Kubo-transformed real-time correlation functions involving operators that are linear functions of positions or momenta, the RPMD and (adiabatic) CMD approximations differ only in the choice of the artificial mass matrix of the system of ring polymer beads that is employed in PIMD. The obvious ansatz for a general method of this type is therefore to regard the elements of the PIMD (or Parrinello-Rahman) mass matrix as an adjustable set of parameters that can be chosen to improve the accuracy of the resulting approximation. We show here that this ansatz leads uniquely to the RPMD approximation when the criterion that is used to select the mass matrix is the short-time accuracy of the Kubo-transformed correlation function. In particular, we show that the leading error in the RPMD position autocorrelation function is O(t(8)) and the error in the velocity autocorrelation function is O(t(6)), for a general anharmonic potential. The corresponding errors in the CMD approximation are O(t(6)) and O(t(4)), respectively.

In order to predict the properties of materials from first principles simulations, both the electronic structures and the movement of nuclei should be accurately described. Till now, the electronic structure calculation methods have been relatively well developed, like the density-functional theory, the quantum chemical methods, and the diffusion quantum monte carlo method, etc. By using them, we can get the electronic structures of materials with different levels of accuracy. But the movements of nuclei are still often treated at a clasical level. More and more experimental and theoretical works have found that, the nuclear quantum effects (NQEs) are playing an important role in some phase transition porcesses. Melting line of hydrogen at high pressures were reported through coexistance simulations based on ab initio PIMD methods. We determined the melting temperature as a function of pressure and finded an atomic solid phase from 500 to 800 GPa, which melts at less than 200 K. Beyond this and up to 1200 GPa, a metallic atomic liquid is stable at temperatures as low as 50 K. The quantum motion of the protons is the key to the low melting temperature reported, as simulations with classical nuclei lead to considerably higher melting temperatures of about 300 K across the entire pressure range considered. Then, a combination of state-of-the-art theoretical methods were used to obtain an atomic level picture of classical and quantum ordering of protons in cold high-pressure solid hydrogen. Besides, Using a self-developed combination of the thermodynamic integration and the ab initio PIMD methods, we quantitatively studied the influence of NQEs on the melting of dense lithium at 45 GPa. We find that although the NQEs significantly change the free-energies of the competing solid and liquid phases, the melting temperature is lowered by only about 15 K, with values obtained using both classical and quantum nuclei in close proximity to a new experiment. So, by combining the path-integral molecular dynamics (PIMD) method with techniques like coexistance simulation and free-energy calculation method, the NQEs in many phase transition processes can be well understood.

The rate coefficients of the barrierless O(1D) + CH4 reaction are determined both theoretically and experimentally at 50-296 K. For the calculations, ring polymer molecular dynamics (RPMD) simulations are performed on the basis of a new neural network potential energy surface (PES) in the reactant asymptotic part. Only the reactant asymptotic part of the PES is constructed because of its barrierless and exothermic properties. Experimentally, the reaction rate coefficients are measured using a supersonic flow reactor. Pulsed laser photolysis of O3 molecules is used as the source of O(1D) atoms, which are detected directly through vacuum ultraviolet laser induced fluorescence at 115 nm. The branching ratio for H atom production is measured by comparing the H atom yields of the O(1D) + CH4 and O(1D) + H2 reactions. At T ≥ 75 K, good agreement between theoretical and experimental rate coefficients is found, while at 50 K, the larger difference is discussed in detail.