Path Integral Molecular Dynamics Research Articles

Isotope Quantum Effects in the Metallization Transition in Liquid Hydrogen.

Quantum effects in condensed matter normally only occur at low temperatures. Here we show a large quantum effect in high-pressure liquid hydrogen at thousands of Kelvins. We show that the metallization transition in hydrogen is subject to a very large isotope effect, occurring hundreds of degrees lower than the equivalent transition in deuterium. We examined this using path integral molecular dynamics simulations which identify a liquid-liquid transition involving atomization, metallization, and changes in viscosity, specific heat, and compressibility. The difference between H_{2} and D_{2} is a quantum mechanical effect that can be associated with the larger zero-point energy in H_{2} weakening the covalent bond. Our results mean that experimental results on deuterium must be corrected before they are relevant to understanding hydrogen at planetary conditions.

Multi-Level Monte Carlo Path Integral Molecular Dynamics for Thermal Average Calculation in the Nonadiabatic Regime

With the path integral approach, the thermal average in a multi-electronic-state quantum systems can be approximated by the ring polymer representation on an extended configuration space, where the additional degrees of freedom are associated with the surface index of each bead. The primary goal of this work is to propose a more efficient sampling algorithm for the calculation of such thermal averages. We reformulate the extended ring polymer approximation according to the configurations of the surface indexes, and by introducing a proper reference measure, the reformulation is recast as a ratio of two expectations of function expansions. By quantitatively estimating the sub-estimators, and minimizing the total variance of the sampled average, we propose a multi-level Monte Carlo path integral molecular dynamics method (MLMC-PIMD) to achieve an optimal balance of computational cost and accuracy.

Isotopic effects in chair graphane

Graphane is a layered material consisting of a sheet of hydrogenated graphene, with a C:H ratio of 1:1. We study isotopic effects in the properties of chair graphane, where H atoms alternate in a chairlike arrangement on both sides of the carbon layer. We use path-integral molecular dynamics simulations, which allows one to analyze the influence of nuclear quantum effects on equilibrium variables of materials. Finite-temperature properties of graphane are studied in the range 50–1500 K as functions of the isotopic mass of the constituent atoms, using an efficient tight-bonding potential. Results are presented for kinetic and internal energy, atomic mean-square displacements, fluctuations in the C–H bond direction, plus interatomic distances and layer area. At low temperature, substituting 13C for 12C gives a fractional change of −2.6 × 10−4 in C–C distance and −3.9 × 10−4 in the graphane layer area. Replacing 2H for 1H causes a larger fractional change in the C–H bond of −5.7 × 10−3. The isotopic effect in C–C bond distance increases (decreases) by applying a tensile (compressive) in-plane stress. These results are interpreted in terms of a quasiharmonic approximation for the vibrational modes. Similarities and differences with isotopic effects in graphene are discussed.

Non-adiabatic ring polymer molecular dynamics with spin mapping variables.

We present a new non-adiabatic ring polymer molecular dynamics (NRPMD) method based on the spin mapping formalism, which we refer to as the spin mapping NRPMD (SM-NRPMD) approach. We derive the path-integral partition function expression using the spin coherent state basis for the electronic states and the ring polymer formalism for the nuclear degrees of freedom. This partition function provides an efficient sampling of the quantum statistics. Using the basic properties of the Stratonovich-Weyl transformation, we further justify a Hamiltonian that we propose for the dynamical propagation of the coupled spin mapping variables and the nuclear ring polymer. The accuracy of the SM-NRPMD method is numerically demonstrated by computing the nuclear position and population auto-correlation functions of non-adiabatic model systems. The results obtained using the SM-NRPMD method agree very well with the numerically exact results. The main advantage of using the spin mapping variables over the harmonic oscillator mapping variables is numerically demonstrated, where the former provides nearly time-independent expectation values of physical observables for systems under thermal equilibrium. We also explicitly demonstrate that SM-NRPMD provides invariant dynamics upon various ways of partitioning the state-dependent and state-independent potentials.

Nuclear quantum effects in thermal conductivity from centroid molecular dynamics.

