Communication: On the consistency of approximate quantum dynamics simulation methods for vibrational spectra in the condensed phase

The simulation of nuclear quantum dynamics in the condensed phase remains a challenge even for state of the art methods. The main problem is that the dimensionality of these systems is too large for a fully “first-principles” quantum calculation of the nuclear dynamics, leaving the whole field without a reliable benchmark. Therefore, it is hard to assess how well different approximations perform even for exhaustively studied systems like bulk water. There is a dire need to find a reliable framework in which the quantum dynamics of these systems can be treated, so that the field can advance further.

One can identify at least two classes of popular methods to model quantum dynamics in the condensed phase. The first treats a subset of the system's degrees of freedom fully quantum mechanically, embedded in the electrostatic environment of the rest of the system. Examples of these methods are the local monomer approximation (LMon),^1,2 which treats quantum mechanically a small subset of all intramolecular degrees of freedom in a static environment, and a mixed quantum-classical (MQC) model based on the semi-classical theory of line shape that treats a single degree of freedom (e.g., an OH stretch) “on-the-fly” in the presence of perturbing electric field generated by the surrounding classical solvent.^3–5 The second class of methods treats the whole system on the same footings, but using (physically motivated) ad hoc approximations to the quantum dynamics. These methods are often based on classical trajectories initiated from quantum distributions,^6–8 or are inspired by imaginary time path integral molecular dynamics (PIMD), as in the cases of centroid molecular dynamics (CMD)^9,10 and (thermostatted) ring polymer molecular dynamics (TRPMD).^11–13 Since these two classes of approaches are grounded on very different approximations, a comparison between their results is perhaps the closest one can get to assess their reliability in the absence of an absolute, exact benchmark.^14,15

Here, we perform such cross-validation, comparing the behavior of MQC, LMon, CMD, and TRPMD when modeling a physical observable that is particularly sensitive to the description of nuclear quantum dynamics—the infra-red absorption spectra of hydrogen-containing compounds.^{13,14,16–18} We focused on prototypical hydrogen-bonded condensed phase systems—HOD in D₂O and neat water—describing them with the inexpensive q-TIP4P/F potential,¹⁹ and using the corresponding linear dipole moment surface (DMS) to evaluate the IR absorption. The rationale for this simplistic choice is that we are not so much interested in comparing with experimental data, but only to nail down differences between methods for condensed-phase quantum dynamics. The simplicity of the potential and the linear DMS make it possible to guarantee thorough statistical sampling and to focus on the problem of modeling the dynamics on a complex, anharmonic potential energy surface, rather than on the fact that CMD and (T)RPMD are much harder to justify when computing correlation functions of nonlinear operators.^20–22 It also allows us to test extensively the dependence of TRPMD, CMD, and LMon on the precise details of the calculations. To validate the two philosophies over a broad range of conditions, and draw conclusions that are not purely anecdotal, we considered three very different thermodynamic regimes: ice Ih at 150 K, liquid water at 300 K, and the hot liquid at 600 K at the experimental liquid/vapor coexistence density. To assess the importance of these approximate nuclear quantum dynamics, we also compare our results to classical MD simulations.

CMD and TRPMD simulations were performed using the i-PI code,²³ and LAMMPS as the force back-end.²⁴ The details of the q-TIP4P/F potential are available elsewhere,¹⁹ as well as the technical details of the integration and thermostatting of PIMD.²⁵ We also performed classical MD reference calculations, using a time step of 0.5 fs and a weak global thermostat.²⁶ For both TRPMD and CMD we applied a Langevin thermostat to the internal modes of the ring polymer. In TRPMD, we set the friction coefficients equal to the free ring-polymer frequencies: γ_k = ω_k. We made this choice as a compromise between the damping of spurious peaks in the spectrum and the broadening of the physical peaks,¹³ and used a time step of 0.25 fs. In our CMD calculations we scaled the dynamical masses of the internal modes of the ring polymer to shift all the frequencies to the common value Ω = 16 000 cm⁻¹, and used a time step of 0.025 fs, which makes CMD simulations ten times more demanding than TRPMD for the same level of statistics. CMD has better sampling efficiency than plain RPMD due to the thermostatting of the ring polymer normal modes,²⁷ but the efficiency of TRPMD should be comparable since it uses an analogous thermostat. In running the CMD calculations we noticed that the details of the thermostat have a significant impact on the centroid dynamics. As discussed in the supplementary material,²⁸ the commonly-adopted prescription of optimally coupling either Langevin or Nosé-Hoover thermostats with the shifted frequency interferes with adiabatic decoupling at least up to Ω ⩽ 16 000 cm⁻¹ in the cases studied here. We found that in the limit of weak thermostat coupling (λ = 0.01), proper adiabatic separation can be reached already at Ω < 8000 cm⁻¹, which could probably enable using longer time steps. This should be considered when performing CMD calculations and when comparing with existing results.

