How to remove the spurious resonances from ring polymer molecular dynamics

Hydrogen-bonding and nuclear quantum effects in clays

Hydrogen bonds are of paramount importance in the chemistry of clays, mediating the interaction between the clay surface and water, and for some materials between separate layers. It is well-established that the accuracy of a computational model for clays depends on the level of theory at which the electronic structure is treated. However, for hydrogen-bonded systems, the motion of light H nuclei on the electronic potential energy surface is often affected by quantum delocalization. Using path integral molecular dynamics, we show that nuclear quantum effects lead to a relatively small change in the structure of clays, but one that is comparable to the variation incurred by treating the clay at different levels of electronic structure theory. Accounting for quantum effects weakens the hydrogen bonds in clays, with H-bonds between different layers of the clay affected more than those within the same layer; this is ascribed to the fact that the confinement of an H atom inside a layer is independent of its participation in hydrogen-bonding. More importantly, the weakening of hydrogen bonds by nuclear quantum effects causes changes in the vibrational spectra of these systems, significantly shifting the O–H stretching peaks and meaning that in order to fully understand these spectra by computational modeling, both electronic and nuclear quantum effects must be included. We show that after reparameterization of the popular clay forcefield CLAYFF, the O–H stretching region of their vibrational spectra better matches the experimental one, with no detriment to the model’s agreement with other experimental properties.

Molecular Dynamics with Constrained Nuclear Electronic Orbital Density Functional Theory: Accurate Vibrational Spectra from Efficient Incorporation of Nuclear Quantum Effects

Nuclear quantum effects play a crucial role in many chemical and biological systems involving hydrogen atoms yet are difficult to include in practical molecular simulations. In this paper, we combine our recently developed methods of constrained nuclear-electronic orbital density functional theory (cNEO-DFT) and constrained minimized energy surface molecular dynamics (CMES-MD) to create a new method for accurately and efficiently describing nuclear quantum effects in molecular simulations. By use of this new method, dubbed cNEO-MD, the vibrational spectra of a set of small molecules are calculated and compared with those from conventional ab initio molecular dynamics (AIMD) as well as from experiments. With the same formal scaling, cNEO-MD greatly outperforms AIMD in describing the vibrational modes with significant hydrogen motion characters, demonstrating the promise of cNEO-MD for simulating chemical and biological systems with significant nuclear quantum effects.

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.

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.

Resolving the Structural Debate for the Hydrated Excess Proton in Water

It has long been proposed that the hydrated excess proton in water (aka the solvated "hydronium" cation) likely has two limiting forms, that of the Eigen cation (H₉O₄⁺) and that of the Zundel cation (H₅O₂⁺). There has been debate over which of these two is the more dominant species and/or whether intermediate (or "distorted") structures between these two limits are the more realistic representation. Spectroscopy experiments have recently provided further results regarding the excess proton. These experiments show that the hydrated proton has an anisotropy reorientation time scale on the order of 1-2 ps. This time scale has been suggested to possibly contradict the picture of the more rapid "special pair dance" phenomenon for the hydrated excess proton, which is a signature of a distorted Eigen cation. The special pair dance was predicted from prior computational studies in which the hydrated central core hydronium structure continually switches (O-H···O)* special pair hydrogen-bond partners with the closest three water molecules, yielding on average a distorted Eigen cation with three equivalent and dynamically exchanging distortions. Through state-of-art simulations it is shown here that anisotropy reorientation time scales of the same magnitude are obtained that also include structural reorientations associated with the special pair dance, leading to a reinterpretation of the experimental results. These results and additional analyses point to a distorted and dynamic Eigen cation as the most prevalent hydrated proton species in aqueous acid solutions of dilute to moderate concentration, as opposed to a stabilized or a distorted (but not "dancing") Zundel cation.

Testing the quasicentroid molecular dynamics method on gas-phase ammonia

Quasicentroid molecular dynamics (QCMD) is a path-integral method for approximating nuclear quantum effects in dynamics simulations, which has given promising results for gas- and condensed-phase water. In this work, by simulating the infrared spectrum of gas-phase ammonia, we test the feasibility of extending QCMD beyond water. Overall, QCMD works as well for ammonia as for water, reducing or eliminating blue shifts from the classical spectrum without introducing the artificial red shifts or broadening associated with other imaginary-time path-integral methods. However, QCMD gives only a modest improvement over the classical spectrum for the position of the symmetric bend mode, which is highly anharmonic (since it correlates with the inversion pathway). We expect QCMD to have similar problems with large-amplitude degrees of freedom in other molecules but otherwise to work as well as for water.

Helium Nanodroplet Infrared Action Spectroscopy of the Proton-Bound Dimer of Hydrogen Sulfate and Formate: Examining Nuclear Quantum Effects

