Path Integral Molecular Dynamics Simulations Research Articles

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.

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.

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.

Quantum nuclear motion in silicene: Assessing structural and vibrational properties through path-integral simulations

This paper explores the interplay between quantum nuclear motion and anharmonicity, which causes nontrivial effects on the structural and dynamical characteristics of silicene, a two-dimensional (2D) allotrope of silicon with interesting electronic and mechanical properties. Employing path-integral molecular dynamics (PIMD) simulations, we investigate the quantum delocalization of nuclei, unraveling its impact on the behavior of silicene at the atomic scale. Our study reveals that this delocalization induces significant deviations in the structural parameters of silicene, influencing in-plane surface area, bond lengths, angles, compressibility, and overall lattice dynamics. Through extensive simulations, we delve into the temperature-dependent behavior between 25 and 1200 K, unveiling the role of quantum nuclear fluctuations in dictating thermal expansion and phonon spectra. The extent of nuclear quantum effects is assessed by comparing results of PIMD simulations using an efficient tight-binding Hamiltonian, with those obtained from classical molecular dynamics simulations. The observed quantum effects showcase non-negligible deviations from classical predictions, emphasizing the need for accurate quantum treatments in understanding the material’s behavior at finite temperatures. At low T, the 2D compression modulus of silicene decreases by a 14% due to quantum nuclear motion. We compare the magnitude of quantum effects in this material with those in other related 2D crystalline solids, such as graphene and SiC monolayers.

Generating interstitial water within the persisting tetrahedral H-bond network explains density increase upon compressing liquid water

Despite its ubiquitous nature, the atomic structure of water in its liquid state is still controversially debated. We use a combination of X-ray Raman scattering spectroscopy in conjunction with ab initio and path integral molecular dynamics simulations to study the local atomic and electronic structure of water under high pressure conditions. Systematically increasing fingerprints of non-hydrogen-bonded H[Formula: see text]O molecules in the first hydration shell are identified in the experimental and computational oxygen K-edge excitation spectra. This provides evidence for a compaction mechanism in terms of a continuous collapse of the second hydration shell with increasing pressure via generation of interstitial water within locally tetrahedral hydrogen-bonding environments.

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.

Two liquid states of distinguishable helium-4: The existence of another non-superfluid frozen by heating.

We show that two liquid states can exist in distinguishable helium-4 (4He) obeying Boltzmann statistics by path integral centroid molecular dynamics (CMD) simulations. This is an indication of quantum liquid polyamorphism induced by the nuclear quantum effect. For 0.08-3.3K and 1-500bar, we extensively conducted the isothermal-isobaric CMD simulations to explore not only possible states and state diagram but also the state characteristics. The distinguishable 4He below 25bar does not freeze down to 0.1K even though it includes no Bosonic exchange effect and, therefore, no Bose condensation. One liquid state, low quantum-dispersion liquid (LQDL), is nearly identical to normal liquid He-I of real 4He. The other is high quantum-dispersion liquid (HQDL) consisting of atoms with longer quantum wavelength. This is another non-superfluid existing below 0.5K or the temperatures of LQDL. The HQDL is also a low-entropy and fragile liquid to exhibit, unlike conventional liquids, rather gas-like relaxation of velocity autocorrelation function, while there the atoms diffuse without noticeable contribution from quantum tunneling. The LQDL-HQDL transition is not a thermodynamic phase transition but a continuous crossover accompanied by the change in the expansion factor of quantum wavelength. Freezing of HQDL into the low quantum-dispersion amorphous solid occurs by heating from 0.2 to 0.3K at 40-50bar, while this P-T condition coincides with the Kim-Chan normal-supersolid phase boundary of real 4He. The obtained state diagram was compared to that of the confined subnano-scale 4He systems, where Bosonic correlation is considerably suppressed.

Nuclear quantum effects on the vibrational dynamics of the water-air interface.

We have applied path-integral molecular dynamics simulations to investigate the impact of nuclear quantum effects on the vibrational dynamics of water molecules at the water-air interface. The instantaneous fluctuations in the frequencies of the O-H stretch modes are calculated using the wavelet method of time series analysis, while the time scales of vibrational spectral diffusion are determined from frequency-time correlation functions and joint probability distributions. We find that the inclusion of nuclear quantum effects leads not only to a redshift in the vibrational frequency distribution by about 120 cm-1 for both the bulk and interfacial water molecules but also to an acceleration of the vibrational dynamics at the water-air interface by as much as 35%. In addition, a blueshift of about 45 cm-1 is seen in the vibrational frequency distribution of interfacial water molecules compared to that of the bulk. Furthermore, the dynamics of water molecules beyond the topmost molecular layer was found to be rather similar to that of bulk water.

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.

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.

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.

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.