Path Integral Simulations Research Articles

Path Integral Molecular Dynamics of Liquid Water in a Mean-Field Particle Reservoir.

We present a simulation scheme for path integral simulation of molecular liquids where a small open region is embedded in a large reservoir of non interacting point‐particles. The scheme is based on the latest development of the adaptive resolution technique AdResS and allows for the space‐dependent change of molecular resolution from a path integral representation with 120 degrees of freedom to a point particle that does not interact with other molecules and vice versa. The method is applied to liquid water and implies a sizable gain regarding the request of computational resources compared to full path integral simulations. Given the role of water as universal solvent with a specific hydrogen bonding network, the path integral treatment of water molecules is important to describe the quantum effects of hydrogen atoms’ delocalization in space on the hydrogen bonding network. The method presented here implies feasible computational efforts compared to full path integral simulations of liquid water which, on large scales, are often prohibitive.

Diamagnetic susceptibility from a nonadiabatic path-integral simulation of few-electron systems

Diamagnetism is the response of dynamical compositions of charged particles, electrons, and nuclei, to an incident magnetic field. In this paper, we study how the finite temperature and finite nuclear masses affect the diamagnetic susceptibilities of selected small atoms and molecules, as limiting cases of dilute gas. We use nonrelativistic path-integral Monte Carlo simulation (PIMC), where both electrons and nuclei are treated on equal footing at finite temperatures in sampling exact Coulomb pair density matrices. The PIMC estimator of diamagnetic susceptibility has been briefly introduced in Ceperley [D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995)], but here we present a comprehensive derivation, discussion of practical effects, and proof-of-concept results for selected few-body Coulomb systems. Our results are in perfect agreement with high-accuracy literature references, where available, but also demonstrate additional thermal effects of the diamagnetic response of low-mass systems.

Nuclear Quantum Effects on the Charge-Density Wave Transition in NbX2 (X = S, Se).

Understanding the origin of charge-density wave (CDW) instability is important for manipulating novel collective electronic states. Many layered transition metal dichalcogenides (TMDs) share similarity in the structural and electronic instability, giving rise to diverse CDW phases and superconductivity. It is still puzzling that even isostructural and isoelectronic TMDs show distinct CDW features. For instance, bulk NbSe2 exhibits CDW order at low temperature, while bulk NbS2 displays no CDW instability. The CDW transitions in single-layer NbS2 and NbSe2 are also different. In the classic limit, we investigate the electron correlation effects on the dimensionality dependence of the CDW ordering. By performing ab initio path integral molecular dynamics simulations and comparative analyses, we further revealed significant nuclear quantum effects in these systems. Specifically, the quantum motion of sulfur anions significantly reduces the CDW transition temperature in both bulk and single-layer NbS2, resulting in distinct CDW features in the NbS2 and NbSe2 systems.

Using Machine Learning to Greatly Accelerate Path Integral Ab Initio Molecular Dynamics.

Ab initio molecular dynamics (AIMD) has become one of the most popular and robust approaches for modeling complicated chemical, liquid, and material systems. However, the formidable computational cost often limits its widespread application in simulations of the largest-scale systems. The situation becomes even more severe in cases where the hydrogen nuclei may be better described as quantized particles using a path integral representation. Here, we present a computational approach that combines machine learning with recent advances in path integral contraction schemes, and we achieve a 2 orders of magnitude acceleration over direct path integral AIMD simulation while at the same time maintaining its accuracy.

Path integral simulations of confined parahydrogen molecules within clathrate hydrates: Merging low temperature dynamics with the zero-temperature limit.

Clathrate hydrates, or cages comprised solely of water molecules, have long been investigated as a clean storage facility for hydrogen molecules. A breakthrough occurred when hydrogen molecules were experimentally placed within a structure-II clathrate hydrate, which sparked much interest to determine their feasibility for energy storage [Mao et al., Science 297, 2247-2249 (2002)]. We use Path Integral Molecular Dynamics (PIMD) and Langevin equation Path Integral Ground State (LePIGS) for finite temperature and zero-temperature studies, respectively, to determine parahydrogen occupancy properties in the small dodecahedral (512) and large hexakaidecahedral (51264) sized cages that comprise the structure-II unit cell. We look at energetic and structural properties of small clusters of hydrogen, treated as point-like particles, confined within each of the different sized clathrates, and treated as rigid, to determine energetic and structural properties in the zero-temperature limit. Our predicted hydrogen occupancy within these two cage sizes is consistent with previous literature values. We then calculate the energies as a function of temperature and merge the low temperature results calculated using finite temperature PIMD with the zero-temperature results using LePIGS, demonstrating that the two methods are compatible.