The proton-bound dimer of hydrogen sulfate and formate is an archetypal structure for ionic hydrogen-bonding complexes that contribute to biogenic aerosol nucleation. Of central importance for the structure and properties of this complex is the location of the bridging proton connecting the two conjugate base moieties. The potential energy surface for bridging proton translocation features two local minima, with the proton localized at either the formate or hydrogen sulfate moiety. However, electronic structure methods reveal a shallow potential energy surface governing proton translocation, with a barrier on the order of the zero-point energy. This shallow potential complicates structural assignment and necessitates a consideration of nuclear quantum effects. In this work, we probe the structure of this complex and its isotopologues, utilizing infrared (IR) action spectroscopy of ions captured in helium nanodroplets. The IR spectra indicate a structure in which a proton is shared between the hydrogen sulfate and formate moieties, HSO₄^-···H⁺···^-OOCH. However, because of the nuclear quantum effects and vibrational anharmonicities associated with the shallow potential for proton translocation, the extent of proton displacement from the formate moiety remains unclear, requiring further experiments or more advanced theoretical treatments for additional insight.

Efficient Quantum Vibrational Spectroscopy of Water with High-Order Path Integrals: From Bulk to Interfaces

Vibrational spectroscopy is key in probing the interplay between the structure and dynamics of aqueous systems. To map different regions of experimental spectra to the microscopic structure of a system, it is important to combine them with first-principles atomistic simulations that incorporate the quantum nature of nuclei. Here we show that the large cost of calculating the quantum vibrational spectra of aqueous systems can be dramatically reduced compared with standard path integral methods by using approximate quantum dynamics based on high-order path integrals. Together with state-of-the-art machine-learned electronic properties, our approach gives an excellent description not only of the infrared and Raman spectra of bulk water but also of the 2D correlation and the more challenging sum-frequency generation spectra of the water-air interface. This paves the way for understanding complex interfaces such as water encapsulated between or in contact with hydrophobic and hydrophilic materials through robust and inexpensive surface-sensitive and multidimensional spectra with first-principles accuracy.

On the "Matsubara heating" of overtone intensities and Fermi splittings

Classical molecular dynamics (MD) and imaginary-time path-integral dynamics methods underestimate the infrared absorption intensities of overtone and combination bands by typically an order of magnitude. Plé et al. [J. Chem. Phys. 155, 2863 (2021)] have shown that this is because such methods fail to describe the coupling of the centroid to the Matsubara dynamics of the fluctuation modes; classical first-order perturbation theory (PT) applied to the Matsubara dynamics is sufficient to recover most of the lost intensity in simple models and gives identical results to quantum (Rayleigh-Schrödinger) PT. Here, we show numerically that the results of this analysis can be used as post-processing correction factors, which can be applied to realistic (classical MD or path-integral dynamics) simulations of infrared spectra. We find that the correction factors recover most of the lost intensity in the overtone and combination bands of gas-phase water and ammonia and much of it for liquid water. We then re-derive and confirm the earlier PT analysis by applying canonical PT to Matsubara dynamics, which has the advantage of avoiding secular terms and gives a simple picture of the perturbed Matsubara dynamics in terms of action-angle variables. Collectively, these variables "Matsubara heat" the amplitudes of the overtone and combination vibrations of the centroid to what they would be in a classical system with the oscillators (of frequency Ω_i) held at their quantum effective temperatures [of ℏΩ_i coth(βℏΩ_i/2)/2k_B]. Numerical calculations show that a similar neglect of "Matsubara heating" causes path-integral methods to underestimate Fermi resonance splittings.

Anharmonic spectral features via trajectory-based quantum dynamics: A perturbative analysis of the interplay between dynamics and sampling

The performance of different approximate algorithms for computing anharmonic features in vibrational spectra is analyzed and compared on model and more realistic systems that present relevant nuclear quantum effects. The methods considered combine approximate sampling of the quantum thermal distribution with classical time propagation and include Matsubara dynamics, path integral dynamics approaches, linearized initial value representation, and the recently introduced adaptive quantum thermal bath. A perturbative analysis of these different methods enables us to account for the observed numerical performance on prototypes for overtones and combination bands and to draw qualitatively correct trends for the numerical results obtained for Fermi resonances. Our results prove that the unequal performances of these approaches often derive from the method employed to sample initial conditions and not, as usually assumed, from the lack of coherence in the time propagation. Furthermore, as confirmed by the analysis reported in Benson and Althorpe, J. Chem. Phys. 130, 194510 (2021), we demonstrate, both via the perturbative approach and numerically, that path integral dynamics methods fail to reproduce the intensities of these anharmonic features and follow purely classical trends with respect to their temperature behavior. Finally, the remarkably accurate performance of the adaptive quantum thermal bath approach is documented and motivated.

Nuclear Quantum Effects in Liquid Water at Near Classical Computational Cost Using the Adaptive Quantum Thermal Bath

