Path Integral Molecular Dynamics Research Articles

Multiple-time scale integration method based on an interpolated potential energy surface for abinitio path integral molecular dynamics.

A new multiple-time scale integration method is presented that propagates abinitio path integral molecular dynamics (PIMD). This method uses a large time step to generate an approximate geometrical configuration whose energy and gradient are evaluated at the level of an abinitio method, and then, a more precise integration scheme, e.g., the Bulirsch-Stoer method or velocity Verlet integration with a smaller time step, is used to integrate from the previous step using the computationally efficient interpolated potential energy surface constructed from two consecutive points. This method makes the integration of PIMD more efficient and accurate compared with the velocity Verlet integration. A Nosé-Hoover chain thermostat combined with this new multiple-time scale method has good energy conservation even with a large time step, which is usually challenging in velocity Verlet integration for PIMD due to the very small chain mass when a large number of beads are used. The new method is used to calculate infrared spectra and free energy profiles to demonstrate its accuracy and capabilities.

Microscopic self-dynamics in liquid Ne-D_{2} mixtures: Quantum features and itinerant oscillators reexamined.

In this paper, we report the results of a centroid molecular dynamics (CMD) study of the canonical velocity autocorrelation functions (VACFs) in liquid Ne-D_{2} mixtures at a temperature of T=30K and in the full D_{2}-concentration range (0%≤x_{D_{2}}≤100%). This binary system was selected because of its moderate, although sizable, quantum effects which, as far as its equilibrium properties are concerned, are fully described by the path integral Monte Carlo (PIMC) simulations that have been also implemented. A comprehensive test of the VACF spectral moments carried out using three physical quantities (namely, mean kinetic energy, Einstein frequency, and mean-squared force) obtained from PIMC was performed revealing the potentialities, as well as the limitations, of the CMD approach to the single-particle dynamics in these low-T liquid mixtures. Additional physical information was extracted from the canonical VACFs by fitting their spectra via two distinct methods: the Levesque-Verlet model (LV, very flexible but highly heuristic) and the itinerant oscillator model (IO, based on the physical ground of a single particle rattling inside a short-lived diffusing pseudocage). Both provided good fits of the CMD outputs, with LV being always more adequate than IO in the case of the Ne VACFs, while, as for the D_{2} VACFs, the LV superiority is evident only at high x_{D_{2}} values. However, a peculiar and systematic effect was found after analyzing the IO-fitted parameters: the estimated pseudocage masses turned out to be at least one order of magnitude lower than the corresponding values inferred from the PIMC simulations. This outcome concerns both the Ne and the D_{2} rattling molecule and, as we also discovered, had already been observed (but promptly forgotten) in purely classical simulations of liquid Ar. The possible physical origins of this finding have been finally discussed in some detail, also in connection with the result of the more recent exponential expansion theory (EET), which manages to shed more light on the concept of single particles rattling inside short-lived pseudocages, ultimately demonstrating its untenability.

Quantum symmetrization transition in superconducting sulfur hydride from quantum Monte Carlo and path integral molecular dynamics

AbstractWe study the structural phase transition, originally associated with the highest superconducting critical temperature Tc measured in high-pressure sulfur hydride. A quantitative description of its pressure dependence has been elusive for any ab initio theory attempted so far, raising questions on the actual mechanism leading to the maximum of Tc. Here, we estimate the critical pressure of the hydrogen bond symmetrization in the Im$$\bar{3}$$ 3 ¯ m structure, by combining density functional theory and quantum Monte Carlo simulations for electrons with path integral molecular dynamics for quantum nuclei. We find that the Tc maximum corresponds to pressures where local dipole moments dynamically form on the hydrogen sites, as precursors of the ferroelectric Im$$\bar{3}$$ 3 ¯ m-R3m transition, happening at lower pressures. For comparison, we also apply the self-consistent harmonic approximation, whose ferroelectric critical pressure lies in between the ferroelectric transition estimated by path integral molecular dynamics and the local dipole formation. Nuclear quantum effects play a major role in a significant reduction (≈50 GPa) of the classical ferroelectric transition pressure at 200 K and in a large isotope shift (≈25 GPa) upon hydrogen-to-deuterium substitution of the local dipole formation pressure, in agreement with the corresponding change in the Tc maximum location.

Nuclear Quantum Effects Accelerate Charge Recombination but Boost the Stability of Inorganic Perovskites in Mild Humidity.