Structures of liquid and aqueous water isotopologues at ambient temperature from ab initio path integral simulations

The heavy hydrogen isotopes D and T are found in trace amounts in water, and when their concentration increases they can play an intricate role in modulating the physical properties of the liquid. We present an analysis of the microscopic structures of ambient light water (H2O(l)), heavy water (D2O(l)), T2O(l), HDO(aq) and HTO(aq) studied by ab initio path integral molecular dynamics (PIMD). Unlike previous ab initio PIMD investigations of H2O(l) and D2O(l) [Chen et al., Phys. Rev. Lett., 2003, 91, 215503] [Machida et al., J. Chem. Phys., 2017, 148, 102324] we find that D2O(l) is more structured than H2O(l), as is predicted by the experiment. The agreement between the experiment and our simulation for H2O(l) and D2O(l) allows us to accurately predict the intra- and intermolecular structures of T2O(l) HDO(aq) and HTO(aq). T2O(l) is found to have a similar intermolecular structure to that of D2O(l), while the intramolecular structure is more compact, giving rise to a smaller dipole moment than those of H2O(l) and D2O(l). For the mixed isotope species, HDO(aq) and HTO(aq), we find smaller dipole moments and fewer hydrogen bonds when compared with the pure species H2O and D2O. We can attribute this effect to the relative compactness of the mixed isotope species, which results in a lower dipole moment than that of the pure species.

Deuterium Isotope Effects on Acid-Base Equilibrium of Organic Compounds.

Deuterium isotope effects on acid–base equilibrium have been investigated using a combined path integral and free-energy perturbation simulation method. To understand the origin of the linear free-energy relationship of versus , we examined two theoretical models for computing the deuterium isotope effects. In Model 1, only the intrinsic isotope exchange effect of the acid itself in water was included by replacing the titratable protons with deuterons. Here, the dominant contribution is due to the difference in zero-point energy between the two isotopologues. In Model 2, the medium isotope effects are considered, in which the free energy change as a result of replacing H2O by D2O in solute–solvent hydrogen-bonding complexes is determined. Although the average change from Model 1 was found to be in reasonable agreement with the experimental average result, the dependence of the solvent isotope effects is absent. A linear free-energy relationship is obtained by including the medium effect in Model 2, and the main factor is due to solvent isotope effects in the anion–water complexes. The present study highlights the significant roles of both the intrinsic isotope exchange effect and the medium solvent isotope effect.

Fast quasi-centroid molecular dynamics.

We describe a fast implementation of the quasi-centroid molecular dynamics (QCMD) method in which the quasi-centroid potential of mean force is approximated as a separable correction to the classical interaction potential. This correction is obtained by first calculating quasi-centroid radial and angular distribution functions in a short path integral molecular dynamics simulation and then using iterative Boltzmann inversion to obtain an effective classical potential that reproduces these distribution functions in a classical NVT simulation. We illustrate this approach with example applications to the vibrational spectra of gas phase molecules, obtaining excellent agreement with QCMD reference calculations for water and ammonia and good agreement with the quantum mechanical vibrational spectrum of methane.

Ab initio study of nuclear quantum effects on sub- and supercritical water.

The structures of water in the ambient, subcritical, and supercritical conditions at various densities were studied systematically by ab initio path integral molecular dynamics simulations. It was found that the nuclear quantum effects (NQEs) have a significant impact on the structure of hydrogen bonds in close contact, not only in the ambient condition but also in the sub- and supercritical conditions. The NQEs on the structure beyond the hydrogen bond contact are important in ambient water, but not much for water in the sub- and supercritical conditions. The NQEs are furthermore important for determining the number of hydrogen bonds in the ambient conditions, and this role is, however, diminished in the sub- and supercritical conditions. The NQEs do, nevertheless, show their importance in determining the intramolecular structure of water and the close contact structures of the hydrogen bonds, even at sub- and supercritical conditions. Using the RPBE-D3 functional, the computed radial distribution functions for ambient water are in excellent agreement with experimental data, upgrading our previous results using the BLYP-D2 functional [Machida et al., J. Chem. Phys. 148, 102324 (2018)]. The computed radial distribution functions for water in the sub- and supercritical conditions were carefully compared with experiment. In particular, we found that the first peak in hydrogen pair distribution functions matches only when the NQEs are taken into account.

