Path Integral Molecular Dynamics Research Articles

Selecting Initial Conditions for Trajectory-Based Nonadiabatic Simulations.

ConspectusPhotochemical reactions have always been the source of a great deal of mystery. While classified as a type of chemical reaction, no doubts are allowed that the general tenets of ground-state chemistry do not directly apply to photochemical reactions. For a typical chemical reaction, understanding the critical points of the ground-state potential (free) energy surface and embedding them in a thermodynamics framework is often enough to infer reaction yields or characteristic time scales. A general working principle is that the energy profile along the minimum energy paths provides the key information to characterize the reaction. These well-developed concepts, unfortunately, rarely stretch to processes involving the formation of a nonstationary state for a molecular system after light absorption.Upon photoexcitation, a molecule is likely to undergo internal conversion processes, that is, changes of electronic states mediated by couplings between nuclear and electronic motion, precisely what the celebrated Born-Oppenheimer approximation neglects. These coupled electron-nuclear processes, coined nonadiabatic processes, allow for the molecule to decay from one electronic state to the other nonradiatively. Understanding the intricate nonadiabatic dynamics is pivotal to rationalizing and predicting the outcome of a molecular photoexcitation and providing insights for experiments conducted, for example, in advanced light sources such as free-electron lasers.Nowadays, most simulations in nonadiabatic molecular dynamics are based on approximations that invoke a near-classical depiction of the nuclei. This reliance is due to practical constraints, and the classical equations of motion for the nuclei must be supplemented by techniques such as surface hopping to account for nonadiabatic transitions between electronic states. A critical but often overlooked aspect of these simulations is the selection of initial conditions, specifically the choice of initial nuclear positions and momenta for the nonadiabatic dynamics, which can significantly influence how well the simulations mimic real quantum systems across various experimental scenarios. The conventional approach for generating initial conditions for nonadiabatic dynamics typically maps the initial state onto a nuclear phase space using a Wigner quasiprobability function within a harmonic approximation, followed by a second approximation where the molecule undergoes a sudden excitation.In this Account, we aim to warn the experienced or potential user of nonadiabatic molecular dynamics about the possible limitations of this strategy for initial-condition generation and its inability to accurately describe the photoexcitation of a molecule. More specifically, we argue that the initial phase-space distribution can be more accurately represented through molecular dynamics simulations by using a quantum thermostat. This method offers a robust framework that can be applied to large, flexible, or even solvated molecular systems. Furthermore, the reliability of this strategy can be benchmarked against more rigorous approaches such as path integral molecular dynamics. Additionally, the commonly used sudden approximation, which assumes a vertical and sudden excitation of a molecule, rarely reflects the excitation triggered by laser pulses used in actual photochemical and spectroscopic experiments. We discuss here a more general approach that can generate initial conditions for any type of laser pulse. We also discuss strategies to tackle excitation triggered by a continuous-wave laser.

Full-dimensional accurate potential energy surface and dynamics for the unimolecular isomerization reaction CH3NC ⇌ CH3CN.

The reaction CH3NC ⇌ CH3CN, a model reaction for the study of unimolecular isomerization, is important in astronomy and atmospheric chemistry and has long been studied by numerous experiments and theories. In this work, we report the first full-dimensional accurate potential energy surface (PES) of this reaction by the permutation invariant polynomial-neural network method based on 30 974 points, whose energies are calculated at the CCSD(T)-F12a/AVTZ level. Then, ring polymer molecular dynamics is used to derive the free energy barrier of the reaction at the experimental temperature range of 472.55-532.92K. Reaction kinetics are studied at the high-pressure limit and in the fall-off region by standard transition state theory and the master equation, respectively. The calculated temperature- and pressure-dependent rate coefficients are in good agreement with previous experimental and theoretical results. Furthermore, quasi-classical trajectory simulations are performed on this PES to study the intramolecular energy transfer dynamics at initial vibrational energies of 4.336, 5.204, and 6.505eV.

Quadratic scaling path integral molecular dynamics for fictitious identical particles and its application to fermion systems

Quadratic scaling path integral molecular dynamics for fictitious identical particles and its application to fermion systems

