On computing spectral densities from classical, semiclassical, and quantum simulations.

The Caldeira-Leggett model provides a compact characterization of a thermal environment in terms of a spectral density function, which has led to a variety of numerically exact quantum methods for reduced density matrix propagation. Since spectral densities are often computed from classical molecular dynamics simulations, we investigate in this paper whether quantum effects should be accounted for in the calculations. Therefore, we reformulate the recently developed Fourier method for spectral density calculations from semiclassical simulations which approximately allow for quantum effects. We propose two possible protocols based on either correlation functions or expectation values. These protocols are tested on a generic Calderra-Leggett model for the linearized semiclassical initial-value representation (LSC-IVR), the thawed Gaussian wave packet dynamics (TGWD), and hybrid schemes combining the two with the more accurate Herman-Kluk formula. Surprisingly, spectral densities from the LSC-IVR method, which treats the dynamics completely classically, are found to be extremely accurate, even in the quantum regime, where this method does not give a correct description of the correlation functions and expectation values. In contrast, the TGWD method turns out as too inaccurate for spectral density calculations, and the hybrid schemes perform well only if the system is close to the classical regime. This implies that, if the bath has a Caldeira-Leggett form, spectral densities are insensitive to quantum effects and any effort to approximately account for them rather leads to errors. Hence, in this case, spectral densities can be computed from classical simulations and used in a reduced quantum simulation as well.

Simulating vibronic spectra via Matsubara-like dynamics: Coping with the sign problem.

Measuring vibronic spectra probes dynamical processes in molecular systems. When interpreted via suitable theoretical tools, the experimental data provides comprehensive information about the system in question. For complex many-body problems, such an approach usually requires the formulation of proper classical-like approximations, which is particularly challenging if multiple electronic states are involved. In this work, we express the imaginary-time shifted time correlation function and, thus, the vibronic spectrum in terms of the so-called Matsubara dynamics, which combines quantum statistics and classical-like dynamics. By applying the Matsubara approximation in the adiabatic limit, we derive a formal generalization of the existing Matsubara dynamics formalism to multiple potential energy surfaces (PESs), which, however, does not feature all the defining properties of its single-PES counterpart though suffering equally from the sign problem. The mathematical analysis for two shifted harmonic oscillators suggests a new modified method to practically simulate the standard correlation function via Matsubara-like dynamics. Importantly, this modified method samples the thermal Wigner function without suffering from the sign problem and yields an accurate approximation to the vibronic absorption spectrum, not only for the harmonic system but also for the anharmonic one.

Quasi-classical approaches to vibronic spectra revisited.

The framework to approach quasi-classical dynamics in the electronic ground state is well established and is based on the Kubo-transformed time correlation function (TCF), being the most classical-like quantum TCF. Here we discuss whether the choice of the Kubo-transformed TCF as a starting point for simulating vibronic spectra is as unambiguous as it is for vibrational ones. Employing imaginary-time path integral techniques in combination with the interaction representation allowed us to formulate a method for simulating vibronic spectra in the adiabatic regime that takes nuclear quantum effects and dynamics on multiple potential energy surfaces into account. Further, a generalized quantum TCF is proposed that contains many well-established TCFs, including the Kubo one, as particular cases. Importantly, it also provides a framework to construct new quantum TCFs. Applying the developed methodology to the generalized TCF leads to a plethora of simulation protocols, which are based on the well-known TCFs as well as on new ones. Their performance is investigated on 1D anharmonic model systems at finite temperatures. It is shown that the protocols based on the new TCFs may lead to superior results with respect to those based on the common ones. The strategies to find the optimal approach are discussed.

A time-correlation function approach to nuclear dynamical effects in X-ray spectroscopy.

Modern X-ray spectroscopy has proven itself as a robust tool for probing the electronic structure of atoms in complex environments. Despite working on energy scales that are much larger than those corresponding to nuclear motions, taking nuclear dynamics and the associated nuclear correlations into account may be of importance for X-ray spectroscopy. Recently, we have developed an efficient protocol to account for nuclear dynamics in X-ray absorption and resonant inelastic X-ray scattering spectra [Karsten et al., J. Phys. Chem. Lett. 8, 992 (2017)], based on ground state molecular dynamics accompanied with state-of-the-art calculations of electronic excitation energies and transition dipoles. Here, we present an alternative derivation of the formalism and elaborate on the developed simulation protocol using gas phase and bulk water as examples. The specific spectroscopic features stemming from the nuclear motions are analyzed and traced down to the dynamics of electronic energy gaps and transition dipole correlation functions. The observed tendencies are explained on the basis of a simple harmonic model, and the involved approximations are discussed. The method represents a step forward over the conventional approaches that treat the system in full complexity and provides a reasonable starting point for further improvements.