We demonstrate the accuracy and efficiency of a recently introduced approach to account for nuclear quantum effects (NQEs) in molecular simulations: the adaptive quantum thermal bath (adQTB). In this method, zero-point energy is introduced through a generalized Langevin thermostat designed to precisely enforce the quantum fluctuation-dissipation theorem. We propose a refined adQTB algorithm with improved accuracy and report adQTB simulations of liquid water. Through extensive comparison with reference path integral calculations, we demonstrate that it provides excellent accuracy for a broad range of structural and thermodynamic observables as well as infrared vibrational spectra. The adQTB has a computational cost comparable to that of classical molecular dynamics, enabling simulations of up to millions of degrees of freedom.

Comparative studies of IR spectra of deprotonated serine with classical and thermostated ring polymer molecular dynamics simulations

Here we report the vibrational spectra of deprotonated serine calculated from the classical molecular dynamics (MD) simulations and thermostated ring-polymer molecular dynamics (TRPMD) simulation with third-order density-functional tight-binding. In our earlier study [Inakollu and Yu, "A systematic benchmarking of computational vibrational spectroscopy with DFTB3: Normal mode analysis and fast Fourier transform dipole autocorrelation function," J. Comput. Chem. 39, 2067 (2018)] of deprotonated serine, we observed a significant difference in the vibrational spectra with the classical MD simulations compared to the infrared multiple photon dissociation spectra. It was postulated that this is due to neglecting the nuclear quantum effects (NQEs). In this work, NQEs are considered in spectral calculation using the TRPMD simulations. With the help of potential of mean force calculations, the conformational space of deprotonated serine is analyzed and used to understand the difference in the spectra of classical MD and TRPMD simulations at 298.15 and 100 K. The high-frequency vibrational bands in the spectra are characterized using Fourier transform localized vibrational mode (FT-ν_N AC) and interatomic distance histograms. At room temperature, the quantum effects are less significant, and the free energy profiles in the classical MD and the TRPMD simulations are very similar. However, the hydrogen bond between the hydroxyl-carboxyl bond is slightly stronger in TRPMD simulations. At 100 K, the quantum effects are more prominent, especially in the 2600-3600 cm^-1, and the free energy profile slightly differs between the classical MD and TRPMD simulations. Using the FT-ν_N AC and the interatomic distance histograms, the high-frequency vibrational bands are discussed in detail.

Nuclear Quantum Effect and Its Temperature Dependence in Liquid Water from Random Phase Approximation via Artificial Neural Network

We report structural and dynamical properties of liquid water described by the random phase approximation (RPA) correlation together with the exact exchange energy (EXX) within density functional theory. By utilizing thermostated ring polymer molecular dynamics, we examine the nuclear quantum effects and their temperature dependence. We circumvent the computational limitation of performing direct first-principles molecular dynamics simulation at this high level of electronic structure theory by adapting an artificial neural network model. We show that the EXX+RPA level of theory accurately describes liquid water in terms of both dynamical and structural properties.

Machine learning phase space quantum dynamics approaches

Derived from phase space expressions of the quantum Liouville theorem, equilibrium continuity dynamics is a category of trajectory-based phase space dynamics methods, which satisfies the two critical fundamental criteria: conservation of the quantum Boltzmann distribution for the thermal equilibrium system and being exact for any thermal correlation functions (even of nonlinear operators) in the classical and harmonic limits. The effective force and effective mass matrix are important elements in the equations of motion of equilibrium continuity dynamics, where only the zeroth term of an exact series expansion of the phase space propagator is involved. We introduce a machine learning approach for fitting these elements in quantum phase space, leading to a much more efficient integration of the equations of motion. Proof-of-concept applications to realistic molecules demonstrate that machine learning phase space dynamics approaches are possible as well as competent in producing reasonably accurate results with a modest computation effort.

Progress and challenges in ab initio simulations of quantum nuclei in weakly bonded systems

Atomistic simulations based on the first-principles of quantum mechanics are reaching unprecedented length scales. This progress is due to the growth in computational power allied with the development of new methodologies that allow the treatment of electrons and nuclei as quantum particles. In the realm of materials science, where the quest for desirable emergent properties relies increasingly on soft weakly bonded materials, such methods have become indispensable. In this Perspective, an overview of simulation methods that are applicable for large system sizes and that can capture the quantum nature of electrons and nuclei in the adiabatic approximation is given. In addition, the remaining challenges are discussed, especially regarding the inclusion of nuclear quantum effects (NQEs) beyond a harmonic or perturbative treatment, the impact of NQEs on electronic properties of weakly bonded systems, and how different first-principles potential energy surfaces can change the impact of NQEs on the atomic structure and dynamics of weakly bonded systems.

Dynamic Nature of High-Pressure Ice VII

Starting from Shannon's definition of dynamic entropy, we propose a theory to describe the rare-event-determined dynamic states in condensed matter and their transitions and apply it to high-pressure ice VII. A dynamic intensive quantity named dynamic field, rather than the conventional thermodynamic intensive quantities such as temperature and pressure, is taken as the controlling variable. The dynamic entropy versus dynamic field curve demonstrates two dynamic states in the stability region of ice VII and dynamic ice VII. Their microscopic differences were assigned to the dynamic patterns of proton transfer. This study puts a similar dynamical theory used in earlier studies of glass models on a simpler and more fundamental basis, which could be applied to describe the dynamic states of more realistic condensed matter systems.