Nuclear quantum effects on glassy water under pressure: Vitrification and pressure-induced transformations.

We perform classical molecular dynamics (MD) and path-integral MD (PIMD) simulations of H2O and D2O using the q-TIP4P/F model over a wide range of temperatures and pressures to study the nuclear quantum effects (NQEs) on (i) the vitrification of liquid water upon isobaric cooling at different pressures and (ii) pressure-induced transformations at constant temperature between low-density amorphous and high-density amorphous ice (LDA and HDA) and hexagonal ice Ih and HDA. Upon isobaric cooling, classical and quantum H2O and D2O vitrify into a continuum of intermediate amorphous ices (IA), with densities in-between those of LDA and HDA (depending on pressure). Importantly, the density of the IA varies considerably if NQEs are included (similar conclusions hold for ice Ih at all pressures studied). While the structure of the IA is not very sensitive to NQE, the geometry of the hydrogen-bond (HB) is. NQE leads to longer and less linear HB in LDA, HDA, and ice Ih than found in the classical case. Interestingly, the delocalization of the H/D atoms is non-negligible and identical in LDA, HDA, and ice Ih at all pressures studied. Our isothermal compression/decompression MD/PIMD simulations show that classical and quantum H2O and D2O all exhibit LDA-HDA and ice Ih-HDA transformations, consistent with experiments. The inclusion of NQE leads to a softer HB-network, which lowers slightly the LDA/ice Ih-to-HDA transformation pressures. Interestingly, the HB in HDA is longer and less linear than in LDA, which is counterintuitive given that HDA is ≈25% denser than LDA. Overall, our results show that, while classical computer simulations provide the correct qualitative phenomenology of ice and glassy water, NQEs are necessary for a quantitative description.

A novel mechanism and model of hydrogen isotope exchange in vacancy considering nuclear quantum effects

Abstract Tritium (T) retention in plasma-facing materials (PFMs) raises significant radiological safety concerns and adversely affects the self-sustained burning of T in fusion reactors. Therefore, the removal of retained T from PFMs has become an urgent task. Hydrogen isotope (HI) exchange has proven to be an effective method for T removal. However, the microscopic mechanisms, particularly the role of isotope effects, remain insufficiently understood. For the first time, this work employs path integral molecular dynamics to account for the nuclear quantum effects of HIs and explore the microscopic mechanisms of HI exchange in tungsten (W) vacancies. Atomic-scale simulations reveal that the fundamental principle of HI exchange is the reduction in binding energy between HIs and vacancies as the defect filling level increases, involving two key processes: de-trapping and trapping. These processes are influenced by isotope effects, with pronounced differences observed at low temperatures. Notably, T exhibits a higher probability of de-trapping from vacancies compared to hydrogen (H), while vacancies demonstrate a stronger affinity for trapping H. In contrast, the isotope effect is less pronounced between deuterium (D) and H, leading to different exchange efficiencies between H and T, as well as between D and T. Based on these isotope effects, we proposed a detailed microscopic mechanism for HI exchange in vacancies and developed a T evolution model that accounts for variations in HI concentrations. This work advances our understanding of HI exchange for T removal applications and offers a more accurate assessment of T retention in future D–T environments.

Potential energy landscape formalism for quantum molecular liquids

The potential energy landscape (PEL) formalism is a powerful tool within statistical mechanics to study the thermodynamic properties of classical low-temperature liquids and glasses. Recently, the PEL formalism has been extended to liquids/glasses that obey quantum mechanics, but applications have been limited to atomistic model liquids. In this work, we extend the PEL formalism to liquid/glassy water using path-integral molecular dynamics (PIMD) simulations, where nuclear quantum effects (NQE) are included. Our PIMD simulations, based on the q-TIP4P/F water model, show that the PEL of quantum water is both Gaussian and anharmonic. Importantly, the ring-polymers associated to the O/H atoms in the PIMD simulations, collapse at the local minima of the PEL (inherent structures, IS) for both liquid and glassy states. This allows us to calculate, analytically, the IS vibrational density of states (IS-VDOS) of the ring-polymer system using the IS-VDOS of classical water (obtained from classical MD simulations). The role of NQE on the structural properties of liquid/glassy water at various pressures are discussed in detail. Overall, our results demonstrate that the PEL formalism can effectively describe the behavior of molecular liquids at low temperatures and in the glass states, regardless of whether the liquid/glass obeys classical or quantum mechanics.