We show that the centroid molecular dynamics (CMD) method provides a realistic way to calculate the thermal diffusivity a = λ/ρcV of a quantum mechanical liquid such as para-hydrogen. Once a has been calculated, the thermal conductivity can be obtained from λ = ρcVa, where ρ is the density of the liquid and cV is the constant-volume heat capacity. The use of this formula requires an accurate quantum mechanical heat capacity cV, which can be obtained from a path integral molecular dynamics simulation. The thermal diffusivity can be calculated either from the decay of the equilibrium density fluctuations in the liquid or by using the Green-Kubo relation to calculate the CMD approximation to λ and then dividing this by the corresponding approximation to ρcV. We show that both approaches give the same results for liquid para-hydrogen and that these results are in good agreement with the experimental measurements of the thermal conductivity over a wide temperature range. In particular, they correctly predict a decrease in the thermal conductivity at low temperatures-an effect that stems from the decrease in the quantum mechanical heat capacity and has eluded previous para-hydrogen simulations. We also show that the method gives equally good agreement with the experimental measurements for the thermal conductivity of normal liquid helium.

Relationship between Hydrogen-Bonding Motifs and the 1b1 Splitting in the X-ray Emission Spectrum of Liquid Water.

The split of the 1b1 peak observed in the X-ray emission (XE) spectrum of liquid water has been the focus of intense research. Although several hypotheses have been proposed to explain the origin of this split, a consensus has not yet been reached. Here, we introduce a novel theoretical/computation approach which, combining path-integral molecular dynamics simulations with the MB-pol model and time-dependent density functional theory calculations, predicts the 1b1 splitting in liquid water and not in crystalline ice, in agreement with the experimental observations. A systematic analysis of the underlying local structure of liquid water at ambient conditions indicates that several different hydrogen-bonding motifs contribute to the overall XE line shape in the energy range corresponding to emissions from the 1b1 orbitals. This suggests that it is not possible to unambiguously attribute the split of the 1b1 peak to only two specific structural arrangements of the underlying hydrogen-bonding network.

Gibbs Energy of Hydrogen Adsorption on Pt Surface by Machine Learning Potential and Metadynamics

A machine learning potential that can describe the ethylene hydrogenation reaction on Pt surface has been developed. Using the obtained potential, path integral molecular dynamics of hydrogen on Pt...

Non-adiabatic Matsubara dynamics and non-adiabatic ring-polymer molecular dynamics.

We present the non-adiabatic Matsubara dynamics, a general framework for computing the time-correlation function (TCF) of electronically non-adiabatic systems. This new formalism is derived based on the generalized Kubo-transformed TCF using the Wigner representation for both the nuclear degrees of freedom and the electronic mapping variables. By dropping the non-Matsubara nuclear normal modes in the quantum Liouvillian and explicitly integrating these modes out from the expression of the TCF, we derived the non-adiabatic Matsubara dynamics approach. Further making the approximation to drop the imaginary part of the Matsubara Liouvillian and enforce the nuclear momentum integral to be real, we arrived at the non-adiabatic ring-polymer molecular dynamics (NRPMD) approach. We have further justified the capability of NRPMD for simulating the non-equilibrium TCF. This work provides the rigorous theoretical foundation for several recently proposed state-dependent RPMD approaches and offers a general framework for developing new non-adiabatic quantum dynamics methods in the future.

Assessing parameters for ring polymer molecular dynamics simulations at low temperatures: DH + H chemical reaction

Ring polymer molecular dynamics (RPMD) is an accurate method for calculating thermal chemical reaction rates. It has recently been discovered that low-temperature calculations are strongly affected by the simulation parameters. Here, for the thermally activated reaction DH + H → D + H2, we calculate the RPMD rate constants at T=50, 100, and 300 K and demonstrate that for T≥100 K the standard input parameters yield accurate results, but at low temperatures (e.g., 50 K) one must increase the asymptotic distance and force constant, and decrease the umbrella integration step.

Advances in methods and applications of nonadiabatic quantum dynamics simulation of condensed matters