All our LMon calculations were performed using the MULTIMODE code.^29–31 The LMon method assumes it is possible to take a snapshot configuration of a condensed-phase system, and solves exactly the Schrödinger equation for a subset of its degrees of freedom, evaluating the changes of the potential as a function of deformations of a molecule along a subset of its “lormal” modes. A detailed explanation of this approach can be found, for instance, in Ref. 1. The LMon approximation makes it possible to progressively include more physics by increasing the dimensionality of the quantum subspace. Using just the three intra-molecular “lormal” modes as the quantum-mechanical subspace (LMon-3) is generally sufficient to describe the stretching and bending regions of the spectrum of water.³² A recent study of vibrational energy relaxation (VER) of dilute HOD in ice Ih considered also three additional intermolecular modes (LMon-6),³³ obtaining a description of the intermolecular vibrational spectrum and an improved accuracy in the intra-molecular region, but at a much higher computational cost. Here, we are mostly interested in the strongly-quantized intra-molecular vibrations, and so we attempted a simplified approach that uses a total of four modes (LMon-4). Hundreds of monomer+environment configurations were obtained from individual replicas of a series of PIMD simulations. For each snapshot, the three intramolecular modes are complemented with the three highest-frequency inter-molecular modes, that are however considered one at a time: for each monomer, 3 sets of LMon-4 calculations are performed separately. A smooth spectrum is obtained by combining Gaussian line shapes centered around the energies of transitions from the ground state to each of the excited states, weighed by the transition dipole moment. Further numerical details about all the calculations can be found in the supplementary material.²⁸ 

MQC calculations were performed with an in-house code. Following Ref. 14, several independent CMD trajectories were performed starting from equilibrated PIMD configurations, from which the time-dependence of the vibrational frequency and transition dipole moment for the OH stretch were computed. The semiclassical vibrational line shape was then evaluated as discussed in Ref. 34, including non-Condon corrections as defined, e.g., in Ref. 35. A more detailed description of the simulations, as well as a discussion of the importance of non-Condon corrections, is reported in the supplementary material.²⁸ From the point of view of computational cost, both MQC and LMon-4 are relatively inexpensive, the most demanding aspect being the preliminary (PI)MD simulation. Traditionally the environment is sampled classically—which particularly at low temperature can be orders of magnitude less demanding—but as we will see results are influenced by the classical or quantum treatment of the environment. Then, the evaluation of the stretching frequency in MQC is trivial, for such a simple potential. Each LMon-3 calculation takes a few seconds on a single core, and about 1 min for LMon-4. The cost would increase steeply as a function of the size of the quantum subspace, but up to LMon-6 the cost will most likely be dominated by the preliminary sampling of bulk configurations.

Let us start by considering a single HOD molecule solvated in D₂O. This system (and its counterpart, HOD in H₂O) is commonly used as a probe for studying the local structure of liquid water theoretically and experimentally.^14,36 The OH stretch is dynamically uncoupled with the environment, making this system well-suited for both MQC and LMon calculations.³³ 

The results for classical MD, CMD, TRPMD, MQC, and LMon are shown in Figure 1. Starting from the highest temperature, 600 K, with a hot and compressed liquid, we observe that classical MD is blue-shifted by around 50 cm⁻¹ with respect to all quantum methods. Even at this high temperature nuclear quantum effects affect the calculation of dynamical properties. The blue shift of the classical simulation gets more pronounced by lowering the temperature, as expected, amounting to approximately 100 cm⁻¹ at 150 K. At 600 K all quantum methods show a large blue-shift, increased band width and asymmetric or structured line shapes, compared to those at 150 and 300 K. The LMon-4 band intensity falls-off faster than other line shapes at the high-frequency edge of the band, and shows a sharper maximum which is 30–50 cm⁻¹ to the “blue” of the CMD and TRPMD maxima, which are not sharp. Considering that the influence of low-frequency intermolecular modes, e.g., hot bands, additional dipole variation, vibrational relaxation, etc., becomes more important with increasing temperature and that LMon-4 is only partially describing this coupling, the limitations in the LMon- 4 theory are also expected to be become more apparent. The importance of inter-molecular couplings and of the dynamics of the environment at the higher temperature is also suggested by the large difference we observe between the MQC line shape, the Condon approximation to the MQC line shape, and MQC vibrational density of states (VDOS, see the supplementary material²⁸).