Nuclear Dynamical Correlation Effects in X-ray Spectroscopy from a Theoretical Time-Domain Perspective.

To date X-ray spectroscopy has become a routine tool that can reveal highly local and element-specific information on the electronic structure of atoms in complex environments. Here, we focus on nuclear dynamical correlation effects in X-ray spectra and develop a rigorous time-correlation function method employing ground state classical molecular dynamics simulations. The importance of nuclear correlation phenomena is demonstrated by comparison against the results from the conventional sampling approach performed on the same data set for gas phase water. In contrast to the first-order absorption, second-order resonant inelastic scattering spectra exhibit pronounced fingerprints of nuclear motions. The developed methodology is not biased to a particular electronic structure method and, owing to its generality, can be applied to, e.g., X-ray photoelectron and Auger spectroscopies.

Vibrational spectroscopy via the Caldeira-Leggett model with anharmonic system potentials.

The Caldeira-Leggett (CL) model, which describes a system bi-linearly coupled to a harmonic bath, has enjoyed popularity in condensed phase spectroscopy owing to its utmost simplicity. However, the applicability of the model to cases with anharmonic system potentials, as it is required for the description of realistic systems in solution, is questionable due to the presence of the invertibility problem [F. Gottwald et al., J. Phys. Chem. Lett. 6, 2722 (2015)] unless the system itself resembles the CL model form. This might well be the case at surfaces or in the solid regime, which we here confirm for a particular example of an iodine molecule in the atomic argon environment under high pressure. For this purpose we extend the recently proposed Fourier method for parameterizing linear generalized Langevin dynamics [F. Gottwald et al., J. Chem. Phys. 142, 244110 (2015)] to the non-linear case based on the CL model and perform an extensive error analysis. In order to judge on the applicability of this model in advance, we give practical empirical criteria and discuss the effect of the potential renormalization term. The obtained results provide evidence that the CL model can be used for describing a potentially broad class of systems.

Quantum free energy landscapes from ab initio path integral metadynamics: Double proton transfer in the formic acid dimer is concerted but not correlated.

With the goal of computing quantum free energy landscapes of reactive (bio)chemical systems in multi-dimensional space, we combine the metadynamics technique for sampling potential energy surfaces with the ab initio path integral approach to treating nuclear quantum motion. This unified method is applied to the double proton transfer process in the formic acid dimer (FAD), in order to study the nuclear quantum effects at finite temperatures without imposing a one-dimensional reaction coordinate or reducing the dimensionality. Importantly, the ab initio path integral metadynamics technique allows one to treat the hydrogen bonds and concomitant proton transfers in FAD strictly independently and thus provides direct access to the much discussed issue of whether the double proton transfer proceeds via a stepwise or concerted mechanism. The quantum free energy landscape we compute for this H-bonded molecular complex reveals that the two protons move in a concerted fashion from initial to product state, yet world-line analysis of the quantum correlations demonstrates that the protons are as quantum-uncorrelated at the transition state as they are when close to the equilibrium structure.

Applicability of the Caldeira-Leggett Model to Vibrational Spectroscopy in Solution.

Formulating a rigorous system-bath partitioning approach remains an open issue. In this context, the famous Caldeira-Leggett model that enables quantum and classical treatment of Brownian motion on equal footing has enjoyed popularity. Although this model is by any means a useful theoretical tool, its ability to describe anharmonic dynamics of real systems is often taken for granted. In this Letter, we show that the mapping between a molecular system under study and the model cannot be established in a self-consistent way, unless the system part of the potential is taken effectively harmonic. Mathematically, this implies that the mapping is not invertible. This "invertibility problem" is not dependent on the peculiarities of particular molecular systems and is rooted in the anharmonicity of the system part of the Caldeira-Leggett model potential.