Quasiparticle Nature of the Bose Polaron at Finite Temperature.

The Bose polaron has attracted theoretical and experimental interest because the mobile impurity is surrounded by a bath that undergoes a superfluid-to-normal phase transition. Although many theoretical works have studied this system in its ground state, only a few analyze its behavior at finite temperature. We have studied the effect of temperature on a Bose polaron system performing abinitio path integral MonteCarlo simulations. This method is able to approach the critical temperature without losing accuracy, in contrast with perturbative approximations. We have calculated the polaron energy for the repulsive and attractive branches and we have observed an asymmetric behavior between the two branches. When the potential is repulsive, the polaron energy decreases when the temperature increases, and contrariwise for the attractive branch. Our results for the effective mass and the dynamical structure factor of the polaron show unambiguously that its quasiparticle nature disappears close to the critical temperature, in agreement with recent experimental findings. Finally, we have also estimated the fraction of bosons in the condensate as well as the superfluid fraction, and we have concluded that the impurity hinders the condensation of the rest of bosons.

Adsorption of helium on a charged propeller molecule: hexaphenylbenzene

Physisorption on planar or curved graphitic surfaces or aromatic rings has been investigated by various research groups, but in these studies, the substrate was usually strictly rigid. Here, we report a combined experimental and theoretical study of helium adsorption on cationic hexaphenylbenzene (HPB), a propeller-shaped molecule. The orientation of its propeller blades is known to be sensitive to the environment, with substantial differences between the molecule in the gas phase and in the crystalline solid. Mass spectra of He_{n}HPB^{+}, synthesized in helium nanodroplets, indicate enhanced stability for ions containing n = 2, 4, 14, 28, 42, 44, or 46 helium atoms. Path-integral molecular dynamics simulations reveal a significant dependence of the dissociation energy on the details of the HPB geometry. Good agreement between the experimental data and calculated dissociation energies is obtained, provided that the symmetry of HPB^{+} is reduced from D_{6} to D_{2}, such a lower symmetry being suggested from quantum chemical calculations as arising upon electron removal.Graphic

Adsorption of Helium on Small Cationic PAHs: Influence of Hydrocarbon Structure on the Microsolvation Pattern.

The adsorption of up to ∼100 helium atoms on cations of the planar polycyclic aromatic hydrocarbons (PAHs) anthracene, phenanthrene, fluoranthene, and pyrene was studied by combining helium nanodroplet mass spectrometry with classical and quantum computational methods. Recorded time-of-flight mass spectra reveal a unique set of structural features in the ion abundance as a function of the number of attached helium atoms for each of the investigated PAHs. Path-integral molecular dynamics simulations were used with a polarizable potential to determine the underlying adsorption patterns of helium around the studied PAH cations and in good general agreement with the experimental data. The calculated structures of the helium–PAH complexes indicate that the arrangement of adsorbed helium atoms is highly sensitive toward the structure of the solvated PAH cation. Closures of the first solvation shell around the studied PAH cations are suggested to lie between 29 and 37 adsorbed helium atoms depending on the specific PAH cation. Helium atoms are found to preferentially adsorb on these PAHs following the commensurate pattern common for graphitic surfaces, in contrast to larger carbonaceous molecules like corannulene, coronene, and fullerenes that exhibit a 1 × 1 commensurate phase.

Quantum effects in the structural and elastic properties of graphite: Path-integral simulations

Graphite, as a well-known carbon-based solid, is a paradigmatic example of the so-called van der Waals layered materials, which display a large anisotropy in their physical properties. Here we study quantum effects in structural and elastic properties of graphite by using path-integral molecular dynamics simulations in the temperature range from 50 to 1500~K. This method takes into account quantization and anharmonicity of vibrational modes in the material. Our results are compared with those found by using classical molecular dynamics simulations. We analyze the volume and in-plane area as functions of temperature and external stress. The quantum motion is essential to correctly describe the in-plane and out-of-plane thermal expansion. Quantum effects cause also changes in the elastic properties of graphite with respect to a classical model. At low temperature we find an appreciable decrease in the linear elastic constants, mainly in $C_{12}$ and $C_{44}$. Quantum corrections in stiffness constants can be in some cases even larger than 20\%. The bulk modulus and Poisson's ratio are reduced in a 4\% and 19\%, respectively, due to zero-point motion of the C atoms. These quantum effects in structural and elastic properties of graphite are nonnegligible up to temperatures higher than 300~K.