Experiments have demonstrated that mild humidity can enhance the stability of the CsPbBr3 perovskite, though the underlying mechanism remains unclear. Utilizing ab initio molecular dynamics, ring polymer molecular dynamics, and non-adiabatic molecular dynamics, our study reveals that nuclear quantum effects (NQEs) play a crucial role in stabilizing the lattice rigidity of the perovskite while simultaneously shortening the charge carrier lifetime. NQEs reduce the extent of geometric disorder and the number of atomic fluctuations, diminish the extent of hole localization, and thereby improve the electron-hole overlap and non-adiabatic coupling. Concurrently, these effects significantly suppress phonon modes and slow decoherence. As a result, these factors collectively accelerate charge recombination by a factor of 1.42 compared to that in scenarios excluding NQEs. The resulting sub-10 ns recombination time scales align remarkably well with experimental findings. This research offers novel insight into how moisture resistance impacts the stability and charge carrier lifetime in all-inorganic perovskites.

Location of the Shared Proton in Proton-Bound Dimer Compound of Hydrogen Sulfate and Formate: Path Integral Molecular Dynamics Study.

The structure of the proton-bound dimer compound of hydrogen sulfate and formate has been studied by considering nuclear quantum effects (NQEs) using the path integral molecular dynamics method. This study unveiled the location of the shared proton and answered the following question: "Is the shared proton localized on either an anion or located around the center of two anions?" We have elucidated that the shared proton is distributed in the region beyond the transition state due to the NQEs, even though the shared proton did not completely overcome the transition state for the proton shuttle.

Accurate Reaction Probabilities for Translational Energies on Both Sides of the Barrier of Dissociative Chemisorption on Metal Surfaces.

Molecular dynamics simulations are essential for a better understanding of dissociative chemisorption on metal surfaces, which is often the rate-controlling step in heterogeneous and plasma catalysis. The workhorse quasi-classical trajectory approach ubiquitous in molecular dynamics is able to accurately predict reactivity only for high translational and low vibrational energies. In contrast, catalytically relevant conditions generally involve low translational and elevated vibrational energies. Existing quantum dynamics approaches are intractable or approximate as a result of the large number of degrees of freedom present in molecule-metal surface reactions. Here, we extend a ring polymer molecular dynamics approach to fully include, for the first time, the degrees of freedom of a moving metal surface. With this approach, experimental sticking probabilities for the dissociative chemisorption of methane on Pt(111) are reproduced for a large range of translational and vibrational energies by including nuclear quantum effects and employing full-dimensional simulations.

Modeling Temperature-Dependent Electron Thermal Diffuse Scattering via Machine-Learned Interatomic Potentials and Path-Integral Molecular Dynamics.

Electron thermal diffuse scattering is shown to be sensitive to subtle changes in atomic vibrations and shows promise in assessing lattice dynamics at nanometer resolution. Here, we demonstrate that machine-learned interatomic potentials (MLIPs) and path-integral molecular dynamics can accurately capture the potential energy landscape and lattice dynamics needed to describe electron thermal diffuse scattering. Using SrTiO_{3} as a test bed at cryogenic and room temperatures, we compare electron thermal diffuse scattering simulations using different approximations to incorporate thermal motion. Only when the simulations are based on quantum mechanically accurate MLIPs in combination with path-integral molecular dynamics that include nuclear quantum effects is there excellent agreement with experiments.

Electronic friction near metal surface: Incorporating nuclear quantum effect with ring polymer molecular dynamics.

The molecular dynamics with electronic friction (MDEF) approach can accurately describe nonadiabatic effects at metal surfaces in the weakly nonadiabatic limit. That being said, the MDEF approach treats nuclear motion classically such that the nuclear quantum effects are completely missing in the approach. To address this limitation, we combine Electronic Friction with Ring Polymer Molecular Dynamics (EF-RPMD). In particular, we apply the averaged electronic friction from the metal surface to the centroid mode of the ring polymer. We benchmark our approach against quantum dynamics to show that EF-RPMD can accurately capture zero-point energy as well as transition dynamics. In addition, we show that EF-RPMD can correctly predict the electronic transfer rate near metal surfaces in the tunneling limit as well as the barrier crossing limit. We expect that our approach will be very useful to study nonadiabatic dynamics near metal surfaces when nuclear quantum effects become essential.

Neural Network Corrections to Intermolecular Interaction Terms of a Molecular Force Field Capture Nuclear Quantum Effects in Calculations of Liquid Thermodynamic Properties.