Predictive evaluation of hydrogen diffusion coefficient on Pd(111) surface by path integral simulations using neural network potential

Quantum diffusion of hydrogen (H) on palladium (Pd) (111) was studied using two types of path integral simulations: quantum transition state theory (QTST) and ring polymer molecular dynamics (RPMD). The use of an artificial neural network potential trained by density functional theory calculations has made it feasible to perform path integral simulations considering nuclear quantum effects (NQEs) of a many-body Pd-H system with accuracy. The QTST result has shown a clear non-Arrhenius behavior in the diffusion coefficient (D) of H below the temperature of 150–200 K due to the NQEs. Comparing the D on Pd surface and bulk Pd, it was found that surface and bulk diffusion are competitive at high temperature. Consistent with this, it was observed in the RPMD simulation at 800 K that a part of the quantum trajectories of H on the surface bifurcates to the Pd subsurface. As the temperature decreases, the surface diffusion becomes much faster than the bulk diffusion. Furthermore, distinct quantum behaviors of H were identified at the surface and within the bulk. The D values obtained from this study using a “flexible” Pd surface model differed considerably from those previously reported using “frozen” Pd surface models, indicating an important contribution from Pd vibrations coupled with H diffusion. Published by the American Physical Society 2024

Nuclear Quantum Effects Enhance Structural Stability but Accelerate Charge Carrier Recombination in MHyPbBr3 Perovskite.

Hybrid organic-inorganic perovskites exhibit significant nuclear quantum effects (NQEs) due to their light hydrogen atoms. By performing ring polymer molecular dynamics, ab initio molecular dynamics, and nonadiabatic molecular dynamics simulations on the MHyPbBr3 (MHy+ = CH3NH2NH2+) perovskites, we demonstrate that NQEs stabilize the lattice by suppressing atomic motions and accelerate nonradiative charge recombination. This stabilization arises from the synergistic effects of the Pb-N coordination bonds and N-H···Br hydrogen bonds, which enhance organic-inorganic interactions. As a result, Pb-Br octahedra, particularly [Pb(1)Br6]4- octahedra supporting electron and hole, are well-preserved, promoting electronic wavefunction delocalization and increasing electron-hole overlap. These effects enhance nonadiabatic coupling by overcoming the reduced atomic motions. Overall, this and the prolonged decoherence time accelerate the nonradiative electron-hole recombination due to NQEs. Our study highlights the unique influence of NQEs on the geometrical stability and charge carrier dynamics in MHyPbBr3, offering fundamental insights for future material design.

Hydrogen diffusion on Ni(100): A combined machine-learning, ring polymer molecular dynamics, and kinetic Monte Carlo study.

We introduce a methodological framework coupling machine-learning potentials, ring polymer molecular dynamics (RPMD), and kinetic Monte Carlo (kMC) to draw a comprehensive physical picture of the collective diffusion of hydrogen atoms on metal surfaces. For the benchmark case of hydrogen diffusion on a Ni(100) surface, the hydrogen adsorption and diffusion energetics and its dependence on the local coverage is described via a neural-network potential, where the training data are computed via periodic density functional theory (DFT) and include all relevant optimized diffusion and desorption paths, sampled by nudged elastic band optimizations and molecular dynamics simulations. Nuclear quantum effects, being crucial for processes involving hydrogen at low temperatures, are treated by RPMD. The diffusion rate constants are calculated with a combination of umbrella samplings employed to map the free energy profile and separate samplings of recrossing trajectories to obtain the transmission coefficient. The calculated diffusion rates for different temperatures and local environments are then combined and fitted into a kMC model allowing access to larger time and length scales. Our results demonstrate an outstanding performance for the trained neural network potential in reproducing reference DFT energies and forces. We report the effective diffusion rates for different temperatures and hydrogen surface coverages obtained via this recipe in good agreement with the experimental results. The method combination proposed in this study can be instrumental for a wide range of applications in materials science.