First-principles molecular dynamics simulations, which can reveal the microscopic mechanisms and the macroscopic properties of molecules and condensed systems, have become a significant driving force in the development of physics. In the past decades, the ground state electronic properties of various materials can be accurately described using the first-principles method based on the Born-Oppenheimer approximation. However, in nonadiabatic processes, such as the ultrafast excited-state dynamics where laser interacts with molecules and condensed matters, the Born-Oppenheimer approximation evolving nuclear wave function on the ground state potential energy plane is no longer valid because of the nonnegligible nonadiabatic couplings between electrons and nuclei. In order to investigate the novel physical phenomena in nonadiabatic processes, many nonadiabatic molecular dynamics methods have been proposed, for instance, full quantum dynamics and mixed quantum-classical dynamics. When the nuclear quantum effects are trivial, the mixed quantum-classical dynamics is effective, which solves the time-dependent electron Schrodinger equation and describes the motion of the nuclei in the form of Newtonian mechanics. Two methods based on the mixed quantum-classical dynamics, the fewest-switches trajectory surface-hopping method (FSSH) and Ehrenfest dynamics, have accepted much attention since they could be easily implemented and combined with real-time time-dependent density functional theory (rt-TDDFT) for high-precision calculations. The combination of quantum and classical dynamics reduces the computational cost and allows the applications for real materials. At low temperatures, the nuclear quantum effects, including quantum tunneling and zero-point vibrations, cannot be ignored. Ring-polymer molecular dynamics (RPMD) based on the imaginary-time path integral is widely used to consider the nuclear quantum effects, which describes the quantum behavior of atomic nuclei by sampling a larger number of quantum configurations and paths. By combining RPMD and mixed quantum-classical nonadiabatic dynamics, one can describe quantum motions of nuclei and the excited electrons in quantum nonadiabatic dynamics simulations. Here, we summarize the recent progress of the rt-TDDFT methods based on numerical atomic orbital basis and plane- wave basis set, as well as highly efficient RP-TDAP method combining rt-TDDFT and RPMD towards a quantum description of electronic-nuclear dynamics. Finally, we show several representative applications employing these methods, including high harmonic generation modulated by two-color light, charge density wave dynamics, photocatalytic water decomposition as well as the quantum evolution of nuclear wave packets. These applications are in good agreement with experimental measurements, which demonstrates the reliability of these methods and their advantages in the corresponding research areas. These developments and applications are in a significant step from ground state methods to full quantum simulation of the coupled motion of electrons and nuclei, providing a new perspective for understanding and predicting the quantum interactions and dynamical behavior of condensed materials in the degree of atomic spatial scales and attosecond time scales.

Persistent Homology Metrics Reveal Quantum Fluctuations and Reactive Atoms in Path Integral Dynamics.

Nuclear quantum effects (NQEs) are known to impact a number of features associated with chemical reactivity and physicochemical properties, particularly for light atoms and at low temperatures. In the imaginary time path integral formalism, each atom is mapped onto a “ring polymer” whose spread is related to the quantum mechanical uncertainty in the particle’s position, i.e., its thermal wavelength. A number of metrics have previously been used to investigate and characterize this spread and explain effects arising from quantum delocalization, zero-point energy, and tunneling. Many of these shape metrics consider just the instantaneous structure of the ring polymers. However, given the significant interest in methods such as centroid molecular dynamics and ring polymer molecular dynamics that link the molecular dynamics of these ring polymers to real time properties, there exists significant opportunity to exploit metrics that also allow for the study of the fluctuations of the atom delocalization in time. Here we consider the ring polymer delocalization from the perspective of computational topology, specifically persistent homology, which describes the 3-dimensional arrangement of point cloud data, (i.e. atomic positions). We employ the Betti sequence probability distribution to define the ensemble of shapes adopted by the ring polymer. The Wasserstein distances of Betti sequences adjacent in time are used to characterize fluctuations in shape, where the Fourier transform and associated principal components provides added information differentiating atoms with different NQEs based on their dynamic properties. We demonstrate this methodology on two representative systems, a glassy system consisting of two atom types with dramatically different de Broglie thermal wavelengths, and ab initio molecular dynamics simulation of an aqueous 4 M HCl solution where the H-atoms are differentiated based on their participation in proton transfer reactions.

Neural network potential energy surface for the low temperature ring polymer molecular dynamics of the H2CO + OH reaction.

A new potential energy surface (PES) and dynamical study of the reactive process of H2CO + OH toward the formation of HCO + H2O and HCOOH + H are presented. In this work, a source of spurious long range interactions in symmetry adapted neural network (NN) schemes is identified, which prevents their direct application for low temperature dynamical studies. For this reason, a partition of the PES into a diabatic matrix plus a NN many-body term has been used, fitted with a novel artificial neural network scheme that prevents spurious asymptotic interactions. Quasi-classical trajectory (QCT) and ring polymer molecular dynamics (RPMD) studies have been carried on this PES to evaluate the rate constant temperature dependence for the different reactive processes, showing good agreement with the available experimental data. Of special interest is the analysis of the previously identified trapping mechanism in the RPMD study, which can be attributed to spurious resonances associated with excitations of the normal modes of the ring polymer.

Dynamics of the isotope exchange reaction of D with H3 +, H2D+, and D2H.