Parametrizing linear generalized Langevin dynamics from explicit molecular dynamics simulations.

Fundamental understanding of complex dynamics in many-particle systems on the atomistic level is of utmost importance. Often the systems of interest are of macroscopic size but can be partitioned into a few important degrees of freedom which are treated most accurately and others which constitute a thermal bath. Particular attention in this respect attracts the linear generalized Langevin equation, which can be rigorously derived by means of a linear projection technique. Within this framework, a complicated interaction with the bath can be reduced to a single memory kernel. This memory kernel in turn is parametrized for a particular system studied, usually by means of time-domain methods based on explicit molecular dynamics data. Here, we discuss that this task is more naturally achieved in frequency domain and develop a Fourier-based parametrization method that outperforms its time-domain analogues. Very surprisingly, the widely used rigid bond method turns out to be inappropriate in general. Importantly, we show that the rigid bond approach leads to a systematic overestimation of relaxation times, unless the system under study consists of a harmonic bath bi-linearly coupled to the relevant degrees of freedom.

Theoretical spectroscopy using molecular dynamics: theory and application to CH5(+) and its isotopologues.

Infrared spectroscopy is a powerful technique to unravel the structure and dynamics of molecular systems of ever increasing complexity. For isolated molecules in the gas phase theoretical approaches that directly rely on solving the Schrödinger equation, either approximately or quasi-exactly, are well established. A distinctly different approach to compute infrared spectra can be based on advanced molecular dynamics, itself being based on classical Newtonian dynamics, in conjunction with concurrent first principles electronic structure calculations. At variance with traditional methods, which are formulated in terms of the Schrödinger representation of quantum mechanics, the molecular dynamics approach stems from Heisenberg's representation and thus relies on computing thermal expectation values of time-correlation functions. Crucial in addition to generating the spectra themselves is their decomposition in terms of modes, which can be assigned to correlated atomic motion. This ab initio molecular dynamics route to compute infrared spectra, and its recent extension to quasiclassical techniques relying on approximate path integral dynamics, is covered in the review part of this Perspective. The usefulness of this unconventional approach, which can be generalized beyond infrared spectroscopy, is demonstrated in detail by applying the full machinery in computing and assigning the infrared spectra of protonated methane and its isotopologues. This particular molecule is often considered to be the most prominent member of the class of floppy or fluxional molecules. CH5(+) has been a longstanding challenge for theoretical infrared spectroscopy because it undergoes intricate large-amplitude motion, which is also reviewed. Molecular dynamics based infrared spectroscopy is general and can be applied to diverse systems such as molecular complexes in the gas phase, chromophores in biomolecular environments, and solute-solvent systems in the liquid phase.

Microsolvation-induced quantum localization in protonated methane.

Nuclear quantum effects are responsible for vivid large-amplitude motion in protonated methane CH(5)(+), which enables so-called "hydrogen scrambling" that leads to the dynamical equivalence of all five protons even at low temperatures. But what is the impact of external perturbations on hydrogen scrambling of CH(5)(+) in this quantum fluxional ground state? We report ab initio path integral simulations of CH(5)(+)(H(2))(n), n = 1, 2, 3 that demonstrate cessation of hydrogen scrambling at low temperatures (20 K), but only slowdown at moderate temperatures (110 K). Importantly, different and unexpected mechanisms that are responsible for freezing the scrambling dynamics are revealed and traced back to distinct microsolvation patterns.

Quantum-induced symmetry breaking explains infrared spectra of CH(5)(+) isotopologues.

For decades, protonated methane, CH(5)(+), has provided new surprises and challenges for both experimentalists and theoreticians. This is because of the correlated large-amplitude motion of its five protons around the carbon nucleus, which leads to so-called hydrogen scrambling and causes a fluxional molecular structure. Here, the infrared spectra of all its H/D isotopologues have been measured using the 'Laser Induced Reactions' technique. Their shapes are found to be extremely dissimilar and depend strongly on the level of deuteration (only CD(5)(+) is similar to CH(5)(+)). All the spectra can be reproduced and assigned based on ab initio quantum simulations. The occupation of the topologically different sites by protons and deuterons is found to be strongly non-combinatorial and thus non-classical. This purely quantum-statistical effect implies a breaking of the classical symmetry of the site occupations induced by zero-point fluctuations, and this phenomenon is key to understanding the spectral changes studied here.

