Path Integral Molecular Dynamics Simulations Research Articles

Neutron scattering and neural-network quantum molecular dynamics investigation of the vibrations of ammonia along the solid-to-liquid transition

Vibrational spectroscopy allows us to understand complex physical and chemical interactions of molecular crystals and liquids such as ammonia, which has recently emerged as a strong hydrogen fuel candidate to support a sustainable society. We report inelastic neutron scattering measurement of vibrational properties of ammonia along the solid-to-liquid phase transition with high enough resolution for direct comparisons to ab-initio simulations. Theoretical analysis reveals the essential role of nuclear quantum effects (NQEs) for correctly describing the intermolecular spectrum as well as high energy intramolecular N-H stretching modes. This is achieved by training neural network models using ab-initio path-integral molecular dynamics (PIMD) simulations, thereby encompassing large spatiotemporal trajectories required to resolve low energy dynamics while retaining NQEs. Our results not only establish the role of NQEs in ammonia but also provide general computational frameworks to study complex molecular systems with NQEs.

Water at negative pressure: nuclear quantum effects

Various condensed phases of water, spanning from the liquid state to multiple ice phases, have been systematically investigated under extreme conditions of pressure and temperature to delineate their stability boundaries. This study focuses on probing the mechanical stability of liquid water through path-integral molecular dynamics simulations, employing the q-TIP4P/F potential to model interatomic interactions in flexible water molecules. Temperature and pressure conditions ranging from 250 to 375 K and − 0.3 to 1 GPa, respectively, are considered. This comprehensive approach enables a thorough exploration of nuclear quantum effects on various physical properties of water through direct comparisons with classical molecular dynamics results employing the same potential model. Key properties such as molar volume, intramolecular bond length, H–O–H angle, internal and kinetic energy are analysed, with a specific focus on the effect of tensile stress. Particular attention is devoted to the liquid-gas spinodal pressure, representing the limit of mechanical stability for the liquid phase, at several temperatures. The quantum simulations reveal a spinodal pressure for water of − 286 and − 236 MPa at temperatures of 250 and 300 K, respectively. At these temperatures, the discernible shifts induced by nuclear quantum motion are quantified at 15 and 10 MPa, respectively. These findings contribute valuable insights into the interplay of quantum effects on the stability of liquid water under diverse thermodynamic conditions.

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.

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.

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.

Quadratic scaling bosonic path integral molecular dynamics.

Bosonic exchange symmetry leads to fascinating quantum phenomena, from exciton condensation in quantum materials to the superfluidity of liquid 4He. Unfortunately, path integral molecular dynamics (PIMD) simulations of bosons are computationally prohibitive beyond ∼100 particles, due to a cubic scaling with the system size. We present an algorithm that reduces the complexity from cubic to quadratic, allowing the first simulations of thousands of bosons using PIMD. Our method is orders of magnitude faster, with a speedup that scales linearly with the number of particles and the number of imaginary time slices (beads). Simulations that would have otherwise taken decades can now be done in days. In practice, the new algorithm eliminates most of the added computational cost of including bosonic exchange effects, making them almost as accessible as PIMD simulations of distinguishable particles.

A path integral molecular dynamics study on the muoniated xanthene-thione molecule.

A positive Mu is a useful tool for investigating the spin density of radical species. The theoretical estimation of its behavior in a molecule requires the inclusion of a quantum effect due to the small mass of muonium. Herein, we performed abinitio a path integral molecular dynamics (PIMD) simulation, which accurately included a multi-dimensional quantum effect, for muoniated 9H-xanthene-9-thione (μXT). Our results showed that the quantum effect significantly increased the hyperfine coupling constant (HFCC) value of μXT, which qualitatively improved the calculated HFCC value, compared to the experimental one. In the PIMD simulation, the bond length between muonium and sulfur in μXT is longer than that between hydrogen and sulfur in a hydrogenated 9H-xanthene-9-thione (HXT), leading to a spin density transfer from XT (9H-xanthene-9-thione) to muonium due to neutral dissociations. Additionally, we found that the S-Mu bond in μXT prefers a structure perpendicular to the molecular plane, where the interaction between Mu and the singly occupied molecular orbital of μXT is the strongest. These structural changes resulted in a larger HFCC value in the PIMD simulation of μXT.

Intra-cluster Charge Migration upon Hydration of Protonated Formic Acid Revealed by Anharmonic Analysis of Cold Ion Vibrational Spectra.