At 300 K the agreement between different quantum techniques is perhaps even better. MQC and CMD overlap almost perfectly. The LMon peak is some 20–40 cm⁻¹ higher in frequency, and the difference between the LMon and MQC line shapes is much less dramatic than at 600 K. At 300 K, inhomogeneous broadening effects (i.e., effects not related to vibrational relaxation or motional narrowing) are dominant. The TRPMD peak is further blue-shifted (by 20–40 cm⁻¹, depending on whether one considers the maximum or the mean position of the peak) and artificially broadened by the strong thermostatting of non-centroid modes.

In ice Ih at 150 K TRPMD and LMon still are in remarkable—although perhaps fortuitous—agreement. CMD shows a pronounced red shift of 180 cm⁻¹, which should probably be attributed to the curvature problem.^13,17 The MQC peak is red-shifted by about 60 cm⁻¹ relative to LMon and TRPMD. One possible explanation for this shift is the description of the environment based on centroid configurations that are very close to classical. Contrary to higher temperatures, the LMon spectrum also depends on whether the environment configurations are obtained from the beads of the PIMD simulation, or from the centroid. In the latter case, as when using configurations from classical MD, the OH stretch peak is red-shifted by 25 cm⁻¹ (see the supplementary material²⁸). It is arguable which choice is more physically justified. On one hand, bead positions provide a statistically accurate snapshot of the quantum environment. On the other hand, centroids can be seen as a mean-field average of the quantum fluctuations of the neighboring molecules, closer perhaps to the spirit of a quantum-classical model. In the absence of a rigorous justification we can see this discrepancy as a sign of the break-down of the classical model for the environment, and consider the difference between the two spectra as an estimate of the reliability of LMon (and MQC) in this low temperature regime.

Having compared all methods for the HOD:D₂O benchmark, we consider in Fig. 2 the spectra of water in the three thermodynamic regimes discussed above. We did not perform MQC simulations, since the isolated chromophore assumption is less justified in H₂O. Let us start by discussing the OH band that shows similar trends to those observed for HOD in heavy water. At 600 K, CMD and TRPMD are in near-perfect agreement, and the peak maxima are slightly red-shifted relative to LMon-4 by about 50 cm⁻¹. This is consistent with our observations in HOD case, and the discrepancy can be attributed to the lack of homogeneous broadening and dynamical couplings in LMon. At 300 K TRPMD and LMon-4 are close to each other, while CMD shows a small red shift of less than 50 cm⁻¹, which is consistent with what we observed in HOD:D₂O. In the case of low-temperature ice the CMD peak is red shifted by 150 cm⁻¹ compared to both LMon-4 and TRPMD that agree well with each other—even though the TRPMD peak is considerably broader. The intensity of the stretching band is much under-estimated relative to experiments. This is due to the linear DMS of q-TIP4P/F, as evidenced by comparison with more sophisticated models,³⁷ but is irrelevant for our comparison of approximate methods for quantum dynamics.

Moving on to the bend, we observe good agreement in the peak position and width among all the quantum methods. All predict a slight red shift with increasing temperature (in qualitative agreement with experiment). TRPMD and CMD give peaks positions at 1640 cm⁻¹ at 150 K and about 1610 cm⁻¹ in at 600 K. LMon-4 gives 1640 cm⁻¹ at 150 K and 1590 at 600 K. The peak at 600 K is further lowered to 1570 cm⁻¹ when using LMon-3 (see the supplementary material²⁸), indicating the growing importance of inter-molecular coupling at high temperature, and suggesting that this slight discrepancy between LMon-4 and PIMD-based methods could be resolved by increasing further the dimensionality of the quantum subspace. Finally, in the low-frequency region CMD and TRPMD are almost identical. In this region one cannot expect LMon-4 to yield quantitative accuracy. The quantum sub-space does not contain the collective modes of the hydrogen-bond network, nor the translation modes of individual monomers.

Within the LMon scheme it is straightforward to treat effects beyond linear absorption, and to give a clear physical attribution of specific features of the spectrum. For instance, all of the spectra display a distinct bump or shoulder around 3200 cm⁻¹, which corresponds to the first overtone of the bend. An interesting feature captured by LMon-4 (but not by LMon-3) is the small peak around 2350 cm⁻¹, which is evident at 150 K and becomes less clear-cut with raising temperature. This feature is due to the combination band of bending and librational modes, and demonstrates how increasing the dimensionality of the quantum subspace progressively includes additional physical effects in LMon calculations. Both CMD and TRPMD display non-zero absorption in this region, but further analyses would be needed to attribute that spectral density to a precise physical origin.

