Machine learning-accelerated path integral molecular dynamics simulations of reactive organic electrolytes.

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.

The impact of quantum nuclear effects on hydrogen (H-) bond strength has been inferred in earlier work from bond lengths obtained from path integral molecular dynamics (PIMD) simulations. To obtain a direct quantitative assessment of such effects, we use constrained-centroid PIMD simulations to calculate the free energy changes upon breaking the H-bonds in dimers of HF and water. Comparing ab initio simulations performed using PIMD and classical nucleus molecular dynamics (MD), we find smaller dissociation free energies with the PIMD method. Specifically, at 50 K, the H-bond in (HF)(2) is about 30% weaker when quantum nuclear effects are included, while that in (H(2)O)(2) is about 15% weaker. In a complementary set of simulations, we compare unconstrained PIMD and classical nucleus MD simulations to assess the influence of quantum nuclei on the structures of these systems. We find increased heavy atom distances, indicating weakening of the H-bond consistent with that observed by direct calculation of the free energies of dissociation.

        D/H ratios serve as a key toolkit for tracing water and organics in forms of gases, liquids and ices. Therefore, accurate estimations for H/D isotope effects during the formation and transfer for these H-incorporated compounds at super-cold conditions (around 10 to 200 K) with different phases (gas, liquid and solid) are urged. We have developed theoretical protocols for gaseous molecules in a previous work [1] by introducing the so-called full partition function ratios (FPFRs). However, such methods are rather complex and unpractical for condensed systems such like water ices, liquid methane and so on. In this study, we developed new protocols for accurate estimations of these equilibrium H/D isotope fractionations especially for condensed systems (e.g., water, ices) by combining path integral molecular dynamics (PIMD) simulations (modified for improving its convergence at low temperatures) with deep potentials (DP) obtained with machine learning. Preliminary results will be introduced and discussed.  

The dynamics of the K atom attached to the 300-atom helium cluster upon photoexcitation has been studied theoretically based on the quantum molecular dynamics method. With the use of a hybrid method between time-dependent quantum dynamics and semiclassical path integral centroid molecular dynamics, the calculation includes the electronic state mixing for the K(2P) states as well as the quantum motions of helium atoms. It was found from the calculated trajectories that the K atom mostly desorbs by itself from the surface of the helium cluster without depending on the types of initial electronic excited states. It was also found that the K*He exciplex can be formed only as a minor product, following the initial excitation to the Π-type electronic state alone. Helium solvation structures of the K*Hen exciplexes were also investigated by path integral molecular dynamics simulations. The result implies that stable exciplexes can be formed up to n = 6 on the lowest excited state until the first solvation shell around the 4p orbital of K* is completed. Meanwhile, it was found that only two helium atoms can be bound to the K(2P) atom on the second lowest excited state.

A path integral centroid molecular dynamics method for Bose and Fermi statistics

Equilibrium magnesium isotope fractionation between aqueous Mg2+ and carbonate minerals: Insights from path integral molecular dynamics

We consider a steady-state (but transient) situation in which a warm dense aggregate is a two-temperature system with equilibrium electrons at temperature Te, ions at Ti, and Te ≠ Ti. Such states are achievable by pump–probe experiments. For warm dense hydrogen in such a two-temperature situation, we investigate nuclear quantum effects (NQEs) on structure and thermodynamic properties, thereby delineating the limitations of ordinary ab initio molecular dynamics. We use path integral molecular dynamics (PIMD) simulations driven by orbital-free density functional theory (OFDFT) calculations with state-of-the-art noninteracting free-energy and exchange-correlation functionals for the explicit temperature dependence. We calibrate the OFDFT calculations against conventional (explicit orbitals) Kohn–Sham DFT. We find that when the ratio of the ionic thermal de Broglie wavelength to the mean interionic distance is larger than about 0.30, the ionic radial distribution function is meaningfully affected by the inclusion of NQEs. Moreover, NQEs induce a substantial increase in both the ionic and electronic pressures. This confirms the importance of NQEs for highly accurate equation-of-state data on highly driven hydrogen. For Te > 20 kK, increasing Te in the warm dense hydrogen has slight effects on the ionic radial distribution function and equation of state in the range of densities considered. In addition, we confirm that compared with thermostatted ring-polymer molecular dynamics, the primitive PIMD algorithm overestimates electronic pressures, a consequence of the overly localized ionic description from the primitive scheme.

Thermal properties of graphene monolayers are studied by path-integral molecular dynamics simulations, which take into account the quantization of vibrational modes in the crystalline membrane and allow one to consider anharmonic effects in these properties. This system was studied at temperatures in the range from 12 to 2000 K and zero external stress, by describing the interatomic interactions through the LCBOPII effective potential. We analyze the internal energy and specific heat and compare the results derived from the simulations with those yielded by a harmonic approximation for the vibrational modes. This approximation turns out to be rather precise up to temperatures of about 400 K. At higher temperatures, we observe an influence of the elastic energy due to the thermal expansion of the graphene sheet. Zero-point and thermal effects on the in-plane and "real" surface of graphene are discussed. The thermal expansion coefficient α of the real area is found to be positive at all temperatures, in contrast to the expansion coefficient αp of the in-plane area, which is negative at low temperatures and becomes positive for T ≳ 1000 K.