We have measured the merged-beams rate coefficient for the titular isotope exchange reactions as a function of the relative collision energy in the range of ∼3 meV-10 eV. The results appear to scale with the number of available sites for deuteration. We have performed extensive theoretical calculations to characterize the zero-point energy corrected reaction path. Vibrationally adiabatic minimum energy paths were obtained using a combination of unrestricted quadratic configuration interaction of single and double excitations and internally contracted multireference configuration interaction calculations. The resulting barrier height, ranging from 68 meV to 89 meV, together with the various asymptotes that may be reached in the collision, was used in a classical over-the-barrier model. All competing endoergic reaction channels were taken into account using a flux reduction factor. This model reproduces all three experimental sets quite satisfactorily. In order to generate thermal rate coefficients down to 10 K, the internal excitation energy distribution of each H3 + isotopologue is evaluated level by level using available line lists and accurate spectroscopic parameters. Tunneling is accounted for by a direct inclusion of the exact quantum tunneling probability in the evaluation of the cross section. We derive a thermal rate coefficient of <1×10-12 cm3 s-1 for temperatures below 44 K, 86 K, and 139 K for the reaction of D with H3 +, H2D+, and D2H+, respectively, with tunneling effects included. The derived thermal rate coefficients exceed the ring polymer molecular dynamics prediction of Bulut et al. [J. Phys. Chem. A 123, 8766 (2019)] at all temperatures.

Nuclear quantum effects on autoionization of water isotopologs studied by ab initio path integral molecular dynamics.

In this study, we investigate the nuclear quantum effects (NQEs) on the acidity constant (pKA) of liquid water isotopologs under the ambient condition by path integral molecular dynamics (PIMD) simulations. We compared simulations using a fully explicit solvent model with a classical polarizable force field, density functional tight binding, and ab initio density functional theory, which correspond to empirical, semiempirical, and ab initio PIMD simulations, respectively. The centroid variable with respect to the proton coordination number of a water molecule was restrained to compute the gradient of the free energy, which measures the reversible work of the proton abstraction for the quantum mechanical system. The free energy curve obtained by thermodynamic integration was used to compute the pKA value based on probabilistic determination. This technique not only reproduces the pKA value of liquid D2O experimentally measured (14.86) but also allows for a theoretical prediction of the pKA values of liquid T2O and aqueous HDO and HTO, which are unknown due to their scarcity. It is also shown that the NQEs on the free energy curve can result in a downshift of 4.5 ± 0.9 pKA units in the case of liquid water, which indicates that the NQEs plays an indispensable role in the absolute determination of pKA. The results of this study can help inform further extensions into the calculation of the acidity constants of isotope substituted species with high accuracy.

Small Nuclear Quantum Effects in Scattering of H and D from Graphene.

We study nuclear quantum effects in H/D sticking to graphene, comparing scattering experiments at near-zero coverage with classical, quantized, and transition-state calculations. The experiment shows H/D sticking probabilities that are indistinguishable from one another and markedly smaller than those expected from a consideration of zero-point energy shifts of the chemisorption transition state. Inclusion of dynamical effects and vibrational anharmonicity via ring-polymer molecular dynamics (RPMD) yields results that are in good agreement with the experimental results. RPMD also reveals that nuclear quantum effects, while modest, arise primarily from carbon and not from H/D motion, confirming the importance of a C atom rehybridization mechanism associated with H/D sticking on graphene.

The metastable structures of anthracene-argon clusters inside helium nanodroplets

Recent experiments have revealed that anthracene forms rather diverse complexes with multiple argon atoms upon embedding in helium nanodroplets. In this work we have modeled such complexes and investigated the relative stability of various isomeric forms using classical force fields reparametrized on electronic structure data at the level of spin-component scaled MP2, nuclear quantum effects being accounted for by means of path-integral molecular dynamics simulations. We generally find that the structures obtained are very sensitive to the details of the potential, but also to the inclusion of zero-point motion, for which the harmonic approximation is found to be generally accurate. Inside He $$_{1000}$$ clusters, the complexes also exhibit different relative stabilities, the quantum solvent usually favoring more asymmetric complexes. No evidence for fluxionality at 1 K in a given cluster is found. Using ring-polymer molecular dynamics trajectories, the stability of an argon dimer initially distant from the anthracene molecule was also investigated. Such metastable configurations are found to have a lifetime of the order of nanoseconds, confirming the possibility of trapping them in experiments.

Quantum Nuclei at Weakly Bonded Interfaces: The Case of Cyclohexane on Rh(111)

