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

Path Integral Simulations of Condensed-Phase Vibrational Spectroscopy

Recent theoretical and algorithmic developments have improved the accuracy with which path integral dynamics methods can include nuclear quantum effects in simulations of condensed-phase vibrational spectra. Such methods are now understood to be approximations to the delocalized classical Matsubara dynamics of smooth Feynman paths, which dominate the dynamics of systems such as liquid water at room temperature. Focusing mainly on simulations of liquid water and hexagonal ice, we explain how the recently developed quasicentroid molecular dynamics (QCMD), fast-QCMD, and temperature-elevated path integral coarse-graining simulations (Te PIGS) methods generate classical dynamics on potentials of mean force obtained by averaging over quantum thermal fluctuations. These new methods give very close agreement with one another, and the Te PIGS method has recently yielded excellent agreement with experimentally measured vibrational spectra for liquid water, ice, and the liquid-air interface. We also discuss the limitations of such methods.

Operator Formulation of Feynman Path Centroid Dynamics for Rotations

An operator formulation of centroid molecular dynamics (CMD) for rotational degrees of freedom is presented. The quasi-density operator concept was introduced by Jang and Voth [J. Chem. Phys 111, 2357 (1999)] and is used to obtain a phase-space mapping without the need for discretized path integrals. The approach allows the calculation of approximate Kubo-transformed time correlation functions. The particle on a ring is chosen as an illustrative example. Numerical results demonstrate that the proposed approach leads to accurate results when compared with exact diagonalization calculations for linear operators. At very low temperatures, it is found that rotational CMD yields results that are in exact agreement with the quantum dynamics of a spin-1 system.

Vibrational Spectra of Highly Anharmonic Water Clusters: Molecular Dynamics and Harmonic Analysis Revisited with Constrained Nuclear-Electronic Orbital Methods

Vibrational spectroscopy is widely used to gain insights into structural and dynamic properties of chemical, biological, and materials systems. Thus, an efficient and accurate method to simulate vibrational spectra is desired. In this paper, we justify and employ a microcanonical molecular simulation scheme to calculate the vibrational spectra of three challenging water clusters: the neutral water dimer (H4O2), the protonated water trimer (H7O3+), and the protonated water tetramer (H9O4+). We find that with the accurate description of quantum nuclear delocalization effects through the constrained nuclear-electronic orbital framework, including vibrational mode coupling effects through molecular dynamics simulations can additionally improve the vibrational spectrum calculations. In contrast, without the quantum nuclear delocalization picture, conventional ab initio molecular dynamics may even lead to less accurate results than harmonic analysis.

Fast Quasi-Centroid Molecular Dynamics for Water and Ice

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.

Fully First-Principles Surface Spectroscopy with Machine Learning

Our current understanding of the structure and dynamics of aqueous interfaces at the molecular level has grown substantially due to the continuous development of surface-specific spectroscopies, such as vibrational sum-frequency generation (VSFG). As in other vibrational spectroscopies, we must turn to atomistic simulations to extract all of the information encoded in the VSFG spectra. The high computational cost associated with existing methods means that they have limitations in representing systems with complex electronic structure or in achieving statistical convergence. In this work, we combine high-dimensional neural network interatomic potentials and symmetry-adapted Gaussian process regression to overcome these constraints. We show that it is possible to model VSFG signals with fully ab initio accuracy using machine learning and illustrate the versatility of our approach on the water/air interface. Our strategy allows us to identify the main sources of theoretical inaccuracy and establish a clear pathway toward the modeling of surface-sensitive spectroscopy of complex interfaces.

Vibrational strong coupling in liquid water from cavity molecular dynamics

We assess the cavity molecular dynamics method for the calculation of vibrational polariton spectra using liquid water as a specific example. We begin by disputing a recent suggestion that nuclear quantum effects may lead to a broadening of polariton bands, finding instead that they merely result in anharmonic red shifts in the polariton frequencies. We go on to show that our simulated cavity spectra can be reproduced to graphical accuracy with a harmonic model that uses just the cavity-free spectrum and the geometry of the cavity as input. We end by showing that this harmonic model can be combined with the experimental cavity-free spectrum to give results in good agreement with optical cavity measurements. Since the input to our harmonic model is equivalent to the input to the transfer matrix method of applied optics, we conclude that cavity molecular dynamics cannot provide any more insight into the effect of vibrational strong coupling on the absorption spectrum than this transfer matrix method, which is already widely used by experimentalists to corroborate their cavity results.

Path integral Brownian chain molecular dynamics: A simple approximation of quantum vibrational dynamics

