Chen
Qu
and
Joel M.
Bowman
*
Department of Chemistry and Cherry L. Emerson Center for Scientific Computation, Emory University, Atlanta, USA. E-mail: jmbowma@emory.edu
First published on 26th October 2018
The subject of this Perspective is quantum approaches, beyond the harmonic approximation, to vibrational dynamics and IR spectroscopy. We begin with a pedagogical, unifying review of the most widely used quantum approaches. Some of the key details that lead to steep computational scaling of these approaches are reviewed, as well as some effective strategies to overcome or at least mitigate them. Considering in particular the application to IR spectroscopy, we stress the strength and weakness of each approach for spectral features that evolve from “simple” to “complex”. We use the 10-atom formic acid dimer as an ideal example of this evolution. The IR spectrum of this dimer and two isotopologs has been obtained computationally using our software, MULTIMODE, and approaches to obtain accurate, ab initio, full-dimensional potential energy and dipole moment surfaces, also developed by our group. The IR spectra obtained with the widely used “ab initio molecular dynamics” approach are also presented and assessed. The extension of quantum approaches to molecular clusters and even condensed phase applications, where molecular dynamics approaches are typically used, is discussed mainly in the context of the local monomer model. This approach is illustrated for a methane clathrate hydrate, where vibrational energies of the symmetric and asymmetric stretches of methane are given for a number of water cages and compared to experiment. The question about interpretation is also addressed throughout the Perspective.
Although the study of vibrational and ro-vibrational dynamics is pervasive in chemical physics, it is fair to assert that IR spectroscopy is the major area where theory meets experiment. High-resolution ro-vibrational spectroscopy is of course an important part of this field; however, here we consider medium-resolution IR spectroscopy (i.e., not rotationally resolved). As such, the focus is band positions, intensities, bandwidths and band sub-structures. Challenges to theory to reproduce such features range from routine to very difficult. Sharp bands usually indicate that a single vibrationally excited eigenstate is dominantly responsible for the band. Often, a standard double-harmonic analysis is sufficient to assign the band; however, to improve quantitative accuracy vibrational perturbation theory is now typically employed. Sharp bands separated by a few wavenumbers may indicate a relatively simple resonance interaction, e.g., Fermi and/or Darlington–Dennison resonances. These are usually successfully treated by degenerate perturbation theory.^{3} The analysis of a complex band, that is a broad band with sub-structure, is more challenging both qualitatively and quantitatively. In this case, a vibrational configuration interaction (VCI) or vibrational coupled cluster (VCC) approach is almost certainly needed.
These approaches, starting with the harmonic one and ending with the VCI approach (which is “exact”, in principle), provide the critical assignment and interpretation of the spectral bands. While the harmonic approach provides the simplest interpretation, it may not be correct for bands that are inherently not captured by this approach, e.g., combination and overtone bands, etc. A more rigorous approach may provide acceptable agreement for a complex experimental band, but then the matter of interpretation may also become complex, as the title of this Perspective suggests. We endeavor to unravel this complexity, at least to some extent here.
To begin, we list the post-harmonic, quantum methods that are widely used for vibrational (and ro-vibrational) dynamics calculations below. Citations, which include some software, are not comprehensive of course. An overview of methods, that is somewhat dated but detailed, appeared in 2008,^{4} and several recent Perspective articles that are more focused are highly recommended.^{5–10}
• VPT2/QFF^{11–13} (vibrational second order perturbation theory with quartic force fields)
• VSCF^{5,14–17} (vibrational self-consistent field theory)
• VCI and VSCF/VCI^{4,5,7–10,18–25} (vibrational/virtual state configuration interaction theory)
• VSCF/PT2^{5,26} (VSCF + second-order perturbation theory)
• VSCF/VCC^{6,27} (VSCF + vibrational coupled-cluster theory)
• MCTDH^{28–30} (multiconfiguration time-dependent hartree theory)
The most rigorous, and thus computationally intensive, methods in this list are the VCI and MCTDH ones. These methods are used and are being developed by a number of groups, including ours, and with different objectives in mind. These range from developing general methods, using general potential energy surfaces (PESs) that can be applied to molecules with large-amplitude motion, perhaps across multiple minima to ones that are more restrictive in the molecular motion. The trade-off between these methods is essentially the size of the molecules that can be considered, with the former approach currently limited to 4–5 atoms. Our focus has been on the latter approach, using general potential energy surfaces, with the aim to extend it to molecules with 10 or more atoms. We use the Watson Hamiltonian, which is a general Hamiltonian in mass-scaled rectilinear normal modes. The trade-off in using this Hamiltonian is the use of rectilinear normal modes, which are not universally optimum for the broad class of molecular vibrations. They are generally better for high-frequency stretch modes than low-frequency, large amplitude ones, which often are better described using curvilinear coordinates.
By making computationally efficient, but still accurate, representations of the potential, our VCI approach (VSCF/VCI) has been applied to a number of molecules with more than 4 atoms and most recently to the IR spectra of the 10-atom formic acid dimer and isotopologs.^{31–35} We will review the salient aspects of this work below. We also applied the so-called AIMD approach, using the same potential and dipole moment surface and since this approach is very widely used, we also review the spectra obtained with it.
The AIMD approach is classical and so it misses a number of important quantum effects, such as zero-point energy (ZPE). Also, since it generally is applied “on-the-fly”, it must be used with a computationally efficient electronic method, and so DFT is the universal choice. Both the classical treatment of the nuclei and the DFT treatment of the electrons might raise some doubts about the quantitative accuracy of the approach or “getting the right answer for the right reason”. However, with many applications and extensive experience, this approach has evolved into a very useful post-harmonic model with wide applicability. Here post-harmonic has a different meaning relative to the quantum qualifier. Specifically, even at room temperature the classical dynamics may be large amplitude, in fact very large amplitude, where isomers or conformers may be sampled by the MD trajectories. These large amplitude motions generally correspond to low-frequency modes with the high-frequency modes, e.g., X–H stretches, “going along for the ride”. These high-frequency modes generally execute small-amplitude motion because the classical energy in these modes at say room temperature is very small, relative even to the zero-point energy of these stretches.
To go significantly beyond the 10-atom formic acid dimer and still use a VCI approach requires some major assumptions about the potential and/or the way modes are coupled. Here, we illustrate one such strategy, which we have termed the local monomer model. We have applied this model extensively to clusters, condensed-phase ice and liquid water,^{36–43} and several clathrate hydrates.^{44–46} This model will be illustrated by our previous work on the vibrational dynamics of CH_{4} encapsulated in clathrate hydrates.
Finally, at the risk of stating the obvious, it is worth noting that an exact quantum treatment of the nuclear motion with the exact potential will produce highly accurate results that should be in excellent agreement with experiment. The corollary statement is the use of an inaccurate potential with an exact treatment of the nuclear motion will lead to inaccurate results, which would not be in agreement with experiment. Of course, modern computational chemistry must contend with less than “exact”; however, for potentials the current “gold standard” is the CCSD(T) method with at least an aug-cc-pVTZ basis or using an explicitly correlated F12 basis. With this level of theory for the electrons, a very precisely fitted PES and an essentially exact treatment of the vibrational dynamics, experience tells us that agreement with experiment at the level of 5 cm^{−1} or less can be expected. A recent example illustrating this has been reported for DOCO.^{47} This level of accuracy can deteriorate, especially for high frequency X–H stretches if a lower level of electronic theory, e.g., MP2 or DFT is used instead of CCSD(T). The scaling of the computational effort for the CCSD(T) with the number of electrons is much steeper than say DFT method, such that CCSD(T) is infeasible to use in say an AIMD calculation, where DFT is feasible and is widely used for large-scale applications, as noted already. We will give some essential details of the AIMD approach later, but we re-state that since it is a classical approach it cannot be considered rigorous. So even if it used with a highly accurate potential energy surface there is no guarantee that the final results will agree well with experiment or with a quantum approach using the same PES; however, under certain conditions this is seen, as we hope to show for the formic acid dimer application below.
The paper is organized as follows. In the next section, we give a brief pedagogical review of several post-harmonic quantum methods, with the aim of making intellectual connections among them. The VSCF/VCI method, the workhorse of our group is ultimately focused on. Then, a brief discussion of the PES, the “elephant in the room”, is given. This is a critical part of any simulation and is an important aspect of method development. The particular representation of the potential that is central to our approach is described in Section 3. Following that, two examples of VSCF/VCI calculations are given. The IR spectra of the formic acid dimer and two isotopologs, and comparison with experiment are presented in Section 4. IR spectra from the AIMD approach are also given and compared with the quantum and experimental ones. The local monomer model and the application to methane clathrate hydrate is given in Section 5. That section includes comments on localized modes and finally some comments on going beyond MD simulations for the condensed phase. A summary and conclusions are given in Section 6.
H^{(0)}Φ_{K} = E^{(0)}_{K}Φ_{K}, H^{(0)} = T + V^{(0)}. | (1) |
(2) |
For a molecule/molecular cluster of N atoms with 3N − 6 vibrational modes, the zero-order potential is a sum of harmonic-oscillator ones and the wavefunction is a single product of zero-order one-mode functions, Φ_{K} = Π^{3N−6}_{i=1}ϕ_{ki}(Q_{i}). While this is high-dimensional, it is basically trivial because it is written in a simple product form. Note K is a row vector [k_{1}, k_{2},…k_{3N−6}]. Inserting these products into eqn (2) to represent the exact wavefunction has exponential growth with the upper limit for each mode. For example, if the upper limit is a modest value of 10, the total number of terms in the summation is 10^{3N−6}, which is already impractical for N = 4, i.e., a non-linear tetraatomic molecule.
Clearly, this approach cannot be used without major modifications and some highly effective ones will be briefly reviewed in the next section. But, this exponential scaling has motivated less rigorous approaches, which do not have such scaling. One is second-order perturbation theory and the other is self-consistent field theory; these theories have a long history in electronic structure theory. We briefly review these next.
(3) |
This version of VPT2 is now available in a number of popular electronic structure codes and is widely used. It is both relatively easy to implement and is often much more accurate (perhaps surprisingly) than the harmonic approximation, provided resonances are not an issue. One reason for the accuracy is that VPT2 produces the exact eigenvalues of the Morse potential, using a quartic expansion of that potential.^{49} This potential and the quartic approximation to it in terms of a variable Q that is zero at the minimum are given in Fig. 1. As seen, the quartic approximation fails for large positive values of Q. However, while this is of no consequence for VPT2 theory, it should be clear that this can become a numerical issue for a variational approach using a basis or discrete variable variational approach, if the quartic potential is used. Of course, using the full Morse potential in variational approaches does produce exact results. For other potentials, for example double-well potentials, VPT2/QFF can be highly inaccurate or simply not applicable.^{50,51} For these, the variational VCI approaches are generally necessary.
Before proceeding, it is important to note that the example of using a harmonic basis is mainly for pedagogical purposes and to be consistent with the widely used version of VPT2. Clearly, any complete orthonormal basis could be used. In fact, as described in the next section, a more sophisticated optimized basis is used in our code MULTIMODE, which makes use of the Watson Hamiltonian. This basis is described briefly along with comments about applications which are dominated by a single minimum. Applications with a large-amplitude mode that spans two equivalent minima can also be described with this Hamiltonian and an extension of the basis. Molecules with internal rotors or multiple minima that impact the vibrational dynamics are more challenging and special bases and specialized coordinates are generally needed. Such molecules are beyond the scope of this Perspective.
The VSCF approach dates from 1978^{14–16} and has been extensively used ever since. It is a general approach that can be used with any molecular Hamiltonian and potential; however, it is most typically used with normal coordinates. As will become clear, it is an efficient approach to describe anharmonic effects and mode-coupling, in the mean field sense. So, it can be applied to large molecules and this has been done in numerous interesting and important applications by Gerber and co-workers.^{5}
The general VSCF ansatz for an n-mode wavefunction for a general quantum state is the direct product Π_{i}ϕ_{ni}(Q_{i}), where the ϕ_{ni}(Q_{i}) are the vibrational modal functions. The variational principle applied to each modal, subject to the constraint that each has a unit norm, leads to the well-known SCF equations. For a 3-mode example, with a Hamiltonian, given by H = T_{1} + T_{2} + T_{3} + V(Q_{1},Q_{2},Q_{3}) these are
(4) |
For the ground vibrational state, n_{i} = 0 for all modes and each effective potential is for the ground state. These 1d Schroedinger equations give the variationally optimized ground vibrational state and in addition excited vibrational states that are termed the virtual states. These states, including the ground VSCF state, are orthonormal and can be used instead of the harmonic-oscillator ones as the CI basis for the rigorous expansion of the exact wavefunction. We have denoted this as a VCI basis; however, this acronym is often used to simply denote a vibrational CI.
This approach, which we have denoted as VSCF/VCI, has been used in collaboration with Stuart Carter and Nicholas Handy to develop the code MULTIMODE,^{19,52,53} which we describe briefly in the next section. The code has recently been applied to a VSCF/VCI calculation of the IR spectra of the formic acid dimer and two isotopologs and we presents some highlights of those calculations in Section 4. Previously unreported VSCF calculations of vibrational energies are also presented.
We will give more details of the VSCF/VCI approach used in MULTIMODE in the next section and we conclude this section with the general expression for the IR intensity from a quantum calculation of the molecular eigenstates:
(5) |
(6) |
Here, we briefly review a general approach for both, for the case where the internal coordinates are mass-scaled normal coordinates, Q. The corresponding Hamiltonian (for non-linear molecules), known as the Watson Hamiltonian,^{56} is given by
(7) |
In many applications, the vibrational angular momentum terms are neglected. This approximate Hamiltonian cannot lead to exact results, but the errors are generally not large, otherwise it would not be widely used. Nevertheless, for molecules with hydrogen atoms, errors can be of the order of 10 cm^{−1}. A particularly glaring case is the ground-state tunneling splitting in H_{3}O^{+}. Using the full Watson Hamiltonian, with an accurate PES, the calculated splitting is 46 cm^{−1},^{57} in good agreement with experiment, versus 21.9 cm^{−1} using the approximate Hamiltonian with the same PES.^{58} We do not neglect these terms in the MULTIMODE software.
Returning to the issue of exponential scaling, we begin with the evaluation of matrix elements. Two general strategies to deal with this have been proposed. One is the sum-of-products (SOP) already mentioned in the discussion of quartic force fields and the second is the n-mode representation of the potential V(Q). In the former, multidimensional quadratures of the V can be written as the sum of products of 1d-quadratures. This is an extremely effective, albeit not universal strategy. The example of the QFF is relevant here. In early vibrational CI calculations, this form of the SOP was used, e.g., in the code POLYMODE^{59} and ANHAR.^{60} However, re-inspection of Fig. 1, and recalling the accompanying text about the behavior of the quartic representation of the Morse potential, tells us that this approach can (and does) fail in the region where the QFF fails to describe the Morse potential in a CI approach. (Other SOP approaches that are not Taylor series expansions of the potential, can in principle overcome this, e.g., the “POTFIT” method in the MCTDH approach;^{28–30} however, with some inevitable numerical inaccuracies and limitations on the complexity of the true potential.)
A second efficient approach to represent the potential is the hierarchical n-mode representation (nMR) of the potential.^{52} This approach was independently proposed and termed the cut-high-dimensional-model-representation of the potential.^{61} In normal coordinates, this representation is given by^{52}
(8) |
This representation of the potential, up to 6MR, is used in MULTIMODE. The 4MR of the full-dimensional PES was used recently for the formic acid dimer. This demonstrates that this representation can be used effectively for 10-atom molecules, with 24 vibrational modes. Also, the n-mode representation will turn out to be useful for the second example of the methane clathrate hydrate, which illustrates extending vibrational CI methods to molecular clusters and condensed phase.
A detailed discussion of this representation has been given in numerous publications, and a sample of these are ref. 6, 19, 52, 53 and 62–64. Also, it is important to note the pioneering work by Jung and Gerber, who introduced the 2MR of the potential in 1996.^{65} Gerber and co-workers and others have applied this efficient representation to vibrational dynamics of biomolecules and clusters.^{5} Indeed the 2MR does provide useful results which are generally a significant improvement over 1MR. However, it has limitations that need to be kept in mind. While it can describe 2:1 Fermi resonances, it cannot describe 1 + 1:1 resonances; these require a 3MR of the potential. More generally, it has been shown numerous times in the literature that the 2MR of general potentials gives results which are not converged to within 10–100 cm^{−1} even for fundamentals.^{4} However, a 4MR typically closes the gap to within roughly 1–5 cm^{−1} (ref. 4, 54 and 66) and we typically use 4MR with an existing full-dimensional PES. Of course tests with a 5MR of the potential can be done, and this is always recommended.
Some comments on the limitations of rectilinear normal modes and thus the Watson Hamiltonian are in order before considering the further bottlenecks in VCI calculations. First, this is not necessarily a limitation for problems involving two equivalent minima, for example the tunneling splitting of H_{3}O^{+} mentioned already. In this case and analogous ones the reference geometry is the relevant saddle point separating the two minima and the rectilinear imaginary frequency normal mode is the large amplitude coordinate. The basis for the MULTIMODE calculations are numerical ones determined by relaxed 1d cuts through the multidimensional potential. Details can be found elsewhere;^{57,67} however, it is sufficient to note that these are not harmonic oscillator bases. For other problems, with say low energy torsional modes or internal rotors, a more general approach is needed. The reaction path version of MULTIMODE^{68} is able to describe these problems; however it is beyond the scope of this Perspective to delve into this class of molecular vibrational dynamics.
The second major bottleneck to all VCI calculations is the size of the H-matrix, which as noted already can scale exponentially with the number of vibrational modes. There are many strategies to deal with this. Basically, they all limit the size of the excitation space, with many schemes taken from electronic structure theory. For example, the excitation space can be limited by using the hierarchical scheme of single, double, triple, etc. excitations. MULTIMODE uses this among other schemes and can consider up to quintuple excitations. Other schemes to prune the CI basis have been suggested and the reader is referred to reviews^{4–10} for more details.
In the next section, we present selected results of recent VSCF/VCI and AIMD calculations of the IR spectra of the formic acid dimer and isotopologs.^{31–33,35}
Numerous calculations have been reported attempting to reproduce the ground-state tunneling splitting. These generally make use of reduced dimensionality or other approximations.^{82–98} Recently, major progress has been made in the calculation of this splitting in full dimensionality,^{99} using a recent full-dimensional ab initio PES reported by us.^{31} The calculated splitting, 0.014 cm^{−1}, is in good agreement with experimental value of roughly 0.015 cm^{−1}, and this demonstrates the accuracy of our PES. The PES is a permutationally invariant fit to roughly 13475 CCSD(T)-F12a electronic energies. Later we reported a full-dimensional MP2-based dipole moment surface. These were used in VSCF/VCI (and MD) calculations of the IR spectrum of FAD.^{32–35}
The experimental IR spectrum has been reported in jet-cooled and room temperature conditions, notably in the C–O, C–H and O–H stretch regions of the spectrum.^{69,70,72,78,80} These spectra inspired us to perform VSCF/VCI calculations of the IR spectra. We reported these recently in a series of papers, for FAD and two deuterated isotopologs.^{32,33,35} The VSCF/VCI spectra are in good agreement with the experimental ones. We also performed ab initio molecular dynamics calculations of the IR spectra. In the present case, “ab initio” means using the fitted potential and dipole moment surfaces instead of the more usual “on-the-fly” approach. Given that a single CCSD(T)-F12a energy calculation of the FAD using the haTZ basis (aVTZ for C and O, and VTZ for H) takes roughly 40 minutes on an 8-core cpu, the “on-the-fly” approach using this level of electronic theory is not feasible.
The complete set of 24 normal modes have been given in ESI of ref. 33. The mode of most interest here is the IR-active O–H stretch, shown in Fig. 2. As can be immediately surmised, this stretch is strongly perturbed by the double hydrogen bond.
Fig. 2 Formic acid dimer IR active O–H(D) stretch normal mode with indicated harmonic frequency. Reproduced from ref. 35 with permission from the Royal Society of Chemistry. |
As seen below, the experimental IR spectrum of FAD illustrates the full range of complexity of IR bands described in the Introduction. Therefore, reproducing these bands theoretically requires increasing rigor starting from the double harmonic approximation, to a full dimensional VSCF/VCI calculation. The VSCF/VCI calculations were done with the code MULTIMODE and are for J = 0. These calculations considered a single minimum, so they do not describe the proton transfer. Based on our published work, we concluded that this transfer has a negligible effect at the level of resolution of experiment and calculations.
To begin, we show the double-harmonic spectra of FAD and two indicated isotopologs in Fig. 3. We contrast this with the VSCF/VCI AIMD and experimental spectra shown in Fig. 4 and 5. As seen, the intense O–H stretch stick (at 3325 cm^{−1} for FAD) and the less intense C–H stretch stick (at 3097 cm^{−1} for FAD) in the double-harmonic spectra are complex fractionated bands in the experimental spectra. This fractionation is captured by the VSCF/VCI calculations and it suggests that there are many CI states that contribute to these bands. Histogram plots of the squares of relevant CI coefficients, c^{(L)}_{K}, cf.eqn (2), are shown in Fig. 6. We will refer to these as weights. The ones in blue correspond to coefficients associated with the O–H stretch fundamental basis function. This is the virtual state, but for discussion we can assume it is the harmonic wavefunction for the O–H stretch normal mode. The red histograms are the analogs for the C–H stretch. Only coefficients with non-negligible values are used in creating these histograms.
Fig. 3 Double harmonic spectra of formic acid dimer and indicated isotopologs. Reproduced from ref. 35 with permission from the Royal Society of Chemistry. |
Fig. 4 VSCF/VCI spectra of (HCOOH)_{2} and (DCOOH)_{2}, and comparison with jet-cooled (ref. 72 for (HCOOH)_{2} and ref. 70 for (DCOOH)_{2}) and room temperature (ref. 80) experiments. Reproduced from ref. 35 with permission from the Royal Society of Chemistry. |
Fig. 5 Comparison between the VSCF/VCI spectra of (HCOOH)_{2} and (DCOOH)_{2} and jet-cooled experiments (ref. 72 for (HCOOH)_{2} and ref. 70 for (DCOOH)_{2}), for the O–H stretch bands. Reproduced from ref. 35 with permission from the Royal Society of Chemistry. |
Fig. 6 Square of the VCI coefficients of the C–H(D) and O–H(D) stretches as a function of energy across the O–H stretch band for: (a) (HCOOH)_{2}, (b) (DCOOH)_{2}, and (c) (DCOOD)_{2}. Reproduced from ref. 35 with permission from the Royal Society of Chemistry. |
Consider the distribution of weights for the O–H stretch first. Starting with panel (a), we see numerous O–H weights that are highly dispersed over a large energy range, indeed roughly the same as the width of the O–H band in the IR spectrum. Single deuteration at the carbon atom has very little effect on the O–H stretch distribution of weights, as seen in panel (b); however, deuteration at the O does result in a narrower distribution of weights, i.e., less fractionation, as seen in panel (c). The situation for the C–H stretch is analogous to the O–H stretch, with the interesting exception of panel (b) where single deuteration at C results in a nearly pure C–D stretch. This is an indication that the C–H(D) stretch is perturbed by the O–H stretch and even the O–D stretch when the two fundamentals are close in energy. For the case in panel (b) they are well separated and thus not interacting.
These histograms of weights also bear on the variations in the IR intensity for the O–H band. This is easy to understand since the zero-order O–H fundamental is a strong absorber (cf.Fig. 3) and so the variations of the weight associated with that state, i.e., the blue histograms just discussed, translate directly into intensity variations seen in the VSCF/VCI spectra and experiment. A detailed analysis of these variations was reported at a 2018 Faraday discussions meeting.^{34}
Clearly, a vibrational CI approach is needed to describe these complex bands theoretically. This is an example of the strong-coupling case mentioned in the Introduction. A detailed analysis of the mode coupling has been given elsewhere,^{32,33} and the take-home message is that excitation of the intermolecular modes that “disrupt” the double-hydrogen bonds couple strongly to the zero-order O–H and C–H-stretch modes, which also become strongly coupled as well.
Next, we consider spectra obtained with the MD method. These calculations have been reported^{33,35} and so details can be found in that literature. In brief, the MD simulations are standard constant-energy (NVE) ones, where the total energy corresponds to a temperature of 300 K, and additional ones in which harmonic ZPE is given to each of the normal modes. The latter follows the “semi-classically prepared molecular dynamics approach”,^{100} but we denoted it as quasi-classical molecular dynamics (QCMD) in a recent paper on small protonated water clusters.^{101} The procedure in QCMD is quasi-classical, so it is subject to the well-known ZPE leak, which can cause artificial broadening of spectral bands. However, the hope is that this quasi-classical approach is able to capture the anharmonic effects much better than the NVE one. The spectra are obtained by (unfiltered) Fourier transformation of the dipole–dipole autocorrelation functions, averaged over roughly 100 trajectories.
Fig. 7 shows the comparison of spectra calculated by VSCF/VCI, MD NVE at 300 K, and QCMD for FAD and two isotopologs. As seen, at the level of graphical resolution, there is good agreement among these methods, with the exception of the location of the complex bands. The MD approaches do capture the complexity of VCI bands; however, their location is significantly up-shifted from the VCI bands, as well as the corresponding experimental ones. So these MD approaches do not capture the large VCI down-shift of these bands relative to the harmonic results; in fact the MD results are very close to harmonic frequencies. Also, a word of caution is needed about the MD complex bands. In the QCMD calculation of the IR spectrum of HONO the O–H stretch appears complex in contrast to the quantum band, which is sharp.^{100} Therefore, the MD approach, while generally very good for lower-frequency bands (numerical comparisons with experiment are given in Table 1), is not accurate for the high-frequency ones. This is not unexpected, since the energy content at 300 K is much less than the quantum energy associated with the fundamental transition. The QCMD method is, unfortunately, not much better than the MD one probably due to ZPE leak from the O–H stretch, initially prepared with its ZPE.
Fig. 7 Comparison of the spectra obtained from VSCF/VCI, classical NVE MD at 300 K, and QCMD, for FAD and its deuterated isotopologs. Reproduced from ref. 35 with permission from the Royal Society of Chemistry. |
Mode | Harmonic | VSCF | VCI | AIMD^{b} | Exp. |
---|---|---|---|---|---|
a Very strong mixing; the band spans over 200 cm^{−1} range. b NVE trajectories at 300 K. c Jet-cooled FTIR spectroscopy, ref. 72. d Raman spectroscopy, ref. 102 and 103. e Raman spectroscopy, ref. 73. f FTIR spectroscopy, ref. 79. | |||||
1 | 70 | 103 | 96 | 69 | 69 ^{ } |
2 | 167 | 171 | 178 | 151 | 137^{d} |
3 | 170 | 250 | 209 | 168 | 168 ^{ } |
4 | 209 | 204 | 213 | 202 | 194^{e} |
5 | 254 | 303 | 273 | 248 | 230^{d} |
6 | 275 | 277 | 286 | 266 | 268, ^{ } 264 ^{ } |
7 | 693 | 688 | 700 | 686 | 677^{d} |
8 | 716 | 713 | 726 | 708 | 698 ^{ } |
9 | 956 | 981 | 945 | 933 | — |
10 | 970 | 1007 | 975 | 950 | 922 ^{ } |
11 | 1084 | 1095 | 1079 | 1075 | 1060^{d} |
12 | 1100 | 1112 | 1091 | 1089 | 1050 ^{ } |
13 | 1255 | 1246 | 1252 | 1236 | 1214^{d} |
14 | 1258 | 1248 | 1252 | 1243 | 1218 ^{ } |
15 | 1406 | 1381 | 1388 | 1387 | 1364 ^{ } |
16 | 1408 | 1384 | 1391 | 1397 | 1375^{d} |
17 | 1448 | 1421 | 1424 | 1428 | 1454 ^{ } |
18 | 1481 | 1458 | 1457 | 1465 | 1415^{d} |
19 | 1715 | 1681 | 1671 | 1660 | 1670^{d} |
20 | 1780 | 1751 | 1761 | 1768 | 1746 ^{ } |
21 | 3095 | 2969 | 2990 | 3050–3100 | 2949^{d} |
22 | 3097 | 2942 | ∼2960 | 3050–3100 | 2939 ^{ } |
23 | 3232 | 2767 | —^{a} | 3200–3400 | — |
24 | 3326 | 2953 | — ^{ } | 3200–3400 | 3084 ^{ } |
As an important aside, we note that we were also motivated to perform these MD calculations by recent MD simulations of the IR spectra of (HCOOH)_{2} and (DCOOH)_{2}, using model potential and dipole moment surfaces, reported along with the 300 K experimental spectra.^{80} These authors adjusted model PESs, based on some ab initio data, mainly by changing the double-proton transfer barrier height, to produce good agreement with experiment for the complex O–H stretch band. The potential barrier was reported as being between 5 and 7 kcal mol^{−1}, which is significantly lower than the calculated CCSD(T) one^{97} and also the one from our PES.^{31} It appears that the effect of lowering the barrier is to also lower the harmonic frequency of the O–H stretch, which directly results in a lowering of the position of the MD band. However, the results above question the assumption of the accuracy of the MD approach to determine the barrier height.
Finally, we show in Table 1 the 24 calculated fundamentals of FAD from the indicated method, as well as experiment using Raman and IR spectroscopy. The ones labeled VSCF are new calculations. These are actually the virtual-state energies obtained from the ground-vibrational-state VSCF Hamiltonian. The AIMD calculations are the NVE ones at an energy that corresponds to 300 K, and since half of the modes are IR inactive, the correlation function of the coordinates and the power spectrum were computed and used to obtain the frequencies. The numbers in boldface indicate IR active modes (see ESI of ref. 33). Modes 22 and 24 are the IR active C–H and O–H stretch modes. Overall, this table shows that VSCF and VCI are able to capture the down-shift in the high-frequency intramolecular modes, especially the C–H and O–H stretches, and for these modes, VSCF are generally in good agreement with VCI results. However, while VSCF do capture mode anharmonicity and mean-field mode coupling, they cannot describe the strong mode mixing, for which the VCI approach is needed.
Starting with the highest energy mode, mode 24, and working down to mode 12, we see the expected down-shift in the VCI energies relative to the harmonic ones. This positive anharmonicity is expected for high-frequency stretch modes. The magnitude of the down-shift decreases as the mode energies decrease. Note, the VSCF results are generally in accord with the trends seen in the VCI energies, and in general, the VSCF energies are much closer to the VCI ones than to the harmonic energies. The AIMD energies also follow this trend, although the down-shifts are considerably smaller than the VSCF ones for the high-frequency C–H and O–H stretches. For the modes below roughly 1000 cm^{−1} the trends relative to the harmonic energies are not uniform, i.e., there are both down-shifts and up-shifts. Further, the VCI energies are generally up-shifted from the HO ones as are the VSCF ones, although not as generally so. The MD energies, by contrast are uniformly down-shifted from the HO ones, and generally in somewhat better agreement with experiment for these modes than the VCI energies. Note, the average absolute difference in the VCI and MD energies for these modes is only 22 cm^{−1}, which is roughly twice the average absolute difference in the AIMD and experimental energies. Very likely, the slight up-shifts in the VSCF and VCI frequencies of the low-frequency modes (below 1000 cm^{−1}) are due to the use of rectilinear normal coordinates for these intermolecular modes. For these low-frequency modes, the AIMD results are probably more realistic, provided resonances are not a concern and the molecule has only one minimum, or several minima with sufficiently high barrier. Down-shifting of higher-frequency intramolecular modes, on the other hand, cannot be captured using harmonic or AIMD approaches. The VCI results are generally much more accurate for these modes than the MD results. Overall, the agreement with experiment is good, although less than the current state-of-the-art value of around 5 cm^{−1} noted in the Introduction. This is likely due to the some deficiencies in the fitted PES and also the convergence of the VSCF/VCI energies being in the range of 10 or so cm^{−1}.
The brief review of the VSCF/VCI quantum treatment of the vibrational dynamics of 10-atom FAD is a fair indication of the state-of-the-art quantum treatments of molecules in the gas phase. We also direct the reader to recent reviews of rigorous and general methods for more general applications to smaller molecules.^{9,10}
For molecules with tens of atoms, rigorous quantum approaches are years away. However, for a large class of non-covalent molecular clusters, approximate coupled-mode quantum approaches are feasible and have been developed.
One general approach to describe the vibrational dynamics of a large molecule or cluster is to consider a subset of all the vibrational modes. This is a straightforward approach that is used in many field of computational chemistry, but where typically the harmonic approximation is made to describe the modes. We have used this approach as well, but treating mode anharmonicity and mode coupling. This is a standard option in MULTIMODE and some examples of doing this in our group, using VSCF/VCI, are H_{7}^{+},^{104} H_{7}O_{3}^{+} and H_{9}O_{4}^{+}.^{105–107} With this straightforward approach, convergence of vibrational bands of interest can be investigated by systematically adding more modes to the calculation.
Another approach is to consider vibrational motion of a subset of atoms while holding the remaining atoms fixed. An example of this is the vibrations of a molecular adsorbate on a solid surface with the positions of the atoms of the solid held fixed. Another is the vibrations of CH_{4} in a clathrate hydrate, where the surrounding water molecules are held fixed. In fact this example was considered, but at the harmonic level and with model potentials, for the CH_{4} and the 2-body CH_{4}–H_{2}O interaction.^{108,109} Earlier important work along these lines, including a treatment of anharmonicity and mode coupling was reported for small water clusters.^{110} We have termed this approach, with mode coupling, the local monomer method. We describe it below and then present some recent results using it for the vibrations of methane in clathrate hydrates.^{45}
Our work in this area was motivated by water clusters. We developed and tested the local monomer model and used it to obtain IR spectra of waters clusters and ice,^{36–38,40} employing the WHBB suite of potential and dipole moment surfaces.^{113} We briefly review this local monomer model next.
To begin, we start with the full mass-scaled force constant (Hessian) for a large molecule or cluster of N atoms, organized in what we might term a monomer-centric fashion. For example and for simplicity, consider the case of M identical monomers, for example a water cluster, in which case F is given by
(9) |
F can be written as a block-diagonal matrix F_{DIAG} plus an off-diagonal one F_{OFFD} as follows:
(10) |
(11) |
It is important to note that eigenvectors, denoted C_{DIAG} are those of F_{DIAG} and not just the diagonal blocks of F_{DIAG}. Thus, these eigenvectors are full dimensional, but with only non-zero elements corresponding to the blocks in F_{DIAG}. The eigenvector matrix is of the form
(12) |
The local normal mode approach is in a sense, a divide-and-conquer method, which deconstructs a full 3N × 3N Hessian into smaller independent blocks, each of which can be calculated and diagonalized independently. This is a major advantage of the method.
As reviewed below, these local normal modes can be used in a coupled-mode treatment of the monomers in a cluster. However, to describe monomer-monomer coupling requires a method to couple these local monomer modes. We did suggest a promising, albeit somewhat heuristic, Hückel approach in which the coupling elements are obtained by comparing the harmonic frequencies of a local and full normal mode analysis. More details for the interested reader are in the literature.^{37,114}
Another recent strategy to develop localized modes that is, however, based on a full normal mode analysis is to localize full normal modes.^{27,115–124} In general, localization of a subset of k normal modes of similar harmonic frequencies is done. There are several criteria for mode localization, largely borrowed from methods to localize orbitals in electronic structure theory, and the interested reader is referred to the literature for the details. The approaches amount to a rotation of the normal mode eigenvectors to a new set that satisfy a localization condition, i.e., reducing the extent of many-atom displacements. An important point about mode localization is that the transformation of normal-mode eigenvectors produces a non-diagonal mass-scaled force constant matrix. Of course, the original transformation from say mass-scaled Cartesian coordinates to normal modes is made to eliminate off-diagonal bilinear coupling. The benefit, in exchange for re-introducing bilinear coupling, is to localize the coupling to a small group of atoms. In many cases, e.g., CH-stretches, this does appear to be the correct description. Obviously, there are other cases where this localization may be detrimental.
This localized normal mode approach has been used in VSCF and VSCF/VCI calculations of vibrational energies of water clusters and polyatomic molecules, where it has been demonstrated that near-converged results can be achieved with smaller n-mode couplings of the potential than using standard normal modes. However, it should be noted that these approaches generally require a full-normal mode analysis, which itself can be very costly if done using a high-level electronic structure theory. The FALCON approach^{27,118,119} aims to avoid this, using a method somewhat like the local monomer one, but for large molecules instead of non-covalent clusters.
[_{m} + V_{m}(Q_{m}) − E^{(m)}_{i}]ϕ^{(m)}_{i}(Q_{m}) = 0 | (13) |
Finally, we note that we have used local normal modes in VPT2 theory, which we termed local VPT2, and showed the effectiveness of this approach for the water dimer.^{128}
Raman spectroscopy has been widely applied to study the methane hydrates.^{130–136} The magnitudes of the down-shift in the frequency of the methane symmetric stretch depend on the types of cages. Therefore, Raman spectroscopy has been used to identify the hydrate structures and the cage occupancies. The IR spectra of methane in 5^{12} and 5^{12}6^{2} cages have also been reported,^{137,138} and similar down-shift of the frequencies of the asymmetric stretch modes were observed.
Theoretically, the intramolecular vibrations of enclathrated methane have been calculated at the harmonic level using DFT,^{139–141} or using ab initio molecular dynamics simulations.^{142–146} These studies are able to reproduce the experimental trend of the frequency shift but not with quantitative accuracy.
To achieve at least semi-quantitative agreement with experiments, an accurate PES and rigorous (quantum) vibrational dynamics approach are required. We reported such a PES to describe the methane clathrates based on the many-body expansion:^{147,148}
(14) |
In that work we focused on the intramolecular vibrations of methane, and in this case, the local monomer model is expected to work well, because the interaction between methane and the cage is weak and so the off-diagonal force constant matrix F_{OFFD} is expected to be small. This can be verified by comparing the result of local monomer normal-mode analysis for CH_{4}(H_{2}O)_{20} with the full-dimensional normal-mode analysis. Fig. 8 shows the approximate density of states from the two normal-mode analyses. As expected, and as is seen in panel (a) of the figure, the local monomer model is not able to describe the low-frequency collective motions very well. However, for the high-frequency intramolecular vibrations, the agreement between local monomer and full-dimensional normal-mode analysis is excellent, especially for the vibrations of methane.
Fig. 8 Density of states calculated at harmonic level using local monomer and full-dimensional normal mode analysis in: (a) low frequency libration and (b) intramolecular vibration regimes. Adapted with permission from ref. 45. Copyright 2016 American Chemical Society. |
To go beyond the harmonic approximation, the local monomer Schroedinger equation of the embedded methane monomer, cf., eqn (13), was solved. In this case, Q_{m} is the set of nine intramolecular local normal modes of methane, obtained from a local normal mode analysis. The VSCF/VCI approach was used to solve this Schroedinger equation, using the code MULTIMODE with a 4MR of the potential, and the details of the calculations can be found in ref. 45.
Table 2 shows the energy of the symmetric stretch fundamental of methane in four cages, from VSCF/VCI calculations, experiment,^{135} and Car–Parrinello molecular dynamics (CP-MD).^{146} The calculated value of the fundamental frequency of this mode in the gas phase is 2917 cm^{−1}, in excellent agreement with experiment.^{150} The VSCF/VCI calculations capture the down-shift of the frequency when methane is confined in the cage, and also predict a larger down-shift in large (5^{12}6^{2} and 5^{12}6^{4}) cages. The overall agreement with experiment is good; however, our calculations systematically underestimated the frequencies by about 10 cm^{−1}. This small difference between the theory and experiments may be due to the inaccuracy in the PES (the truncation in the many-body expansion, the fitting error, and the inaccuracy in the electronic structure theory) and the approximations made in the local monomer model. Nevertheless, the VSCF/VCI calculations, using local monomer model, achieved an unprecedented level of accuracy to predict the vibrational frequencies of methane inside clathrate hydrate cages.
The frequencies of asymmetric stretches of methane in 5^{12} and 5^{12}6^{2} cages are presented in Table 3, together with Fourier transform infrared spectrum (FTIR)^{137,138} and recent CP-MD results.^{146} The frequency of this triply-degenerate mode in the gas-phase, calculated with MULTIMODE, is 3021 cm^{−1}, which is only 2 cm^{−1} higher than the experiment.^{150} As expected, this degeneracy is lifted in the hydrate cages and the value of the splittings is a sensitive function of the interaction with the surrounding waters. The assignment of the experimental spectra was somewhat ambiguous, so only estimated frequencies of the band centers are listed. Similar to the symmetric stretch, the VSCF/VCI frequencies with local monomer modes are within 10 cm^{−1} compared to the experiments, and this level of accuracy is way beyond the other methods such as the CP-MD shown in the tables. In addition, in the VSCF/VCI calculation, we were able to obtained the splittings of the three degenerate stretches due to the perturbation of the cages.
Cage | MULTIMODE | FTIR^{a} | CP-MD^{b} |
---|---|---|---|
a Ref. 137 and 138. b Ref. 146. | |||
5^{12} | 3010, 3021, 3027 | ∼3016 | 2979 |
5^{12}6^{2} | 2994, 3007, 3014 | ∼3000 | 2963 |
Gas phase | 3021 |
In addition to the VSCF/VCI calculations, diffusion Monte Carlo (DMC)^{152,153} were done to obtain the ground state wavefunction of the confined methane in four different cages, with methane treated in full dimensionality while the water molecules are held fixed. These calculations accurately characterize the vibration–translation–rotation motion of the methane molecule in the ground vibrational state. The unbiased random walk algorithm was applied,^{55,154} and the calculations were done in Cartesian coordinates instead of the normal coordinates. The details of the DMC calculations are given in ref. 45.
It is worth noting that these simulations require ∼10^{9} potential evaluations per DMC trajectory, and each potential evaluation requires ∼20 methane–water two-body and 200–400 methane–water–water three-body potential evaluations (depending on the size of the cage). These calculations would not be possible without the analytical PES. Even with our efficient PES, a single calculation of the CH_{4}(H_{2}O)_{20} (the smallest system) still takes about four days of cpu time on a sixteen-processor computer.
The DMC ground-state wavefunctions of methane in four cages, represented as isosurfaces, are shown in Fig. 9. The isovalue is 40% of the maximum wavefunction amplitude in all the figures. The ball-and-stick models in the figures represent the geometries at the minimum of the potential, while the clouds represent the nuclear wavefunctions. Compared to the VCI calculations, which were restricted to the nine intra-molecular vibrations, the DMC nuclear wave functions provide additional information about the translations and hindered rotations of methane in the cage. It can be clearly seen that the wavefunctions in the large cages (5^{12}6^{2} and 5^{12}6^{4}) are more delocalized than those in the small cages (5^{12} and 4^{3}5^{6}6^{3}).
Fig. 9 Ground state nuclear wave function of methane in 5^{12}, 5^{12}6^{2}, 4^{3}5^{6}6^{3} and 5^{12}6^{4} clathrate cages. Reprinted with permission from ref. 45. Copyright 2016 American Chemical Society. |
The local monomer approach we have used employs snapshots of a standard MD simulation and performs quantum calculations on the collection of monomers in each snapshot and then averages the resulting spectra over the snapshots. This approach cannot describe time-dependent coherence effects; however, at higher temperature these are expected to be small. An interesting comparison of CMD, thermostat RPMD and local monomer calculations of the IR spectrum of liquid water at 298 K, using an empirically adjusted model water potential, shows that these three approaches are in good accord with each other.^{41} The local monomer approach was also applied successfully to intramolecular vibrations of water ice,^{39} and accurate simulations of the IR spectrum of room temperature water and of heavy (deuterated) water.^{42,127} The latter calculations used an ab initio, flexible dipole moment surface,^{126} which is needed to obtain correct band intensities.
Clearly, more theoretical work needs to be done using quantum approaches for the condensed phase and we anticipate progress in this direction in the coming decade.
The IR spectra of the formic acid dimer and isotopologs were presented as an example to illustrate both simple and complex IR bands, seen experimentally, and thus a major challenge for theory to describe fully. Using recent, full-dimensional ab initio potential and dipole moment surfaces, VSCF/VCI calculations of the IR spectra show good agreement with experiment for all bands. The widely used AIMD approach to calculate IR spectra was also examined using the same potential and dipole moment surfaces. The tests found very good accuracy for some modes but significant quantitative inaccuracy for the complex O–H band. However, qualitatively, the fractionation of the band is reproduced by the approach. In this example, the strong coupling of the O–H and C–H stretches to the intermolecular modes that “disrupt” the double hydrogen bond was found to be responsible for the complexity of the band.
The third part of the Perspective dealt with extending the VSCF/VCI method to molecular clusters and towards the condensed phase via the anharmonic local monomer approach. Following a review of the method, an application to the vibrational energies of methane clathrate hydrate was given. Using full-dimensional 2-body CH_{4}–H_{2}O and 3-body CH_{4}–H_{2}O–H_{2}O potentials, the vibrational energies of the caged CH_{4} were obtained and shown to be in good agreement with experiment.
Final comments about recent research on localized modes and entering the condensed phase with quantum methods were made. These are areas of research that we anticipate great strides to be made in the future.
This journal is © the Owner Societies 2019 |