Dynamics & Spectroscopy with Neutrons-Recent Developments & Emerging Opportunities

This work provides an up-to-date overview of recent developments in neutron spectroscopic techniques and associated computational tools to interrogate the structural properties and dynamical behavior of complex and disordered materials, with a focus on those of a soft and polymeric nature. These have and continue to pave the way for new scientific opportunities simply thought unthinkable not so long ago, and have particularly benefited from advances in high-resolution, broadband techniques spanning energy transfers from the meV to the eV. Topical areas include the identification and robust assignment of low-energy modes underpinning functionality in soft solids and supramolecular frameworks, or the quantification in the laboratory of hitherto unexplored nuclear quantum effects dictating thermodynamic properties. In addition to novel classes of materials, we also discuss recent discoveries around water and its phase diagram, which continue to surprise us. All throughout, emphasis is placed on linking these ongoing and exciting experimental and computational developments to specific scientific questions in the context of the discovery of new materials for sustainable technologies.

Non-adiabatic Matsubara dynamics and non-adiabatic ring-polymer molecular dynamics

We present the non-adiabatic Matsubara dynamics, a general framework for computing the time-correlation function (TCF) of electronically non-adiabatic systems. This new formalism is derived based on the generalized Kubo-transformed TCF using the Wigner representation for both the nuclear degrees of freedom and the electronic mapping variables. By dropping the non-Matsubara nuclear normal modes in the quantum Liouvillian and explicitly integrating these modes out from the expression of the TCF, we derived the non-adiabatic Matsubara dynamics approach. Further making the approximation to drop the imaginary part of the Matsubara Liouvillian and enforce the nuclear momentum integral to be real, we arrived at the non-adiabatic ring-polymer molecular dynamics (NRPMD) approach. We have further justified the capability of NRPMD for simulating the non-equilibrium TCF. This work provides the rigorous theoretical foundation for several recently proposed state-dependent RPMD approaches and offers a general framework for developing new non-adiabatic quantum dynamics methods in the future.

Persistent Homology Metrics Reveal Quantum Fluctuations and Reactive Atoms in Path Integral Dynamics

Nuclear quantum effects (NQEs) are known to impact a number of features associated with chemical reactivity and physicochemical properties, particularly for light atoms and at low temperatures. In the imaginary time path integral formalism, each atom is mapped onto a "ring polymer" whose spread is related to the quantum mechanical uncertainty in the particle's position, i.e., its thermal wavelength. A number of metrics have previously been used to investigate and characterize this spread and explain effects arising from quantum delocalization, zero-point energy, and tunneling. Many of these shape metrics consider just the instantaneous structure of the ring polymers. However, given the significant interest in methods such as centroid molecular dynamics and ring polymer molecular dynamics that link the molecular dynamics of these ring polymers to real time properties, there exists significant opportunity to exploit metrics that also allow for the study of the fluctuations of the atom delocalization in time. Here we consider the ring polymer delocalization from the perspective of computational topology, specifically persistent homology, which describes the 3-dimensional arrangement of point cloud data, (i.e. atomic positions). We employ the Betti sequence probability distribution to define the ensemble of shapes adopted by the ring polymer. The Wasserstein distances of Betti sequences adjacent in time are used to characterize fluctuations in shape, where the Fourier transform and associated principal components provides added information differentiating atoms with different NQEs based on their dynamic properties. We demonstrate this methodology on two representative systems, a glassy system consisting of two atom types with dramatically different de Broglie thermal wavelengths, and ab initio molecular dynamics simulation of an aqueous 4 M HCl solution where the H-atoms are differentiated based on their participation in proton transfer reactions.

Nuclear quantum effects on the thermodynamic, structural, and dynamical properties of water

We perform path-integral molecular dynamics (PIMD) simulations of H₂O and D₂O using the q-TIP4P/F model. Simulations are performed at P = 1 bar and over a wide range of temperatures that include the equilibrium (T≥ 273 K) and supercooled (210 ≤T < 273 K) liquid states of water. The densities of both H₂O and D₂O calculated from PIMD simulations are in excellent agreement with experiments in the equilibrium and supercooled regimes. We also evaluate important thermodynamic response functions, specifically, the thermal expansion coefficient α_P(T), isothermal compressibility κ_T(T), isobaric heat capacity C_P(T), and static dielectric constant ε(T). While these properties are in excellent [α_P(T) and κ_T(T)] or semi-quantitative agreement [C_P(T) and ε(T)] with experiments in the equilibrium regime, they are increasingly underestimated upon further cooling. It follows that the inclusion of nuclear quantum effects in PIMD simulations of (q-TIP4P/F) water is not sufficient to reproduce the anomalous large fluctuations in density, entropy, and electric dipole moment characteristic of supercooled water. It has been hypothesized that water may exhibit a liquid-liquid critical point (LLCP) in the supercooled regime at P > 1 bar and that such a LLCP generates a maximum in C_P(T) and κ_T(T) at 1 bar. Consistent with this hypothesis and in particular, with experiments, we find a maximum in the κ_T(T) of q-TIP4P/F light and heavy water at T≈ 230-235 K. No maximum in C_P(T) could be detected down to T≥ 210 K. We also calculate the diffusion coefficient D(T) of H₂O and D₂O using the ring-polymer molecular dynamics (RPMD) technique and find that computer simulations are in remarkable good agreement with experiments at all temperatures studied. The results from RPMD/PIMD simulations are also compared with the corresponding results obtained from classical MD simulations of q-TIP4P/F water where atoms are represented by single interacting sites. Surprisingly, we find minor differences in most of the properties studied, with C_P(T), D(T), and structural properties being the only (expected) exceptions.