Isolated hydrogen, deuterium, and muonium in diamond have been studied by path-integral molecular dynamics simulations in the canonical ensemble. Finite-temperature properties of these point defects were analyzed in the range from 100 to $800\phantom{\rule{0.3em}{0ex}}\mathrm{K}$. Interatomic interactions were modeled by a tight-binding potential fitted to density-functional calculations. The most stable position for these hydrogenic impurities is found at the $\mathrm{C}\text{\ensuremath{-}}\mathrm{C}$ bond center. Vibrational frequencies have been obtained from a linear-response approach, based on correlations of atom displacements at finite temperatures. The results show a large anharmonic effect in impurity vibrations at the bond center site, which hardens the vibrational modes with respect to a harmonic approximation. Zero-point motion causes an appreciable shift of the defect level in the electronic gap, as a consequence of electron-phonon interaction. This defect level goes down by $70\phantom{\rule{0.3em}{0ex}}\mathrm{meV}$ when replacing hydrogen by muonium.

We describe how the fast quasi-centroid molecular dynamics (f-QCMD) method can be applied to condensed-phase systems by approximating the quasi-centroid potential of mean force as a sum of inter- and intramolecular corrections to the classical interaction potential. The corrections are found by using a regularized iterative Boltzmann inversion procedure to recover the inter- and intramolecular quasi-centroid distribution functions obtained from a path integral molecular dynamics simulation. The resulting methodology is found to give good agreement with a previously published QCMD dipole absorption spectrum for liquid water and satisfactory agreement for ice. It also gives good agreement with spectra from a recent implementation of CMD that uses a precomputed elevated temperature potential of mean force. Modern centroid molecular dynamics methods, therefore, appear to be reaching a consensus regarding the impact of nuclear quantum effects on the vibrational spectra of water and ice.

Starting with separated atoms on the surface of helium clusters 4HeN, and as a first step to assess the formation of rubidium dimers Rb2 in the triplet state, we perform Path Integral Molecular Dynamics simulations in the NVT canonical ensemble. Based on an accurate potential energy surface (PES) for the He–Rb2(3Σ u +) interaction [Guillon et al., J. Chem. Phys. 136, 174307 (2012)], the total PES is analytically described as the addition of pair interactions. The i-PI code [Ceriotti et al., Comput. Phys. Commun. 185, 1019 (2014)] was used to perform the simulations. At a temperature of 2 K, clusters containing up to N = 70 helium atoms, with a number up to 200 beads per particle to describe quantum effects, were considered.

Nuclear quantum effect on intramolecular hydrogen bond of hydrogen maleate anion: An ab initio path integral molecular dynamics study

The monoprotonated compound N,N',N''-tris(p-tolyl)azacalix[3](2,6)pyridine (TAPH) contains an intramolecular hydrogen bond that is formed from three N atoms in its cavity. Constrained by the macrocyclic molecular structure, the separations between the N atoms in this bifurcated hydrogen bond are about 2.6 Å, considerably shorter than those typically observed for hydrogen bonded systems in the condensed phases. As such, TAPH exhibits significantly elongated N-H lengths in its hydrogen bond and a downfield 1H NMR chemical shift of 22.1 ppm. In this work, we carry out ab initio molecular dynamics and ab initio path integral molecular dynamics simulations of TAPH in the acetonitrile solution to reveal the geometry and proton sharing conditions of the bifurcated short hydrogen bond and uncover how the interplay of electronic and nuclear quantum effects gives rise to its far downfield 1H chemical shift. Taking a linear short hydrogen bond as a reference, we demonstrate the distinct features of competing quantum effects and electronic shielding effects in the bifurcated hydrogen bond of TAPH. We further use the degree of deshielding on the proton as a measure of the hydrogen bonding interactions and evaluate the strength of the bifurcated short hydrogen bond as compared to its linear counterpart.

Path integral molecular dynamics (PIMD) simulations of the extended simple point charge (SPC/E) model of liquid water (both H2O and D2O) have been performed in order to find the temperature at which the maximum liquid density is obtained. For comparison, purely classical SPC and SPC/E water have also been simulated over long periods (more than 200 ps). Structural properties and various temperature dependent quantities are reported. Special attention has been given to the fluctuations of pressure and volume in simulations run in the canonical (constant NVT) and isothermal–isobaric (constant NPT) ensembles, respectively. Although a density maximum is detected by monitoring energy–volume correlations as a function of temperature, the large statistical uncertainty in the correlations reduces the significance of this finding. The main conclusions of this work are threefold. First, the use of the energy–volume correlation moment to determine the density maximum appears to be a useful method, although a number of very long simulations of between 0.5–1.0 ns over a fairly wide temperature range are necessary to obtain accurate results. Second, neither the SPC nor the SPC/E water models are able to reproduce the density of cold water accurately. Finally, a quantum mechanical treatment of the water model is necessary to correctly represent the structure, energy, and fluctuations of the condensed phase.

Path integral molecular dynamics for hydrogen adsorption site of zeolite-templated carbon with semi-empirical PM3 potential