An approximate approach to quantum vibrational dynamics, "Brownian chain molecular dynamics (BCMD)," is proposed to alleviate the chain resonance and curvature problems in the imaginary time-based path integral (PI) simulation. Here the non-centroid velocity is randomized at each step when solving the equation of motion of path integral molecular dynamics. This leads to a combination of the Newton equation and the overdamped Langevin equation for the centroid and non-centroid variables, respectively. BCMD shares the basic properties of other PI approaches such as centroid and ring polymer molecular dynamics: It gives the correct Kubo-transformed correlation function at short times, conserves the time symmetry, has the correct high-temperature/classical limits, gives exactly the position and velocity autocorrelations of harmonic oscillator systems, and does not have the zero-point leakage problem. Numerical tests were done on simple molecular models and liquid water. On-the-fly ab initio BCMD simulations were performed for the protonated water cluster, H 5 O 2 + $$ {\mathrm{H}}_5{\mathrm{O}}_2^{+} $$ , and its isotopologue, D 5 O 2 + $$ {\mathrm{D}}_5{\mathrm{O}}_2^{+} $$ .

Quantitative molecular simulations

All-atom simulations can provide molecular-level insights into the dynamics of gas-phase, condensed-phase and surface processes. One important requirement is a sufficiently realistic and detailed description of the underlying intermolecular interactions. The present perspective provides an overview of the present status of quantitative atomistic simulations from colleagues' and our own efforts for gas- and solution-phase processes and for the dynamics on surfaces. Particular attention is paid to direct comparison with experiment. An outlook discusses present challenges and future extensions to bring such dynamics simulations even closer to reality.

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 (H₂O(l)), heavy water (D₂O(l)), T₂O(l), HDO(aq) and HTO(aq) studied by ab initio path integral molecular dynamics (PIMD). Unlike previous ab initio PIMD investigations of H₂O(l) and D₂O(l) [Chen et al., Phys. Rev. Lett., 2003, 91, 215503] [Machida et al., J. Chem. Phys., 2017, 148, 102324] we find that D₂O(l) is more structured than H₂O(l), as is predicted by the experiment. The agreement between the experiment and our simulation for H₂O(l) and D₂O(l) allows us to accurately predict the intra- and intermolecular structures of T₂O(l) HDO(aq) and HTO(aq). T₂O(l) is found to have a similar intermolecular structure to that of D₂O(l), while the intramolecular structure is more compact, giving rise to a smaller dipole moment than those of H₂O(l) and D₂O(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 H₂O and D₂O. 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.

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.

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.

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.

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.

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.

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.

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.

Combined Effects from Solvation and Nuclear Quantum Fluctuations on Autoionization Mechanisms in Aqueous Clusters

Using path-integral molecular dynamics simulations, we examine isomerization paths involving collective proton transfers in [H₂O]₅ and [H₂O]₈ clusters under cryogenic conditions. We focused attention on combined effects derived from solvation and nuclear quantum fluctuations on the characteristics of free energy barriers and relative stabilities of reactants and products. In particular, we analyzed two different processes: the first one involves the exchange of donor-acceptor hydrogen bond roles along cyclic moieties, whereas the second one corresponds to charge separation leading to stable [H₃O]⁺[OH]^- ion pairs. In the first case, the explicit incorporation of quantum tunneling introduces important modifications in the classical free energy profile. The resulting quantum profile presents two main contributions, one corresponding to compressions of O-O distances and a second one ascribed to nuclear tunneling of the light protons. Solvation effects promote a moderate polarization of the cyclic structures and a partial loss of concertedness in the collective modes, most notably, at the onset of tunneling. Still, the latter effects are also sufficiently strong to promote the stabilization of ion pairs along the classical trajectories. In contrast, the explicit incorporation of nuclear quantum fluctuations leads to charge separated configurations that are marginally stable. As such, the latter states could also be regarded as short-lived intermediate states along the reactive exchange path.

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.

Probing photodissociation dynamics using ring polymer molecular dynamics

The performance of the ring polymer molecular dynamics (RPMD) approach to simulate typical photodissociation processes is assessed. The correct description of photodissociation requires the calculation of correlation functions or expectation values associated with non-equilibrium initial conditions, which was shown to be possible with RPMD very recently [J. Chem. Phys. 145, 204118 (2016)]. This approach is combined with treatment of the nonadiabatic dynamics employing the ring polymer surface hopping approach (RPSH), which is based on Tully's fewest switches surface hopping (FSSH) approach. The performance is tested using one-dimensional photodissociation models. It is found that RPSH with non-equilibrium initial conditions can well reproduce the time-dependent dissociation probability, and adiabatic and diabatic populations for cases where the crossing point is below and above the Franck-Condon point, respectively, while standard FSSH fails to reproduce the exact quantum dynamics in the latter case. Thus, it is shown that RPSH is an efficient and accurate alternative to standard FSSH, which is one of the most widely employed approaches to study photochemistry.

Approximating Matsubara dynamics using the planetary model: Tests on liquid water and ice

Matsubara dynamics is the quantum-Boltzmann-conserving classical dynamics which remains when real-time coherences are taken out of the exact quantum Liouvillian [T. J. H. Hele et al., J. Chem. Phys. 142, 134103 (2015)]; because of a phase-term, it cannot be used as a practical method without further approximation. Recently, Smith et al. [J. Chem. Phys. 142, 244112 (2015)] developed a "planetary" model dynamics which conserves the Feynman-Kleinert (FK) approximation to the quantum-Boltzmann distribution. Here, we show that for moderately anharmonic potentials, the planetary dynamics gives a good approximation to Matsubara trajectories on the FK potential surface by decoupling the centroid trajectory from the locally harmonic Matsubara fluctuations, which reduce to a single phase-less fluctuation particle (the "planet"). We also show that the FK effective frequency can be approximated by a direct integral over these fluctuations, obviating the need to solve iterative equations. This modification, together with use of thermostatted ring-polymer molecular dynamics, allows us to test the planetary model on water (gas-phase, liquid, and ice) using the q-TIP4P/F potential surface. The "planetary" fluctuations give a poor approximation to the rotational/librational bands in the infrared spectrum, but a good approximation to the bend and stretch bands, where the fluctuation lineshape is found to be motionally narrowed by the vibrations of the centroid.

First-Principles Simulations of Liquid Water Using a Dielectric-Dependent Hybrid Functional

We carried out first-principles simulations of liquid water under ambient conditions using a dielectric-dependent hybrid functional, where the fraction of exact exchange is set equal to the inverse of the high-frequency dielectric constant of the liquid. We found excellent agreement with experiment for the oxygen-oxygen partial correlation function at the experimental equilibrium density and 311 ± 3 K. Other structural and dynamical properties, such as the diffusion coefficient, molecular dipole moments, and vibrational spectra, are also in good agreement with experiment. Our results, together with previous findings on electronic properties of the liquid with the same functional, show that the dielectric-dependent hybrid functional accurately describes both the structural and electronic properties of liquid water.

IR Spectra of (HCOOH)2 and (DCOOH)2: Experiment, VSCF/VCI, and Ab Initio Molecular Dynamics Calculations Using Full-Dimensional Potential and Dipole Moment Surfaces

We report quantum VSCF/VCI and ab initio molecular dynamics (AIMD) calculations of the IR spectra of (HCOOH)₂ and (DCOOH)₂, using full-dimensional, ab initio potential energy and dipole moment surfaces (PES and DMS). These surfaces are fits, using permutationally invariant polynomials, to 13 475 ab initio CCSD(T)-F12a electronic energies and MP2 dipole moments. Here "AIMD" means using these ab initio potential and dipole moment surfaces in the MD calculations. The VSCF/VCI calculations use all (24) normal modes for coupling, with a four-mode representation of the potential. The quantum spectra align well with jet-cooled and room-temperature experimental spectra over the spectral range 600-3600 cm^-1. Analyses of the complex O-H and C-H stretch bands are made based on the mixing of the VSCF/VCI basis functions. The comparisons of the AIMD IR spectra with both experimental and VSCF/VCI ones provide tests of the accuracy of the AIMD approach. These indicate good accuracy for simple bands but not for the complex O-H stretch band, which is upshifted from experimental and VSCF/VCI bands by roughly 300 cm^-1. In addition to testing the AIMD approach, the PES, DMS, and VSCF/VCI calculations for formic acid dimer provide opportunities for testing other methods to represent high-dimensional data and other methods that perform postharmonic vibrational calculations.

The strengths and limitations of effective centroid force models explored by studying isotopic effects in liquid water

The development of effective centroid potentials (ECPs) is explored with both the constrained-centroid and quasi-adiabatic force matching using liquid water as a test system. A trajectory integrated with the ECP is free of statistical noises that would be introduced when the centroid potential is approximated on the fly with a finite number of beads. With the reduced cost of ECP, challenging experimental properties can be studied in the spirit of centroid molecular dynamics. The experimental number density of H2O is 0.38% higher than that of D2O. With the ECP, the H2O number density is predicted to be 0.42% higher, when the dispersion term is not refit. After correction of finite size effects, the diffusion constant of H2O is found to be 21% higher than that of D2O, which is in good agreement with the 29.9% higher diffusivity for H2O observed experimentally. Although the ECP is also able to capture the redshifts of both the OH and OD stretching modes in liquid water, there are a number of properties that a classical simulation with the ECP will not be able to recover. For example, the heat capacities of H2O and D2O are predicted to be almost identical and higher than the experimental values. Such a failure is simply a result of not properly treating quantized vibrational energy levels when the trajectory is propagated with classical mechanics. Several limitations of the ECP based approach without bead population reconstruction are discussed.

Quantum Dynamics and Spectroscopy of Ab Initio Liquid Water: The Interplay of Nuclear and Electronic Quantum Effects

Understanding the reactivity and spectroscopy of aqueous solutions at the atomistic level is crucial for the elucidation and design of chemical processes. However, the simulation of these systems requires addressing the formidable challenges of treating the quantum nature of both the electrons and nuclei. Exploiting our recently developed methods that provide acceleration by up to 2 orders of magnitude, we combine path integral simulations with on-the-fly evaluation of the electronic structure at the hybrid density functional theory level to capture the interplay between nuclear quantum effects and the electronic surface. Here we show that this combination provides accurate structure and dynamics, including the full infrared and Raman spectra of liquid water. This allows us to demonstrate and explain the failings of lower-level density functionals for dynamics and vibrational spectroscopy when the nuclei are treated quantum mechanically. These insights thus provide a foundation for the reliable investigation of spectroscopy and reactivity in aqueous environments.

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.