Neural network potential energy surface for the low temperature ring polymer molecular dynamics of the H2CO + OH reaction

A new potential energy surface (PES) and dynamical study of the reactive process of H₂CO + OH toward the formation of HCO + H₂O and HCOOH + H are presented. In this work, a source of spurious long range interactions in symmetry adapted neural network (NN) schemes is identified, which prevents their direct application for low temperature dynamical studies. For this reason, a partition of the PES into a diabatic matrix plus a NN many-body term has been used, fitted with a novel artificial neural network scheme that prevents spurious asymptotic interactions. Quasi-classical trajectory (QCT) and ring polymer molecular dynamics (RPMD) studies have been carried on this PES to evaluate the rate constant temperature dependence for the different reactive processes, showing good agreement with the available experimental data. Of special interest is the analysis of the previously identified trapping mechanism in the RPMD study, which can be attributed to spurious resonances associated with excitations of the normal modes of the ring polymer.

Simulating the ghost: quantum dynamics of the solvated electron

The nature of the bulk hydrated electron has been a challenge for both experiment and theory due to its short lifetime and high reactivity, and the need for a high-level of electronic structure theory to achieve predictive accuracy. The lack of a classical atomistic structural formula makes it exceedingly difficult to model the solvated electron using conventional empirical force fields, which describe the system in terms of interactions between point particles associated with atomic nuclei. Here we overcome this problem using a machine-learning model, that is sufficiently flexible to describe the effect of the excess electron on the structure of the surrounding water, without including the electron in the model explicitly. The resulting potential is not only able to reproduce the stable cavity structure but also recovers the correct localization dynamics that follow the injection of an electron in neat water. The machine learning model achieves the accuracy of the state-of-the-art correlated wave function method it is trained on. It is sufficiently inexpensive to afford a full quantum statistical and dynamical description and allows us to achieve accurate determination of the structure, diffusion mechanisms, and vibrational spectroscopy of the solvated electron.

A generalized class of strongly stable and dimension-free T-RPMD integrators

Recent work shows that strong stability and dimensionality freedom are essential for robust numerical integration of thermostatted ring-polymer molecular dynamics (T-RPMD) and path-integral molecular dynamics, without which standard integrators exhibit non-ergodicity and other pathologies [R. Korol et al., J. Chem. Phys. 151, 124103 (2019) and R. Korol et al., J. Chem. Phys. 152, 104102 (2020)]. In particular, the BCOCB scheme, obtained via Cayley modification of the standard BAOAB scheme, features a simple reparametrization of the free ring-polymer sub-step that confers strong stability and dimensionality freedom and has been shown to yield excellent numerical accuracy in condensed-phase systems with large time steps. Here, we introduce a broader class of T-RPMD numerical integrators that exhibit strong stability and dimensionality freedom, irrespective of the Ornstein-Uhlenbeck friction schedule. In addition to considering equilibrium accuracy and time step stability as in previous work, we evaluate the integrators on the basis of their rates of convergence to equilibrium and their efficiency at evaluating equilibrium expectation values. Within the generalized class, we find BCOCB to be superior with respect to accuracy and efficiency for various configuration-dependent observables, although other integrators within the generalized class perform better for velocity-dependent quantities. Extensive numerical evidence indicates that the stated performance guarantees hold for the strongly anharmonic case of liquid water. Both analytical and numerical results indicate that BCOCB excels over other known integrators in terms of accuracy, efficiency, and stability with respect to time step for practical applications.

Effective mass path integral simulations of quasiparticles in condensed phases

The quantum many-body problem in condensed phases is often simplified using a quasiparticle description, such as effective mass theory for electron motion in a periodic solid. These approaches are often the basis for understanding many fundamental condensed phase processes, including the molecular mechanisms underlying solar energy harvesting and photocatalysis. Despite the importance of these effective particles, there is still a need for computational methods that can explore their behavior on chemically relevant length and time scales. This is especially true when the interactions between the particles and their environment are important. We introduce an approach for studying quasiparticles in condensed phases by combining effective mass theory with the path integral treatment of quantum particles. This framework incorporates the generally anisotropic electronic band structure of materials into path integral simulation schemes to enable modeling of quasiparticles in quantum confinement, for example. We demonstrate the utility of effective mass path integral simulations by modeling an exciton in solid potassium chloride and electron trapping by a sulfur vacancy in monolayer molybdenum disulfide.