We incorporate nuclear quantum effects (NQE) in condensed matter simulations by introducing short-range neural network (NN) corrections to the ab initio fitted molecular force field ARROW. Force field NN corrections are fitted to average interaction energies and forces of molecular dimers, which are simulated using the Path Integral Molecular Dynamics (PIMD) technique with restrained centroid positions. The NN-corrected force field allows reproduction of the NQE for computed liquid water and methane properties such as density, radial distribution function (RDF), heat of evaporation (HVAP), and solvation free energy. Accounting for NQE through molecular force field corrections circumvents the need for explicit computationally expensive PIMD simulations in accurate calculations of the properties of chemical and biological systems. The accuracy and locality of pairwise NN NQE corrections indicate that this approach could be applicable to complex heterogeneous systems, such as proteins.

Nuclear quantum effects in the acetylene:ammonia plastic co-crystal.

Organic molecular solids can exhibit rich phase diagrams. In addition to structurally unique phases, translational and rotational degrees of freedom can melt at different state points, giving rise to partially disordered solid phases. The structural and dynamic disorder in these materials can have a significant impact on the physical properties of the organic solid, necessitating a thorough understanding of disorder at the atomic scale. When these disordered phases form at low temperatures, especially in crystals with light nuclei, the prediction of material properties can be complicated by the importance of nuclear quantum effects. As an example, we investigate nuclear quantum effects on the structure and dynamics of the orientationally disordered, translationally ordered plastic phase of the acetylene:ammonia (1:1) co-crystal that is expected to exist on the surface of Saturn's moon Titan. Titan's low surface temperature (∼90K) suggests that the quantum mechanical behavior of nuclei may be important in this and other molecular solids in these environments. By using neural network potentials combined with ring polymer molecular dynamics simulations, we show that nuclear quantum effects increase orientational disorder and rotational dynamics within the acetylene:ammonia (1:1) co-crystal by weakening hydrogen bonds. Our results suggest that nuclear quantum effects are important to accurately model molecular solids and their physical properties in low-temperature environments.

Simulations of disordered matter in 3D with the morphological autoregressive protocol (MAP) and convolutional neural networks.

Disordered molecular systems, such as amorphous catalysts, organic thin films, electrolyte solutions, and water, are at the cutting edge of computational exploration at present. Traditional simulations of such systems at length scales relevant to experiments in practice require a compromise between model accuracy and quality of sampling. To address this problem, we have developed an approach based on generative machine learning called the Morphological Autoregressive Protocol (MAP), which provides computational access to mesoscale disordered molecular configurations at linear cost at generation for materials in which structural correlations decay sufficiently rapidly. The algorithm is implemented using an augmented PixelCNN deep learning architecture that, as we previously demonstrated, produces excellent results in 2 dimensions (2D) for mono-elemental molecular systems. Here, we extend our implementation to multi-elemental 3D and demonstrate performance using water as our test system in two scenarios: (1) liquid water and (2) samples conditioned on the presence of pre-selected motifs. We trained the model on small-scale samples of liquid water produced using path-integral molecular dynamics simulations, including nuclear quantum effects under ambient conditions. MAP-generated water configurations are shown to accurately reproduce the properties of the training set and to produce stable trajectories when used as initial conditions in quantum dynamics simulations. We expect our approach to perform equally well on other disordered molecular systems in which structural correlations decay sufficiently fast while offering unique advantages in situations when the disorder is quenched rather than equilibrated.

Rate coefficients for the O + H2 and O + D2 reactions: how well ring polymer molecular dynamics accounts for tunelling.

We present here extensive calculations of the O(3P) + H2 and O(3P) + D2 reaction dynamics spanning the temperature range from 200 K to 2500 K. The calculations have been carried out using fully converged time-independent quantum mechanics (TI QM), quasiclassical trajectories (QCT) and ring polymer molecular dynamics (RPMD) on the two lowest lying adiabatic potential energy surfaces (PESs), 13A' and 13A'', calculated by Zanchet et al. [J. Chem. Phys., 2019, 151, 094307]. TI QM rate coefficients were determined using the cumulative reaction probability formalism on each PES including all of the total angular momenta and the Coriolis coupling and can be considered to be essentially exact within the Born-Oppenheimer approximation. The agreement between the rate coefficients calculated by using QM and RPMD is excellent for the reaction with D2 in almost the whole temperature range. For the reaction with H2, although the agreement is very good above 500 K, the deviations are significant at lower temperatures. In contrast, the QCT calculations largely underestimate the rate coefficients for the two isotopic variants due to their inability to account for tunelling. The differences found in the disagreements between RPMD and QM rate coefficients for the reactions for both the isotopologues are indicative of the ability of the RPMD method to accurately describe systems where tunelling plays a relevant role. Considering that both reactions are dominated by tunelling below 500 K, the present results show that RPMD is a very powerful tool for determining rate coefficients. The present QM rate coefficients calculated on adiabatic PESs slightly underestimate the best global fits of the experimental measurements, which we attribute to the intersystem crossing with the singlet 11A' PES.