The rates of many chemical reactions are accelerated when carried out in micron-sized droplets, but the molecular origin of the rate acceleration remains unclear. One example is the condensation reaction of 1,2-diaminobenzene with formic acid to yield benzimidazole. The observed rate enhancements have been rationalized by invoking enhanced acidity at the surface of methanol solvent droplets with low water content to enable protonation of formic acid to generate a cationic species (protonated formic acid or PFA) formed by attachment of a proton to the neutral acid. Because PFA is the key feature in this reaction mechanism, vibrational spectra of cryogenically cooled, microhydrated PFA·(H2O)n=1-6 were acquired to determine how the extent of charge localization depends on the degree of hydration. Analysis of these highly anharmonic spectra with path integral ab initio molecular dynamics simulations reveals the gradual displacement of the excess proton onto the water network in the microhydration regime at low temperatures with n = 3 as the tipping point for intra-cluster proton transfer.

Structural and Dynamical Properties of H2O and D2O under Confinement.

Water (H2O) is of great societal importance, and there has been a significant amount of research on its fundamental properties and related physical phenomena. Deuterium dioxide (D2O), known as heavy water, also draws much interest as an important medium for medical imaging, nuclear reactors, etc. Although many experimental studies on the fundamental properties of H2O and D2O have been conducted, they have been primarily limited to understanding the differences between H2O and D2O in the bulk state. In this paper, using path integral molecular dynamics simulations, the structural and dynamical properties of H2O and D2O in bulk and under nanoscale confinement in a (14,0) carbon nanotube are studied. We find that in bulk, structural properties such as bond angle and bond length of D2O are slightly smaller than those of H2O while D2O is slightly more structured than H2O. The dipole moment of D2O tends to be 4% higher than that of H2O, and the hydrogen bonding of D2O is also stronger than that of H2O. Under nanoscale confinement in a (14,0) carbon nanotube, H2O and D2O exhibit a smaller bond length and bond angle. The hydrogen bond number decreases, which demonstrates a weakened hydrogen bond interaction. Moreover, confinement results in a lower libration frequency and a higher OH(OD) bond stretching frequency with an almost unchanged HOH(DOD) bending frequency. The D2O-filled (14,0) carbon nanotube is found to have a smaller radial breathing mode than the H2O-filled (14,0) carbon nanotube.

Many-body quantum muon effects and quadrupolar coupling in solids

Strong quantum zero-point motion (ZPM) of light nuclei and other particles is a crucial aspect of many state-of-the-art quantum materials. However, it has only recently begun to be explored from an ab initio perspective, through several competing approximations. Here we develop a unified description of muon and light nucleus ZPM and establish the regimes of anharmonicity and positional quantum entanglement where different approximation schemes apply. Via density functional theory and path-integral molecular dynamics simulations we demonstrate that in solid nitrogen, α–N2, muon ZPM is both strongly anharmonic and many-body in character, with the muon forming an extended electric-dipole polaron around a central, quantum-entangled [N2–μ–N2]+ complex. By combining this quantitative description of quantum muon ZPM with precision muon quadrupolar level-crossing resonance experiments, we independently determine the static 14N nuclear quadrupolar coupling constant of pristine α–N2 to be –5.36(2) MHz, a significant improvement in accuracy over the previously-accepted value of –5.39(5) MHz, and a validation of our unified description of light-particle ZPM.

Effect of Quantum Delocalization on Temperature Dependent Double Proton Transfer in Molecular Crystals of Terephthalic Acid.

Double proton transfers (DPTs) are important for several physical processes, both in molecules and in the condensed phase. While these have been widely studied in biological systems, their study in crystalline environments is rare. In this work, using path integral molecular dynamics simulations, we have studied temperature dependent DPT in molecular crystals of terephthalic acid (TPA). In accordance with experimental reports, we find evidence for a double proton transfer induced order-to-disorder transition that is sensitive to the inclusion of nuclear quantum effects. Our simulations show that in addition to the presence of L and R tautomers of terepthalic acid, there are a small but non-negligible concentration of positive and negatively charged pairs of TPA molecules. At the onset of the transition at low temperatures, DPT likely occurs through a tunneling mechanism while at room temperature, likely involving the dominance of activated hopping. Through an analysis of the electronic structure of the system using Wannier functions, we show that the H atom shuttling between the donor and acceptor O atoms involves a proton.

The Importance of Nuclear Quantum Effects on the Thermodynamic and Structural Properties of Low-Density Amorphous Ice: A Comparison with Hexagonal Ice.

We study the nuclear quantum effects (NQE) on the thermodynamic properties of low-density amorphous ice (LDA) and hexagonal ice (Ih) at P = 0.1 MPa and T ≥ 25 K. Our results are based on path-integral molecular dynamics (PIMD) and classical MD simulations of H2O and D2O using the q-TIP4P/F water model. We show that the inclusion of NQE is necessary to reproduce the experimental properties of LDA and ice Ih. While MD simulations (no NQE) predict that the density ρ(T) of LDA and ice Ih increases monotonically upon cooling, PIMD simulations indicate the presence of a density maximum in LDA and ice Ih. MD and PIMD simulations also predict a qualitatively different T-dependence for the thermal expansion coefficient αP(T) and bulk modulus B(T) of both LDA and ice Ih. Remarkably, the ρ(T), αP(T), and B(T) of LDA are practically identical to those of ice Ih. The origin of the observed NQE is due to the delocalization of the H atoms, which is identical in LDA and ice Ih. H atoms delocalize considerably (over a distance ≈ 20-25% of the OH covalent-bond length) and anisotropically (preferentially perpendicular to the OH covalent bond), leading to less linear hydrogen bonds HB (larger HOO angles and longer OO separations) than observed in classical MD simulations.