Ring-polymer, centroid, and mean-field approximations to multi-time Matsubara dynamics

Based on a recently developed generalization of Matsubara dynamics to the multi-time realm, we present a formal derivation of multi-time generalizations of ring-polymer molecular dynamics, thermostatted ring-polymer molecular dynamics (TRPMD), centroid molecular dynamics (CMD), and mean-field Matsubara dynamics. Additionally, we analyze the short-time accuracy of each methodology. We find that for multi-time correlation functions of linear operators, (T)RPMD is accurate up to order t³, while CMD is only correct up to t, indicating a degradation in the accuracy of these methodologies with respect to the single-time counterparts. The present work provides a firm justification for the use of path-integral-based approximations for the calculation of multi-time correlation functions.

Two-dimensional Raman spectroscopy of Lennard-Jones liquids via ring-polymer molecular dynamics

The simulation of multidimensional vibrational spectroscopy of condensed-phase systems including nuclear quantum effects is challenging since full quantum-mechanical calculations are still intractable for large systems comprising many degrees of freedom. Here, we apply the recently developed double Kubo transform (DKT) methodology in combination with ring-polymer molecular dynamics (RPMD) for evaluating multi-time correlation functions [K. A. Jung et al., J. Chem. Phys. 148, 244105 (2018)], providing a practical method for incorporating nuclear quantum effects in nonlinear spectroscopy of condensed-phase systems. We showcase the DKT approach in the simulation of the fifth-order two-dimensional (2D) Raman spectroscopy of Lennard-Jones liquids as a prototypical example, which involves nontrivial nonlinear spectroscopic observables of systems described by anharmonic potentials. Our results show that the DKT can faithfully reproduce the 2D Raman response of liquid xenon at high temperatures, where the system behaves classically. In contrast, liquid neon at low temperatures exhibits moderate but discernible nuclear quantum effects in the 2D Raman response compared to the responses obtained with classical molecular dynamics approaches. Thus, the DKT formalism in combination with RPMD simulations enables simulations of multidimensional optical spectroscopy of condensed-phase systems that partially account for nuclear quantum effects.

Inexpensive modeling of quantum dynamics using path integral generalized Langevin equation thermostats

The properties of molecules and materials containing light nuclei are affected by their quantum mechanical nature. Accurate modeling of these quantum nuclear effects requires computationally demanding path integral techniques. Considerable success has been achieved in reducing the cost of such simulations by using generalized Langevin dynamics to induce frequency-dependent fluctuations. Path integral generalized Langevin equation methods, however, have this far been limited to the study of static, thermodynamic properties due to the large perturbation to the system's dynamics induced by the aggressive thermostatting. Here, we introduce a post-processing scheme, based on analytical estimates of the dynamical perturbation induced by the generalized Langevin dynamics, which makes it possible to recover meaningful time correlation properties from a thermostatted trajectory. We show that this approach yields spectroscopic observables for model and realistic systems that have an accuracy comparable to much more demanding approximate quantum dynamics techniques based on full path integral simulations.

Dimension-free path-integral molecular dynamics without preconditioning

Convergence with respect to imaginary-time discretization (i.e., the number of ring-polymer beads) is an essential part of any path-integral-based molecular dynamics (MD) calculation. However, an unfortunate property of existing non-preconditioned numerical integration schemes for path-integral molecular dynamics-including essentially all existing ring-polymer molecular dynamics (RPMD) and thermostatted RPMD (T-RPMD) methods-is that for a given MD time step, the overlap between the exact ring-polymer Boltzmann-Gibbs distribution and that sampled using MD becomes zero in the infinite-bead limit. This has clear implications for hybrid Metropolis Monte Carlo/MD sampling schemes, and it also causes the divergence with bead number of the primitive path-integral kinetic-energy expectation value when using standard RPMD or T-RPMD. We show that these and other problems can be avoided through the introduction of "dimension-free" numerical integration schemes for which the sampled ring-polymer position distribution has non-zero overlap with the exact distribution in the infinite-bead limit for the case of a harmonic potential. Most notably, we introduce the BCOCB integration scheme, which achieves dimension freedom via a particular symmetric splitting of the integration time step and a novel implementation of the Cayley modification [R. Korol et al., J. Chem. Phys. 151, 124103 (2019)] for the free ring-polymer half-steps. More generally, we show that dimension freedom can be achieved via mollification of the forces from the external physical potential. The dimension-free path-integral numerical integration schemes introduced here yield finite error bounds for a given MD time step, even as the number of beads is taken to infinity; these conclusions are proven for the case of a harmonic potential and borne out numerically for anharmonic systems that include liquid water. The numerical results for BCOCB are particularly striking, allowing for nearly three-fold increases in the stable time step for liquid water with respect to the Bussi-Parrinello (OBABO) and Leimkuhler (BAOAB) integrators, while introducing negligible errors in the calculated statistical properties and absorption spectrum. Importantly, the dimension-free, non-preconditioned integration schemes introduced here preserve ergodicity and global second-order accuracy, and they remain simple, black-box methods that avoid additional computational costs, tunable parameters, or system-specific implementations.