Computational study of the post-transition state dynamics for the OH + CH3OH reaction probed by photodetachment of the CH3O-(H2O) anion.

Dissociative photodetachment dynamics simulations were conducted to study the CH3O-(H2O) → CH3O + H2O + e- reaction using classical molecular dynamics (MD) and ring-polymer molecular dynamics (RPMD) techniques on two newly formulated neutral potential energy surfaces (PES1 and PES2) by different research groups. While the dissociation dynamics exhibited similarities between classical MD and RPMD, there were noticeable differences in the fluctuation of probability densities for the internal modes due to nuclear quantum effects. Upon comparison of our findings with experimental data concerning the electron binding energy distribution and photofragment relative energy, it suggests that the potential energy landscapes of PES2 are reasonably precise. The time evolution of occupied vibrational states of the H2O photofragment is presented in this study.

Nuclear quantum and H/D isotope effects on aromaticity: path integral molecular dynamics study.

Aromaticity is an important concept in organic chemistry, and thus, many theoretical and experimental studies have been conducted so far. However, the majority of theoretical studies have concentrated on the aromaticity of the stationary point structures. Herein, the influence of nuclear quantum fluctuation (nuclear quantum effects: NQEs) and thermal fluctuation on the aromaticity of benzene have been analyzed by path integral molecular dynamics (PIMD) simulation. The PIMD simulations revealed that the NQEs affected not only the C-H bonds but also the C-C bonds. The HOMA and NICS calculations demonstrated that the aromaticity decreased due to the NQEs of carbon atoms, attributed to an increase in the contribution from specific vibrational modes strongly correlated with benzene's aromaticity.

Water-air interface revisited by means of path-integral ab initio molecular dynamics.

Although nuclear quantum effects (NQEs) have been considered on bulk liquid water, the impact of these latter on the air-water interface has not yet been reported. Herein, by performing and comparing ab initio molecular dynamics (AIMD) and path integral AIMD (PI-AIMD) simulations, we reveal the impact of NQEs on structural, dynamical and electronic properties as well as IR spectra of the air-water interface at room temperature. NQEs, being able to describe a more accurate proton delocalization in H-bonded system than AIMD, reveal a different structural arrangement and dynamical behaviour of both bulk and interfacial water molecules in comparison to AIMD results. A more de-structured and de-bound water arrangement and coordination are identified when the quantum nature of nuclei are considered for both bulk and interfacial water molecules. Structural properties, such as inter-/intra-molecular bond lengths, coordination numbers and H-bonding angles of bulk and interfacial water molecules here calculated, are affected by NQEs mitigating the overstructured description given by AIMD. Further evidences of an AIMD overstructured description of bulk water are in the computed IR spectra, where an increased absorption peak intensity and an increased strength of the hydrogen-bond network are alleviated by NQEs. In addition, NQEs show a valuable impact on the electronic structure of the air-water interface, reducing the total valence bandwidth and the electronic energy band-gap when passing from bulk to interfacial water. This work proves how NQEs significantly affect properties and features of the air-water interface, that are essential to accurately describe H-bonded interfacial systems.

Nuclear quantum effects on the intramolecular hydrogen bonds in biuret and biguanide.

We focus on the unique aspects of biuret and biguanide, which form six-membered ring structures via intramolecular hydrogen bonds. The proton donor and acceptor atoms differ between biuret and biguanide, leading to varying energy barrier heights for proton transfer. We performed path integral molecular dynamics (PIMD) simulations for biuret and biguanide to investigate the correlation between proton transfer and the degree of the delocalization of π-electrons in the six-membered ring framework structure. The results indicate that the π-electrons in the framework structure are delocalized regardless of the ease of intramolecular proton transfer.