The conclusions of our comparison of approximate quantum dynamics simulation methods are overall optimistic: all the methods we considered are generally consistent with each other, while the difference with the position of the OH stretching peak observed in classical molecular dynamics confirms the importance of including nuclear quantum effects to reproduce quantitatively spectroscopic measurements of hydrogen-containing systems. TRPMD and LMon agree within a few tens of cm⁻¹ over a range of thermodynamic conditions going from ice Ih at 150 K to water at 300 K and finally to the hot, compressed liquid at 600 K. So do CMD down to room temperature, and MQC in the cases where we could apply it. Furthermore, this analysis reinforces the notion that performing PIMD-based and quantum-subspace simulations in tandem does not only provide a degree of cross-validation, but also makes it possible to profit simultaneously from the complete (albeit approximate) description of the absorption spectrum given by the former family of methods, and from the interpretation of distinct spectral features that is enabled by an explicit quantum treatment.³⁸ 

The comparison between the different methods also highlights their current limitations, and suggests that results of quantum dynamics should always be interpreted with caution. Extensions of MQC beyond the case of an isolated chromophore require one to model the interaction between the different chromophores, either explicitly or by devising an effective mapping of the molecular Hamiltonian on a collective coordinate, such as the local electrostatic potential.³⁹ Even though LMon can in principle be improved systematically, its computational cost increases exponentially with the size of the quantum subspace, making it challenging to fully assess the importance of inter-molecular couplings. Treating the environment semiclassically, in the same spirit as the MQC model, might be a better way to account for inter-molecular dynamical effects. One has also to consider that the present work represents the first application of LMon to a finite-temperature scenario. More work is needed to include systematically a description of monomer-monomer couplings, hot bands and the dynamical nature of the environment. Particularly at the highest temperature, these effects may affect the line shape of the absorption peaks, as highlighted by the comparison with MQC. We also observe that at 150 K the results of a LMon calculation depends sizably on whether one takes snapshots of the environment that are consistent with quantum statistics (e.g., from the beads in an PIMD simulation) or that are essentially classical (e.g., taken from classical MD, or using the ring polymer centroids).

Methods inspired by PIMD also have their issues. In its straightforward partially-adiabatic implementation, CMD requires a small integration time step, and one should pay attention to the fine details of the thermostatting to avoid interfering with the adiabatic decoupling of the centroid. At low temperatures, a pronounced red-shift relative to all other methods is observed. The precise temperature where this artifact becomes sizable should depend on the details of the system being studied. TRPMD, on the other hand, can be straightforwardly used with the same time step of a PIMD calculation, and does not suffer from systematic low-temperature red shifts. However, the OH stretch mode can be up to 50 cm⁻¹ blue-shifted relative to some of the other techniques. Furthermore, the precise peak position varies depending on the choice of the thermostat coupling, which is somewhat arbitrary. The strong thermostatting that is used to fix the resonance problems of RPMD also leads to a significant broadening of the peaks, so the TRPMD line shape should not be taken too seriously, unless it is dominated by inhomogeneous broadening or is weakly quantized.

Our results give some confidence on the reliability of approximate quantum dynamical methods for simulating the dynamics and the spectroscopy of condensed-phase systems or large molecules. The difference between quantum and classical simulations is comparable to the errors due to imperfect potentials, but the discrepancies between different quantum methods are less important. The best agreement between the various methods is seen for liquid water at room temperature—perhaps unsurprisingly considering that it corresponds to the regime for which most of approximate quantum techniques were designed. Cross-validation between different approaches, particularly at thermodynamic conditions outside of their “comfort zone,” offers the best promise for demonstrating the accuracy of existing methods, for uncovering their errors and shortcomings, and for delivering more reliable, simple, and affordable simulation strategies.

We gratefully acknowledge stimulating discussion with David Manolopoulos, Gregory Voth, and Lu Wang. M.R. thanks the Deutsche Forschungsgemeinschaft (DFG), RO 4637/1, for funding. H.L. and J.M.B. thank the U.S. National Science Foundation, CHE-1145227, for funding. F.P. thanks the U.S. National Science Foundation, CHE-1111364, for funding, and the Extreme Science and Engineering Discovery Environment (OCI-1053575, TG-CHE110009) for computational time.