Temperature dependence of the vibrational spectrum of porphycene: a qualitative failure of classical-nuclei molecular dynamics

Approximate quantum dynamics succeed in predicting a temperature-dependent blue-shift of the high-frequency stretch bands that arise from vibrational coupling between low-frequency thermally activated modes and high-frequency quantized ones. Classical nuclei molecular dynamics fail and instead predict a red-shift.

Which quantum statistics–classical dynamics method is best for water?

We make a systematic comparison of TRPMD, CMD, QCMD and LSC-IVR by calculating the infrared spectrum of water in its gas, liquid and ice phases.

Cayley modification for strongly stable path-integral and ring-polymer molecular dynamics

Path-integral-based molecular dynamics (MD) simulations are widely used for the calculation of numerically exact quantum Boltzmann properties and approximate dynamical quantities. A nearly universal feature of MD numerical integration schemes for equations of motion based on imaginary-time path integrals is the use of harmonic normal modes for the exact evolution of the free ring-polymer positions and momenta. In this work, we demonstrate that this standard practice creates numerical artifacts. In the context of conservative (i.e., microcanonical) equations of motion, it leads to numerical instability. In the context of thermostated (i.e., canonical) equations of motion, it leads to nonergodicity of the sampling. These pathologies are generally proven to arise at integration time steps that depend only on the system temperature and the number of ring-polymer beads, and they are numerically demonstrated for the cases of conventional ring-polymer MD (RPMD) and thermostated RPMD (TRPMD). Furthermore, it is demonstrated that these numerical artifacts are removed via replacement of the exact free ring-polymer evolution with a second-order approximation based on the Cayley transform. The Cayley modification introduced here can immediately be employed with almost every existing integration scheme for path-integral-based MD-including path-integral MD (PIMD), RPMD, TRPMD, and centroid MD-providing strong symplectic stability and ergodicity to the numerical integration, at no penalty in terms of computational cost, algorithmic complexity, or accuracy of the overall MD time step. Furthermore, it is shown that the improved numerical stability of the Cayley modification allows for the use of larger MD time steps. We suspect that the Cayley modification will therefore find useful application in many future path-integral-based MD simulations.

Reduced rovibrational coupling Cartesian dynamics for semiclassical calculations: Application to the spectrum of the Zundel cation

We study the vibrational spectrum of the protonated water dimer, by means of a divide-and-conquer semiclassical initial value representation of the quantum propagator, as a first step in the study of larger protonated water clusters. We use the potential energy surface from the work of Huang et al. [J. Chem. Phys. 122, 044308 (2005)]. To tackle such an anharmonic and floppy molecule, we employ fully Cartesian dynamics and carefully reduce the coupling to global rotations in the definition of normal modes. We apply the time-averaging filter and obtain clean power spectra relative to suitable reference states that highlight the spectral peaks corresponding to the fundamental excitations of the system. Our trajectory-based approach allows for the physical interpretation of the very challenging proton transfer modes. We find that it is important, for such a floppy molecule, to selectively avoid initially exciting lower energy modes, in order to obtain cleaner spectra. The estimated vibrational energies display a mean absolute error (MAE) of ∼29 cm^-1 with respect to available multiconfiguration time-dependent Hartree calculations and MAE ∼ 14 cm^-1 when compared to the optically active experimental excitations of the Ne-tagged Zundel cation. The reasonable scaling in the number of trajectories for Monte Carlo convergence is promising for applications to higher dimensional protonated cluster systems.

Sampling the thermal Wigner density via a generalized Langevin dynamics

The Wigner thermal density is a function of considerable interest in the area of approximate (linearized or semiclassical) quantum dynamics where it is employed to generate initial conditions for the propagation of appropriate sets of classical trajectories. In this paper, we propose an original approach to compute the Wigner density based on a generalized Langevin equation. The stochastic dynamics is nontrivial in that it contains a coordinate-dependent friction coefficient and a generalized force that couples momenta and coordinates. These quantities are, in general, not known analytically and have to be estimated via auxiliary calculations. The performance of the new sampling scheme is tested on standard model systems with highly nonclassical features such as relevant zero point energy effects, correlation between momenta and coordinates, and negative parts of the Wigner density. In its current brute force implementation, the algorithm, whose convergence can be systematically checked, is accurate and has only limited overhead compared to schemes with similar characteristics. We briefly discuss potential ways to further improve its numerical efficiency.

Multi-time formulation of Matsubara dynamics