Communications: On artificial frequency shifts in infrared spectra obtained from centroid molecular dynamics: Quantum liquid water.

Centroid molecular dynamics (CMD) is a popular method to extract approximate quantum dynamics from path integral simulations. Very recently we have shown that CMD gas phase infrared spectra exhibit significant artificial redshifts of stretching peaks, due to the so-called "curvature problem" imprinted by the effective centroid potential. Here we provide evidence that for condensed phases, and in particular for liquid water, CMD produces pronounced artificial redshifts for high-frequency vibrations such as the OH stretching band. This peculiar behavior intrinsic to the CMD method explains part of the unexpectedly large quantum redshifts of the stretching band of liquid water compared to classical frequencies, which is improved after applying a simple and rough "harmonic curvature correction."

On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy.

Centroid molecular dynamics (CMD) and ring polymer molecular dynamics (RPMD) are two conceptually distinct extensions of path integral molecular dynamics that are able to generate approximate quantum dynamics of complex molecular systems. Both methods can be used to compute quasiclassical time correlation functions which have direct application in molecular spectroscopy; in particular, to infrared spectroscopy via dipole autocorrelation functions. The performance of both methods for computing vibrational spectra of several simple but representative molecular model systems is investigated systematically as a function of temperature and isotopic substitution. In this context both CMD and RPMD feature intrinsic problems which are quantified and investigated in detail. Based on the obtained results guidelines for using CMD and RPMD to compute infrared spectra of molecular systems are provided.

Microsolvation of protonated methane: structures and energetics of CH5(+)(H2)n.

Effects of microsolvating CH(5)(+) with up to four H(2) molecules have been investigated in terms of structures and energies. For the smaller complexes, benchmark calculations have been carried out using MP2 and CCSD(T) with basis sets up to aug-cc-pV5Z quality and energies have been extrapolated to the infinite basis set limit. It is found that MP2 calculations using the aug-cc-pVQZ basis set or better yield robust reference data for both structures and energies. More than 30 stationary points including minima and first-order as well as second-order stationary points have been characterized by this method and are discussed in terms of solvation motifs. Finally, the performance of several density functionals has been assessed for this very demanding case. Popular GGA functionals such as BLYP and PBE fail, whereas the TPSS meta-GGA functional captures many structural and energetic aspects of microsolvation satisfactorily.

Simulations of one- and two-electron systems by Bead-Fourier path integral molecular dynamics.

The Bead-Fourier path integral molecular dynamics technique introduced earlier [S. D. Ivanov, A. P. Lyubartsev, and A. Laaksonen, Phys. Rev. E 67 066710 (2003)] is applied for simulation of electrons in the simplest molecules: molecular hydrogen, helium atom, and their ions. Special attention is paid to the correct description of electrons in the core region of a nucleus. In an attempt to smooth the Coulomb potential at small distances, a recipe is suggested. The simulation results are in excellent agreement with the analytical solution for the "harmonic helium atom", as well as with the vibrational potential of the H(2) molecule and He ionization energies. It is demonstrated, that the Bead-Fourier path integral molecular dynamics technique is able to provide the accuracy required for the description of electron structure and chemical bonds in cases when electron exchange effects need not be taken into account.

Bead-Fourier path integral molecular dynamics.

Molecular dynamics formulation of Bead-Fourier path integral method for simulation of quantum systems at finite temperatures is presented. Within this scheme, both the bead coordinates and Fourier coefficients, defining the path representing the quantum particle, are treated as generalized coordinates with corresponding generalized momenta and masses. Introduction of the Fourier harmonics together with the center-of-mass thermostating scheme is shown to remove the ergodicity problem, known to pose serious difficulties in standard path integral molecular dynamics simulations. The method is tested for quantum harmonic oscillator and hydrogen atom (Coulombic potential). The simulation results are compared with the exact analytical solutions available for both these systems. Convergence of the results with respect to the number of beads and Fourier harmonics is analyzed. It was shown that addition of a few Fourier harmonics already improves the simulation results substantially, even for a relatively small number of beads. The proposed Bead-Fourier path integral molecular dynamics is a reliable and efficient alternative to simulations of quantum systems.