A theoretical study on muoniated N-heterocyclic carbenes using path integral molecular dynamics.

Several N-heterocyclic carbenes (NHCs) are experimentally observed upon the addition of muonium (Mu), and the hyperfine coupling constants (HFCCs) of muon are measured. Theoretical investigation of Mu has been challenging due to significant quantum effects. Herein, we performed an abinitio path integral molecular dynamics (PIMD) simulation, which accurately considers multi-dimensional quantum effects, to theoretically investigate muoniated 1,3-dihydro-2H-imidazole-2-ylidene (Mu-IY). Our findings indicate that quantum effects have two contradictory contributions: the quantum effect of bond vibrations increases the HFCC values, whereas that of out-of-plane angular vibrations decreases the HFCC values. Moreover, we show that the HFCC values of other NHCs can be predicted without the PIMD simulations by applying the structural changes caused by the quantum effect derived from the PIMD simulations of Mu-IY.

Generalized Path Integral Energy and Heat Capacity Estimators of Quantum Oscillators and Crystals Using Harmonic Mapping.

Imaginary-time path integral (PI) is a rigorous tool to treat nuclear quantum effects in static properties. However, with its high computational demand, it is crucial to devise precise estimators. We introduce generalized PI estimators for the energy and heat capacity that utilize coordinate mapping. While it can reduce to the standard thermodynamic and centroid virial (CVir) estimators, the formulation can also take advantage of harmonic character of quantum oscillators and crystals to construct a coordinate mapping. The method is not applicable to fluids or systems with imaginary modes such as double-well potentials. This yields harmonically mapped averaging (HMA) estimators, with mappings that decouple (HMAc) or couple (HMAq) the centroid and internal modes. The HMAq is constructed with normal mode coordinates (HMAq-NM) with quadratic scaling of cost or harmonic oscillator staging (HMAq-SG) coordinates with linear scaling. The estimator performance is examined for a 1D anharmonic oscillator and a 3D Lennard-Jones crystal using path integral molecular dynamics (PIMD) simulation. The HMA estimators consistently provide more precise estimates compared to CVir, with the best performance obtained by HMAq-NM, followed by HMAq-SG, and then HMAc. We also examine the effect of anharmonicity (for AO), intrinsic quantumness, and Trotter number. The HMA formulation introduced assumes the availability of forces and Hessian matrix; however, an equally efficient finite difference alternative is possible when these derivatives are inaccessible. The remarkable improvement in precision offered by HMAq estimators provides a framework for efficient PI simulation of more challenging systems, such as those based on ab initio calculations.

Collision of cesium atoms on helium nanodroplets: Unraveling mechanisms for surface capture at experimental velocities.

The collision of cesium atoms on the surface of helium nanodroplets (HNDs) containing 1000 atoms is described by the ZPAD-mPL approach, a zero-point averaged dynamics (ZPAD) method based on a He-He pseudopotential adjusted to better reproduce the total energy of He1000. Four types of collisional patterns were identified depending on the initial projectile speed v0 and impact parameter b. At the lowest speeds (v0 ≲ 250ms-1), Cs atoms are softly captured by the HND surface, while at the highest ones (v0 ≳ 500-600ms-1), Cs atoms can travel through the droplet and move away. In between these two extreme cases, Cs atoms can be temporarily submerged in the HND before being expelled to the surface if b = 0 or cross the HND before being captured on its surface. The possibility for Cs capture at experimental velocities and droplet piercings at the highest ones contrasts with time-dependent density functional theory calculations, which predict Cs capture for velocities lower than 75ms-1, and ring-polymer molecular dynamics (RPMD) or former ZPAD-like methods, which predict soft Cs capture up to 500ms-1. ZPAD-mPL results are attributed to the liquid but non-superfluid nature of the droplet, which favors energy exchanges with the helium environment, and to low He-He binding energy and HND surface tension, which can stimulate helium ejections, especially at high projectile speeds. Despite the use of a pseudopotential to model He-He interactions, the heliophobicity of Cs atoms is maintained as demonstrated by their ability to remain localized on the HND surface or to be expelled to the HND surface after transient submersion in helium.