Enabling large-scale quantum path integral molecular dynamics simulations through the integration of Dcdftbmd and i-PI codes.

A large-scale quantum chemical calculation program, Dcdftbmd, was integrated with a Python-based advanced atomistic simulation program, i-PI. The implementation of a client-server model enabled hierarchical parallelization with respect to replicas and force evaluations. The established framework demonstrated that quantum path integral molecular dynamics simulations can be executed with high efficiency for systems consisting of a few tens of replicas and containing thousands of atoms. The application of the framework to bulk water systems, with and without an excess proton, demonstrated that nuclear quantum effects are significant for intra- and inter-molecular structural properties, including oxygen-hydrogen bond distance and radial distribution function around the hydrated excess proton.

Mechanism of Cations Suppressing Proton Diffusion Kinetics for Electrocatalysis.

Proton transfer is crucial for electrocatalysis. Accumulating cations at electrochemical interfaces can alter the proton transfer rate and then tune electrocatalytic performance. However, the mechanism for regulating proton transfer remains ambiguous. Here, we quantify the cation effect on proton diffusion in solution by hydrogen evolution on microelectrodes, revealing the rate can be suppressed by more than 10 times. Different from the prevalent opinions that proton transport is slowed down by modified electric field, we found water structure imposes a more evident effect on kinetics. FTIR test and path integral molecular dynamics simulation indicate that proton prefers to wander within the hydration shell of cations rather than to hop rapidly along water wires. Low connectivity of water networks disrupted by cations corrupts the fast-moving path in bulk water. This study highlights the promising way for regulating proton kinetics via a modified water structure.

Mechanism of Cations Suppressing Proton Diffusion Kinetics for Electrocatalysis

AbstractProton transfer is crucial for electrocatalysis. Accumulating cations at electrochemical interfaces can alter the proton transfer rate and then tune electrocatalytic performance. However, the mechanism for regulating proton transfer remains ambiguous. Here, we quantify the cation effect on proton diffusion in solution by hydrogen evolution on microelectrodes, revealing the rate can be suppressed by more than 10 times. Different from the prevalent opinions that proton transport is slowed down by modified electric field, we found water structure imposes a more evident effect on kinetics. FTIR test and path integral molecular dynamics simulation indicate that proton prefers to wander within the hydration shell of cations rather than to hop rapidly along water wires. Low connectivity of water networks disrupted by cations corrupts the fast‐moving path in bulk water. This study highlights the promising way for regulating proton kinetics via a modified water structure.

Nuclear quantum effects on zeolite proton hopping kinetics explored with machine learning potentials and path integral molecular dynamics

Proton hopping is a key reactive process within zeolite catalysis. However, the accurate determination of its kinetics poses major challenges both for theoreticians and experimentalists. Nuclear quantum effects (NQEs) are known to influence the structure and dynamics of protons, but their rigorous inclusion through the path integral molecular dynamics (PIMD) formalism was so far beyond reach for zeolite catalyzed processes due to the excessive computational cost of evaluating all forces and energies at the Density Functional Theory (DFT) level. Herein, we overcome this limitation by training first a reactive machine learning potential (MLP) that can reproduce with high fidelity the DFT potential energy surface of proton hopping around the first Al coordination sphere in the H-CHA zeolite. The MLP offers an immense computational speedup, enabling us to derive accurate reaction kinetics beyond standard transition state theory for the proton hopping reaction. Overall, more than 0.6 μs of simulation time was needed, which is far beyond reach of any standard DFT approach. NQEs are found to significantly impact the proton hopping kinetics up to ~473 K. Moreover, PIMD simulations with deuterium can be performed without any additional training to compute kinetic isotope effects over a broad range of temperatures.

Collective Proton Transfers in Cyclic Water–Ammonia Tetramers: A Path Integral Machine-Learning Study