Mode-coupling theory of lattice dynamics for classical and quantum crystals

The dynamical properties of nuclei, carried by the concept of phonon quasiparticles , are central to the field of condensed matter. While the harmonic approximation can reproduce a number of properties observed in real crystals, the inclusion of anharmonicity in lattice dynamics is essential to accurately predict properties such as heat transport or thermal expansion. For highly anharmonic systems, non-perturbative approaches are needed, which result in renormalized theories of lattice dynamics. In this article, we apply the Mori–Zwanzig projector formalism to derive an exact generalized Langevin equation describing the quantum dynamics of nuclei in a crystal. By projecting this equation on quasiparticles in reciprocal space, and with results from linear response theory, we obtain a formulation of vibrational spectra that fully accounts for the anharmonicity. Using a mode-coupling approach, we construct a systematic perturbative expansion in which each new order is built to minimize the following ones. With a truncation to the lowest order, we show how to obtain a set of self-consistent equations that can describe the lineshapes of quasiparticles. The only inputs needed for the resulting set of equations are the static Kubo correlation functions, which can be computed using (fully quantum) path-integral molecular dynamics or approximated with (classical or ab initio) molecular dynamics. We illustrate the theory with an application on fcc 4He, an archetypal quantum crystal with very strong anharmonicity.

Analysis of the structural H/D isotope effect with an excess proton/deuteron in light/heavy water solvent using path integral molecular dynamics simulations

Abstract We analyzed the difference in the structural H/D isotope effect between an excess proton in light water (H-body) and an excess deuteron in heavy water (D-body), including the nuclear quantum effect, using path integral molecular dynamics simulations. We found that the second peak of the H-body is shorter than that of the D-body in the radial distribution function of O*–O, where O* is the oxygen atom of the H3O+/D3O+ fragment. The main reason for this would be the difference in the ratio of the Zundel structure with the sp3-like configuration, where the Zundel structure in the H-body (14.0%) is greater than that in the D-body (12.0%). We also found rare occurrences of double H3O+/D3O+ configurations, mainly including Zundel–Zundel-like structures such as H7O3+/D7O3+ and H9O4+/D9O4+. The ratios of such configurations appearing in our simulations are 0.89% and 0.20% for the H-body and the D-body, respectively.

Exploring the impact of ions on oxygen K-edge X-ray absorption spectroscopy in NaCl solution using the GW-Bethe-Salpeter-equation approach.

X-ray absorption spectroscopy (XAS) is a powerful experimental tool to probe the local structure in materials with the core hole excitations. Here, the oxygen K-edge XAS spectra of the NaCl solution and pure water are computed by using a recently developed GW-Bethe-Salpeter equation approach, based on configurations modeled by path-integral molecular dynamics with the deep-learning technique. The neural network is trained on abinitio data obtained with strongly constrained and appropriately normed density functional theory. The observed changes in the XAS features of the NaCl solution, compared to those of pure water, are in good agreement between experimental and theoretical results. We provided detailed explanations for these spectral changes that occur when NaCl is solvated in pure water. Specifically, the presence of solvating ion pairs leads to localization of electron-hole excitons. Our theoretical XAS results support the theory that the effects of the solvating ions on the H-bond network are mainly confined within the first hydration shell of ions, however beyond the shell the arrangement of water molecules remains to be comparable to that observed in pure water.

Thermal dependence of the hydrated proton and optimal proton transfer in the protonated water hexamer

Water is a key ingredient for life and plays a central role as solvent in many biochemical reactions. However, the intrinsically quantum nature of the hydrogen nucleus, revealing itself in a large variety of physical manifestations, including proton transfer, gives rise to unexpected phenomena whose description is still elusive. Here we study, by a combination of state-of-the-art quantum Monte Carlo methods and path-integral molecular dynamics, the structure and hydrogen-bond dynamics of the protonated water hexamer, the fundamental unit for the hydrated proton. We report a remarkably low thermal expansion of the hydrogen bond from zero temperature up to 300 K, owing to the presence of short-Zundel configurations, characterised by proton delocalisation and favoured by the synergy of nuclear quantum effects and thermal activation. The hydrogen bond strength progressively weakens above 300 K, when localised Eigen-like configurations become relevant. Our analysis, supported by the instanton statistics of shuttling protons, reveals that the near-room-temperature range from 250 K to 300 K is optimal for proton transfer in the protonated water hexamer.