AbstractThe electronic properties of interfaces can depend on their isotopic constitution. One known case is that of cyclohexane physisorbed on Rh(111), in which isotope effects have been measured on the work function change and desorption energies. These effects can only be captured by calculations including nuclear quantum effects (NQE). In this paper, this interface is addressed employing dispersion‐inclusive density‐functional theory coupled to a quasi‐harmonic (QH) approximation for NQE, as well as to fully anharmonic ab initio path integral molecular dynamics (PIMD). The QH approximation is able to capture that deuterated cyclohexane has a smaller adsorption energy and lies about 0.01 Å farther from the Rh(111) surface than its isotopologue, which can be correlated to the isotope effect in the work function change. An investigation of the validity of the QH approximation relying on PIMD simulations, leads to the conclusion that although this interface is highly impacted by anharmonic quantum fluctuations in the molecular layer and at bonding sites, these anharmonic contributions play a minor role when analyzing isotope effects at low temperatures. Nevertheless, anharmonic quantum fluctuations cause an increase in the distance between the molecular layer and Rh(111), a consequent smaller overall work function change, and intricate changes in orbital hybridization.

A generalized class of strongly stable and dimension-free T-RPMD integrators.

Recent work shows that strong stability and dimensionality freedom are essential for robust numerical integration of thermostatted ring-polymer molecular dynamics (T-RPMD) and path-integral molecular dynamics, without which standard integrators exhibit non-ergodicity and other pathologies [R. Korol et al., J. Chem. Phys. 151, 124103 (2019) and R. Korol et al., J. Chem. Phys. 152, 104102 (2020)]. In particular, the BCOCB scheme, obtained via Cayley modification of the standard BAOAB scheme, features a simple reparametrization of the free ring-polymer sub-step that confers strong stability and dimensionality freedom and has been shown to yield excellent numerical accuracy in condensed-phase systems with large time steps. Here, we introduce a broader class of T-RPMD numerical integrators that exhibit strong stability and dimensionality freedom, irrespective of the Ornstein-Uhlenbeck friction schedule. In addition to considering equilibrium accuracy and time step stability as in previous work, we evaluate the integrators on the basis of their rates of convergence to equilibrium and their efficiency at evaluating equilibrium expectation values. Within the generalized class, we find BCOCB to be superior with respect to accuracy and efficiency for various configuration-dependent observables, although other integrators within the generalized class perform better for velocity-dependent quantities. Extensive numerical evidence indicates that the stated performance guarantees hold for the strongly anharmonic case of liquid water. Both analytical and numerical results indicate that BCOCB excels over other known integrators in terms of accuracy, efficiency, and stability with respect to time step for practical applications.

Importance of nuclear quantum effects on the hydration of chloride ion

Path-integral ab initio molecular dynamics (PI-AIMD) calculations have been employed to probe the nature of chloride ion solvation in aqueous solution. Nuclear quantum effects (NQEs) are shown to weaken hydrogen bonding between the chloride anion and the solvation shell of water molecules. As a consequence, the disruptive effect of the anion on the solvent water structure is significantly reduced compared to what is found in the absence of NQEs. The chloride hydration structure obtained from PI-AIMD agrees well with information extracted from neutron scattering data. Inparticular, the observed satellite peak in the hydrogen-chloride-hydrogen triple angular distribution serves as a clear signature of NQEs. The present results suggest that NQEs are likely to play acrucial role in determining the structure of saline solutions.

First-principles equation of state database for warm dense matter computation.

We put together a first-principles equation of state (FPEOS) database for matter at extreme conditions by combining results from path integral Monte Carlo and density functional molecular dynamics simulations of the elements H, He, B, C, N, O, Ne, Na, Mg, Al, and Si as well as the compounds LiF,B_{4}C,BN,CH_{4},CH_{2},C_{2}H_{3},CH,C_{2}H,MgO,andMgSiO_{3}. For all these materials, we provide the pressure and internal energy over a density-temperature range from ∼0.5 to 50 gcm^{-3} and from ∼10^{4} to 10^{9}K, which are based on ∼5000 different first-principles simulations. We compute isobars, adiabats, and shock Hugoniot curves in the regime of L- and K-shell ionization. Invoking the linear mixing approximation, we study the properties of mixtures at high density and temperature. We derive the Hugoniot curves for water and alumina as well as for carbon-oxygen, helium-neon, and CH-silicon mixtures. We predict the maximal shock compression ratios of H_{2}O, H_{2}O_{2}, Al_{2}O_{3},CO,andCO_{2} to be 4.61, 4.64, 4.64, 4.89, and 4.83, respectively. Finally we use the FPEOS database to determine the points of maximum shock compression for all available binary mixtures. We identify mixtures that reach higher shock compression ratios than their end members. We discuss trends common to all mixtures in pressure-temperature and particle-shock velocity spaces. In the Supplemental Material, we provide all FPEOS tables as well as computer codes for interpolation, Hugoniot calculations, and plots of various thermodynamic functions.