We present results from machine-learning-based path integral molecular dynamics simulations that describe isomerization paths articulated via collective proton transfers along mixed, cyclic tetramers combining water and ammonia at cryogenic conditions. The net result of such isomerizations is a reverse of the chirality of the global hydrogen-bonding architecture along the different cyclic moieties. In monocomponent tetramers, the classical free energy profiles associated with these isomerizations present the usual symmetric double-well characteristics whereas the reactive paths exhibit full concertedness among the different intermolecular transfer processes. Contrastingly, in mixed water/ammonia tetramers, the incorporation of a second component introduces imbalances in the strengths of the different hydrogen bonds leading to a partial loss of concertedness, most notably at the vicinity of the transition state. As such, the highest and lowest degrees of progression are registered along OH···N and O···HN coordinations, respectively. These characteristics lead to polarized transition state scenarios akin to solvent-separated ion-pair configurations. The explicit incorporation of nuclear quantum effects promotes drastic depletions in the activation free energies and modifications in the overall shape of the profiles which include central plateau-like stages, indicating the prevalence of deep tunneling regimes. On the other hand, the quantum treatment of the nuclei partially restores the degree of concertedness among the evolutions of the individual transfers.

Path-integral molecular dynamics predictions of equilibrium H and O isotope fractionations between brucite and water

Precise prediction of hydrogen (H) isotope fractionations among different substances is a long-standing challenge in isotope geochemistry, as it needs treatments beyond the harmonic approximation. Path-integral molecular dynamics (PIMD) simulations have recently been proved to be valid in predicting equilibrium isotope fractionations for light elements. However, the lack of reliable force fields hinders the application of PIMD to the condensed phases. In this work, the deep potential models trained on the first-principles molecular dynamics (FPMD) data are applied to PIMD simulations to determine the D/H and 18O/16O isotope fractionation factors between brucite and water. To quantitatively assess the influence of quantum effects, the D/H and 18O/16O isotope fractionations between brucite and water are also determined under the framework of harmonic approximation.By comparing the results of the PIMD and those of harmonic calculations, we find that H and O atoms in brucite and water are strongly affected by the anharmonic effects. The accuracy of the harmonic isotope fractionations mainly depends on the cancellation percentage of the anharmonic effects on the reduced partition function ratios (RPFRs) of the two substances. The isotope fractionations predicted by PIMD simulations are close to the experimental results at high temperatures. However, at low temperatures, the predicted isotope fractionations are different from those of experiments. The discrepancies are attributed to the approximations in the density functional theory (DFT) functionals used in this work, i.e., the inaccuracy in predicting the proton positions in brucite and water. Moreover, we find that the H isotope fractionations are extremely sensitive to the change of pressure. The direction of H isotope fractionation between brucite and water (at 300 K) will be even inversed as the pressure increases to more than 3 GPa. This strong pressure sensitivity may be a common characteristic of hydrous minerals and water systems. Therefore, the extrapolated H isotope fractionations used in the studies of dehydration of subducting slabs or deep Earth H distribution based on the isotope fractionations under low pressures may need to be rechecked.

Data-Efficient Machine Learning Potentials from Transfer Learning of Periodic Correlated Electronic Structure Methods: Liquid Water at AFQMC, CCSD, and CCSD(T) Accuracy.

Obtaining the atomistic structure and dynamics of disordered condensed-phase systems from first-principles remains one of the forefront challenges of chemical theory. Here we exploit recent advances in periodic electronic structure and provide a data-efficient approach to obtain machine-learned condensed-phase potential energy surfaces using AFQMC, CCSD, and CCSD(T) from a very small number (≤200) of energies by leveraging a transfer learning scheme starting from lower-tier electronic structure methods. We demonstrate the effectiveness of this approach for liquid water by performing both classical and path integral molecular dynamics simulations on these machine-learned potential energy surfaces. By doing this, we uncover the interplay of dynamical electron correlation and nuclear quantum effects across the entire liquid range of water while providing a general strategy for efficiently utilizing periodic correlated electronic structure methods to explore disordered condensed-phase systems.

Direct Elucidation of the Vibrationally Averaged Structure of Benzene: A Path Integral Molecular Dynamics Study.

Path integral molecular dynamics (PIMD) simulations for C6H6, C6D6, and C6T6 have been carried out to directly estimate the distribution of projected C-H(D,T) bond lengths onto the principal axis plane. The average values of raw C-H(D,T) bond lengths obtained from PIMD simulations are in the order of ⟨RC-H⟩ > ⟨RC-D⟩ > ⟨RC-T⟩ due to the anharmonicity of the potential energy curve. However, the projected C-H(D,T) bond lengths are almost the same as those reported by Hirano et al. [J. Mol. Struct. 2021, 1243, 130537]. Our PIMD simulations directly and strongly support the explanation by Hirano et al. for the experimental observations that almost the same projected C-H(D) bond lengths are found for C6H6 and C6D6. The PIMD simulations also predicted the same projected bond lengths for C6T6 as those of C6H(D)6. In addition to the previous local mode analysis, the present PIMD simulations predicted, for benzene isotopologues, that the vibrationally averaged structure is planar but non-flat.