Thermal Rates and High-Temperature Tunneling from Surface Reaction Dynamics and First-Principles.

Studying dynamics of the dissociative adsorption and recombinative desorption of hydrogen on copper surfaces has shaped our atomic-scale understanding of surface chemistry, yet experimentally determining the thermal rates for these processes, which dictate the outcome of catalytic reactions, has been impossible so far. In this work, we determine the thermal rate constants for dissociative adsorption and recombinative desorption of hydrogen on Cu(111) between 200 and 1000 K using data from reaction dynamics experiments. Contrary to current understanding, our findings demonstrate the predominant role of quantum tunneling, even at temperatures as high as 400 K. We also provide precise values for the reaction barrier (0.619 ± 0.020 eV) and adsorption energy (0.348 ± 0.026 eV) for H2 on Cu(111). Remarkably, the thermal rate constants are in excellent agreement with a first-principles quantum rate theory based on a new implementation of ring polymer molecular dynamics for reactions on surfaces, paving the way to discovering better catalysts using reliable and efficient computational methods.

Proton quantal delocalization and H/D translocations in (MeOH)nH+ (n = 2, 3).

In this study, we present results from path integral molecular dynamics simulations that describe the characteristics of the quantum spatial delocalizations of protons participating in OH bonds in (MeOH)2H+ and in (MeOH)3H+. The characterization was carried out by examining the overall structures of the corresponding isomorphic polymers. To introduce full flexibility in the force treatment, we have adopted a neural network fitting procedure based on second-order Møller-Plesset perturbation theory predictions. For the dimer case, we found that the spatial extent of the shared connective proton can be portrayed in terms of a prolate-like structure with typical dimensions of ∼0.1 Å. On the other hand, the dangling polymers lie confined within a thin spherical layer, spread over length scales of the order of ∼0.25 Å. In contrast, connective protons in (MeOH)3H+ exhibit larger delocalizations along the O-H bond and more localized ones along perpendicular directions, compared to their dangling counterparts. We also examined the characteristics of the relative propensities of H and D isotopes to be localized in dangling and connective positions. Physical interpretations of the different thermodynamic trends are provided in terms of the local geometrical characteristics and of the strengths of the corresponding intermolecular connectivities.

Accurate and Fast Estimation of the Continuum Limit in Path Integral Simulations of Quantum Oscillators and Crystals.

Convergence of imaginary-time path integral results requires a substantial number of beads (n) when quantum effects are significant. Traditional Trotter scaling approaches estimate the continuum limit (n → ∞) through extrapolation; however, their validity is restricted to the asymptotic domain of large n. We introduce an efficient extrapolation approach for quantum oscillators and crystals with harmonic character. The fitting function is inspired by the analytic solution of the harmonic oscillator (HO), or the Einstein crystal for solids. While the formulation is designed for first derivative properties (energy and pressure), its extension to second derivative properties (e.g., heat capacity) is straightforward. We apply the method to a one-dimensional HO and anharmonic oscillator (AO), as well as a three-dimensional Lennard-Jones crystal. Configurations are sampled using path integral molecular dynamics simulations in the canonical ensemble, at a low temperature of T = 0.1 (simulation units). Compared to Trotter extrapolation approaches, the new method demonstrates a substantial accuracy in estimating the continuum limit, using only a few simulations of relatively small system sizes. This efficiency significantly reduces computational cost, providing a powerful tool to facilitate computations for more complex and challenging systems, such as molecules and crystals modeled using ab initio methods.

Rate Constants of the H + HCF3 → H2 + CF3 Reaction from Ring Polymer Molecular Dynamics on a Highly Accurate Potential Energy Surface.

The reaction between H and HCF3 is the primary consumption pathway of HCF3 in the atmosphere and combustion. In this work, ring polymer molecular dynamics (RPMD) calculations are performed to calculate the rate constants of the reaction on a recently developed accurate potential energy surface. 36, 20, and 8 beads are used to compute the rate constants at 350 K ≤ T < 800 K, 800 K ≤ T ≤ 1000 K, and T > 1000 K, respectively. The obtained RPMD rate constants agree well with the experimental measurements. In addition, a detailed analysis of the free-energy curves and transmission coefficients reveals that the quantum tunneling significantly affects the reaction dynamics, even at high temperatures.