Matsubara dynamics has recently emerged as the most general form of a quantum-Boltzmann-conserving classical dynamics theory for the calculation of single-time correlation functions. Here, we present a generalization of Matsubara dynamics for the evaluation of multitime correlation functions. We first show that the Matsubara approximation can also be used to approximate the two-time symmetrized double Kubo transformed correlation function. By a straightforward extension of these ideas to the multitime realm, a multitime Matsubara dynamics approximation can be obtained for the multitime fully symmetrized Kubo transformed correlation function. Although not a practical method, due to the presence of a phase-term, this multitime formulation of Matsubara dynamics represents a benchmark theory for future development of Boltzmann preserving semiclassical approximations to general higher order multitime correlation functions. It also reveals a connection between imaginary time-ordering in the path integral and the classical dynamics of multitime correlation functions.

Coarse-graining of many-body path integrals: Theory and numerical approximations

Feynman's imaginary time path integral approach to quantum statistical mechanics provides a theoretical formalism for including nuclear quantum effects (NQEs) in simulation of condensed matter systems. Sinitskiy and Voth [J. Chem. Phys. 143, 094104 (2015)] have presented the coarse-grained path integral (CG-PI) theory, which provides a reductionist coarse-grained representation of the imaginary time path integral based on the quantum-classical isomorphism. In this paper, the many-body generalization of the CG-PI theory is presented. It is shown that the N interacting particles obeying quantum Boltzmann statistics can be represented as a system of N pairs of classical-like pseudoparticles coupled to each other analogous to the pseudoparticle pair of the one-body theory. Moreover, we present a numerical CG-PI (n-CG-PI) method applying a simple approximation to the coupling scheme between the pseudoparticles due to numerical challenges of directly implementing the full many-body CG-PI theory. Structural correlations of two liquid systems are investigated to demonstrate the performance of the n-CG-PI method. Both the many-body CG-PI theory and the n-CG-PI method not only present reductionist views of the many-body quantum Boltzmann statistics but also provide theoretical and numerical insight into how to explicitly incorporate NQEs in the representation of condensed matter systems with minimal additional degrees of freedom.

Characterization of a trans-trans Carbonic Acid-Fluoride Complex by Infrared Action Spectroscopy in Helium Nanodroplets

The high Lewis basicity and small ionic radius of fluoride promote the formation of strong ionic hydrogen bonds in the complexation of fluoride with protic molecules. Herein, we report that carbonic acid, a thermodynamically disfavored species that is challenging to investigate experimentally, forms a complex with fluoride in the gas phase. Intriguingly, this complex is highly stable and is observed in abundance upon nanoelectrospray ionization of an aqueous sodium fluoride solution in the presence of gas-phase carbon dioxide. We characterize the structure and properties of the carbonic acid–fluoride complex, F^–(H₂CO₃), and its deuterated isotopologue, F^–(D₂CO₃), by helium nanodroplet infrared action spectroscopy in the photon energy range of 390–2800 cm^–1. The complex adopts a C_2v symmetry structure with the carbonic acid in a planar trans–trans conformation and both OH groups forming ionic hydrogen bonds with the fluoride. Substantial vibrational anharmonic effects are observed in the infrared spectra, most notably a strong blue shift of the symmetric hydrogen stretching fundamental relative to predictions from the harmonic approximation or vibrational second-order perturbation theory. Ab initio thermostated ring-polymer molecular dynamics simulations indicate that this blue shift originates from strong coupling between the hydrogen stretching and bending vibrations, resulting in an effective weakening of the OH···F^– ionic hydrogen bonds.

Elucidating the Nuclear Quantum Dynamics of Intramolecular Double Hydrogen Transfer in Porphycene

We address the double hydrogen transfer (DHT) dynamics of the porphycene molecule, a complex paradigmatic system in which the making and breaking of H-bonds in a highly anharmonic potential energy surface require a quantum mechanical treatment not only of the electrons but also of the nuclei. We combine density functional theory calculations, employing hybrid functionals and van der Waals corrections, with recently proposed and optimized path-integral ring-polymer methods for the approximation of quantum vibrational spectra and reaction rates. Our full-dimensional ring-polymer instanton simulations show that below 100 K the concerted DHT tunneling pathway dominates but between 100 and 300 K there is a competition between concerted and stepwise pathways when nuclear quantum effects are included. We obtain ground-state reaction rates of 2.19 × 10¹¹ s^–1 at 150 K and 0.63 × 10¹¹ s^–1 at 100 K, in good agreement with experiment. We also reproduce the puzzling N–H stretching band of porphycene with very good accuracy from thermostated ring-polymer molecular dynamics simulations. The position and line shape of this peak, centered at around 2600 cm^–1 and spanning 750 cm^–1, stem from a combination of very strong H-bonds, the coupling to low-frequency modes, and the access to cis-like isomeric conformations, which cannot be appropriately captured with classical-nuclei dynamics. These results verify the appropriateness of our general theoretical approach and provide a framework for a deeper physical understanding of hydrogen transfer dynamics in complex systems.

Quantum approaches to vibrational dynamics and spectroscopy: is ease of interpretation sacrificed as rigor increases?

The subject of this Perspective is quantum approaches, beyond the harmonic approximation, to vibrational dynamics and IR spectroscopy.