Two distinct non-Arrhenius behaviors of hydrogen diffusivities in fcc aluminum, silver, and copper determined byab initiopath integral simulations

Nuclear quantum effects (NQEs) are highly important for understanding a host of kinetic processes that occur with the participation of H (e.g., H adsorption, diffusion, permeation, and trapping in materials). In this paper, ab initio path integral molecular dynamics simulations were used to investigate NQEs on the lattice diffusion of H in common face-centered cubic (fcc) metals such as Al, Ag, and Cu over a wide temperature range of 75--1200 K (75--900 K for Al). We determined that the dependence of H diffusivities on temperature in Ag and Cu has a ``reversed-S'' shape on Arrhenius plots, as confirmed for fcc Pd in our recent study [H. Kimizuka et al., Phys. Rev. B 100, 024104 (2019)]. This result illustrates that the phenomenon is common in many fcc metals in which H atoms prefer to occupy octahedral sites. On the other hand, in the case of Al, in which H atoms prefer to occupy tetrahedral sites, the dependence of H diffusivities on temperature exhibits a familiar ``C'' shape. Such counterintuitive behavior is ascribed to differences in the dependence on temperature of the activation barriers for H migration between both types of fcc metals; this is due to the NQEs involving a competition between deceleration of H migration, which becomes effective at high temperatures because of zero-point vibrations, and acceleration of H migration, which becomes effective at low temperatures because of quantum tunneling. The dominance of the two mechanisms is determined by the coupling of the NQEs and the site preference of H depending on the metal. This finding has important implications for the interpretation of kinetic processes involving the crossover from classical to quantum behavior of H atoms jumping between different types of interstitial sites.

Quantum polyamorphism in compressed distinguishable helium-4.

We demonstrate that two amorphous solid states can exist in 4He consisting of distinguishable Boltzmann atoms under compressed conditions. The isothermal compression of normal or supercritical fluid 4He was conducted at 3-25K using the isobaric-isothermal path integral centroid molecular dynamics simulation. The compression of fluid first produced the low-dispersion amorphous (LDA) state possessing modest extension of atomic necklaces. Further isothermal compression up to the order of 10 kbar to 1 Mbar or an isobaric cooling of LDA induced the transition to the high-dispersion amorphous (HDA) state. The HDA was characterized by long quantum wavelengths of atoms extended over several Angstroms and the promotion of atomic residual diffusion. They were related to the quantum tunneling of atoms bestriding the potential saddle points in this glass. The change in pressure or temperature induced the LDA-HDA transition reversibly with hysteresis, while it resembled the coil-globule transition of classical polymers. The HDA had lower kinetic and higher Gibbs free energies than the LDA at close temperature. The HDA was absent at T ≥ 13 K, while the LDA-HDA transition pressure significantly decreased with lowering temperature. The LDA and HDA correspond to the trapped and tunneling regimes proposed by Markland et al. [J. Chem. Phys. 136, 074511 (2012)], respectively. The same reentrant behavior as they found was observed for the expansion factor of the quantum wavelength as well as for atomic diffusivity.

Isotope Quantum Effects in the Metallization Transition in Liquid Hydrogen.

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

Isotopic effects in chair graphane

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

Nuclear quantum effects in thermal conductivity from centroid molecular dynamics.

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

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

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

Quantum and Classical Proton Diffusion in Superconducting Clathrate Hydrides.

The discovery of near room temperature superconductivity in clathrate hydrides has ignited the search for both higher temperature superconductors and deeper understanding of the underlying physical phenomena. In a conventional electron-phonon mediated picture for the superconductivity for these materials, the high critical temperatures predicted and observed can be ascribed to the low mass of the protons, but this also poses nontrivial questions associated with how the proton dynamics affect the superconductivity. Using clathrate superhydride Li_{2}MgH_{16} as an example, we show through abinitio path integral simulations that proton diffusion in this system is remarkably high, with a diffusion coefficient, for example, reaching 6×10^{-6} cm^{2}/s at 300K and 250GPa. The diffusion is achieved primarily through proton transfer among interstitial voids within the otherwise rigid Li_{2}Mg sublattice at these conditions. The findings indicate the coexistence of proton quantum diffusion together with hydrogen-induced superconductivity, with implications for other very-high-temperature superconducting hydrides.