The Hydrogen-Bond Continuum in the Salt/Cocrystal Systems of Quinoline and Chloro-Nitrobenzoic Acids.

This study investigates the hydrogen-bond geometry in six two-component solid systems composed of quinoline and chloro-nitrobenzoic acids. New X-ray diffraction studies were conducted using both the conventional independent-atom model and the more recent Hirshfeld atom-refinement method, with the latter providing precise hydrogen-atom positions. The systems can be divided into salts (the hydrogen atom transferred to the quinoline nitrogen), cocrystals (the hydrogen atom retained by the acid), and intermediate structures. Solid-state NMR experiments corroborated the X-ray diffraction-derived H-N distances. DFT calculations, using five functionals including hybrid B3LYP and PBE0, showed varying energy profiles for the hydrogen bonds, with notable differences across functionals. These calculations revealed different preferences for salt or cocrystal structures, depending on the functional used. Path-integral molecular dynamics simulations incorporating nuclear quantum effects demonstrated significant hydrogen-atom delocalization, forming a hydrogen-bond continuum, and provided average N-H distances in excellent agreement with experimental results. This comprehensive experimental and theoretical approach highlights the complexity of multicomponent solids. The study emphasizes that the classification into salts or cocrystals is frequently inadequate, as the hydrogen atom is often significantly delocalized in the hydrogen bond. This insight is crucial for understanding and predicting the behavior of such systems in pharmaceutical applications.

Application of Deconvolution in Path Integral Simulations.

In path integral molecular dynamics (PIMD) simulations, atoms are represented by several replicas connected with harmonic springs, so additional vibrations appear beyond the physical vibrations because of the normal mode frequencies coming from the springs of the ring polymer. In harmonic approximation, the frequencies of these internal modes can be determined exactly from the physical frequencies. We show that this formal effect of the path integral simulations on the vibrations can be considered as a convolution if we use the square of the frequency as an independent variable. This convolution can be represented as a matrix multiplication. The potential of the formalism is demonstrated in two applications. We present an alternative method to determine the power spectrum of thermostats used in PIMD simulations. We also show that in simple anharmonic model systems, the physical frequencies can be obtained from ring polymer molecular dynamics simulations by deconvolution, even in cases where spurious resonances appear.

Deuteration removes quantum dipolar defects from KDP crystals

Dielectric properties of the hydrogen-bonded ferroelectric crystal KH2PO4 (KDP) differ significantly from those of KD2PO4 (DKDP). It is well established that deuteration affects the interplay of hydrogen-bond switches and heavy ion displacements that underlie the emergence of macroscopic polarization, but a detailed microscopic model is missing. We show that all-atom path integral molecular dynamics simulations can predict the isotope effects, revealing the microscopic mechanism that differentiates KDP and DKDP. Proton tunneling generates phosphate configurations that do not contribute to the polarization. At low temperatures, these quantum dipolar defects are substantial in KDP but negligible in DKDP. These intrinsic defects explain why KDP has lower spontaneous polarization and transition entropy than DKDP. The prominent role of quantum fluctuations in KDP is related to the unusual strength of the hydrogen bonds and should be equally important in other crystals of the KDP family, which exhibit similar isotope effects.

How to simulate dissociative chemisorption of methane on metal surfaces.

The dissociation of methane is not only an important reaction step in catalytic processes, but also of fundamental interest. Dynamical effects during the dissociative chemisorption of methane on metal surfaces cause significant differences in computed reaction rates, compared to what is predicted by typical transition state theory (TST) models. It is clear that for a good understanding of the catalytic activation of methane dynamical simulations are required. In this paper, a general blueprint is provided for performing dynamical simulations of the dissociative chemisorption of methane on metal surfaces, by employing either the quasi-classical trajectory or ring polymer molecular dynamics approach. If the computational setup is constructed with great care-since results can be affected considerably by the setup - chemically accurate predictions are achievable. Although this paper concerns methane dissociation, the provided blueprint is, so far, applicable to the dissociative chemisorption of most molecules.