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

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

Received 5th August 2018 , Accepted 24th October 2018

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.

1 Introduction

This Perspective describes post-harmonic quantum approaches to vibrational dynamics, with an emphasis on infrared (IR) spectroscopy, of molecules and also molecular clusters. This is a very active area of research, as attested to by recent themed issues of Phys. Chem. Chem. Phys.1 and J. Phys. Chem. A.2 In this Perspective, the focus is on the challenges to theory to provide accurate quantum simulations of these vibrational dynamics as well as useful interpretations of the simulations. Although not part of this stated focus, we also remark on the widely-used ab initio molecular dynamics approach (“AIMD”) to the calculation of IR spectra and present recent calculations from our group on the IR spectra the formic acid dimer and isotopologs using this approach.

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/QFF11–13 (vibrational second order perturbation theory with quartic force fields)

• VSCF5,14–17 (vibrational self-consistent field theory)

• VCI and VSCF/VCI4,5,7–10,18–25 (vibrational/virtual state configuration interaction theory)

• VSCF/PT25,26 (VSCF + second-order perturbation theory)

• VSCF/VCC6,27 (VSCF + vibrational coupled-cluster theory)

• MCTDH28–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 CH4 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.

2 Quantum approaches that build on the separable harmonic model

To begin, let H be the exact Hamiltonian, given simply as the sum of kinetic and potential energy operators, T + V. Zero-order eigenfunctions are eigenfunctions of a zero-order Hamiltonian, given by
H(0)ΦK = E(0)KΦK, H(0) = T + V(0).(1)
For example, if V(0) is a separable harmonic-oscillator potential in say normal modes, the zero-order eigenfunctions are harmonic-oscillator functions. These zero-order eigenfunctions can be used to represent the exact eigenfunctions of H by
image file: c8cp04990d-t1.tif(2)
The goal is to obtain the expansion coefficients c(L)K. Once these are known, the problem is “solved”; however, the challenge remains to extract information from these coefficients that can be used to interpret the answer. Obviously, if a single coefficient dominates, e.g., has a magnitude of say 0.9 or greater, the interpretation is essentially the one from the familiar normal mode analysis and the experimental band should be “simple” to interpret. The challenge comes from those interesting cases where there is not one or even several dominant coefficients. In this case, one could conclude that the harmonic normal mode analysis has broken down and that there is strong mode–mode coupling. This is usually manifested experimentally by a broad and complex band, which is perforce difficult to interpret. The goal in this case theoretically is to at least determine which modes are strongly coupled and ultimately to understand the source of the coupling.

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−6i=1ϕki(Qi). While this is high-dimensional, it is basically trivial because it is written in a simple product form. Note K is a row vector [k1, k2,…k3N−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 103N−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.

2.1 VPT2

We begin with second-order perturbation theory, where the perturbation V′ is given by VV(0). The wavefunction including the first-order correction to it (which gives the second-order correction to the energy) is given by
image file: c8cp04990d-t2.tif(3)
We point out that this is still of the form of an expansion in the zero-order basis, but the expansion coefficients are given by approximate analytical expressions. The usual assumption is that the zero-order term in the expansion dominates (and is given the coefficient 1.0) and the other terms are (hopefully) small corrections. There are two issues with this expression. One is computational, i.e., the evaluation of the potential matrix element in the high-dimensional (3N − 6) space of the wavefunction. The second is small denominators, known as resonances (Fermi and Darlington–Dennison resonances being well-known examples).3 These resonances can result in large corrections which invalidate the assumption that these coefficients are small. This breakdown of VPT2 is actually of interest, as it may signal something unusual in the dynamics. Degenerate perturbation theory, i.e., basically a small VCI, is able to handle these resonances.48 This approach, however, still has a limited range of validity. The other issue mentioned, i.e., the computational cost of matrix elements of the perturbation potential is actually not a problem if one follows the conventional approach and expands V′ in a quartic force field (QFF), using perturbation parameters in the formal theory for the third and fourth order terms.11 In fact, the QFF is a Taylor series expansion of the full potential in normal coordinates up to fourth order. The form of the QFF is, using modern jargon, a sum-of-products (SOP). This makes the evaluation of matrix elements just a sum of products of 1d elements, which in a harmonic basis can be done analytically.

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.

image file: c8cp04990d-f1.tif
Fig. 1 Plot of the Morse potential and quartic expansion around the minimum.

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.

2.2 VSCF and VSCF + VCI

Returning to eqn (2), the rigorous approach is to determine the exact expansion coefficients, where we have used the harmonic-oscillator wavefunctions as the example basis. Another, more sophisticated basis consists of the virtual states of a VSCF Hamiltonian, analogous to the basis used in electronic structure theory, and defined below is the one we use. In fact the VSCF method itself is a general post-harmonic method that we briefly review next.

The VSCF approach dates from 197814–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(Qi), where the ϕni(Qi) 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 = T1 + T2 + T3 + V(Q1,Q2,Q3) these are

image file: c8cp04990d-t3.tif(4)
These are coupled, integro-differential equations for the optimum modal wavefunctions. They are solved in the usual iterative fashion to self-consistency. At convergence these are just three 1d Schroedinger equations for each modal with effective potentials that depend on the quantum state of the other modals. These are, in principle, straightforward to solve; however, for an n-mode molecule, the effective potential for each modal requires the evaluation of an (n − 1) dimensional matrix element. (For this pedagogical example, we assume that the kinetic energy operator of each modal is independent of the other modes.) Clearly, the modals depend on the quantum state of the other modals.

For the ground vibrational state, ni = 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:

image file: c8cp04990d-t4.tif(5)
where NA is the Avogadro's number, h is the Planck constant, c is the speed of light, ε0 is the vacuum permittivity, νfi is the wavenumber of the transition, and μα(Q) is the dipole moment surface of the molecule. The factor (ninf) is the difference in the fraction of molecules in the initial and final states, and when the temperature is low, for the transitions originating from the ground vibrational state this term tends to 1. Converted to conventional units (km mol−1) (assuming ninf = 1), this formula becomes
image file: c8cp04990d-t5.tif(6)
with ν in cm−1 and dipole in debyes. This expression is for zero angular momentum and so clearly it is not appropriate for line-list ro-vibrational spectroscopy. MULTIMODE can perform ro-vibrational line-list spectroscopy with the appropriate expression for the intensity; an example is C2H6.54

3 The Watson Hamiltonian and the n-mode representation of potential

Quantum approaches to vibrational dynamics are generally performed using a set of 3N − 6 internal vibrational coordinates, as these can be done for specified values of the total angular momentum. (An important exception is quantum diffusion Monte Carlo calculations, which are routinely done in all the Cartesian coordinates.55) The Hamiltonian operator is transformed to a matrix using a basis representation of the wavefunction, cf., eqn (2) for example. The details depend on the choice of internal coordinates,10 but in a basis of zero-order functions, e.g., harmonic oscillator or virtual states, the potential matrix generally requires multidimensional quadrature. This is one, if not the major, computational bottleneck in applying VCI approaches, for two reasons. First, a strictly numerical quadrature scales exponentially with dimensionality. For illustration, if there are 15 quadrature points per degree of freedom, each potential matrix element effort scales as [scr O, script letter O](153N−6). Second, the size of the H-matrix also scales exponentially with the number of degrees of freedom. For example, if there are 10 basis functions per degree of freedom, the H-matrix is 103N−6 × 103N−6. So, this is basically hopeless even for a nonlinear tetraatomic molecule with 6 degrees of freedom. Clearly, efficient, but computationally rigorous strategies are needed to deal with these two aspects of exponential scaling. Some of these strategies can be found in recent reviews.4–10

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

image file: c8cp04990d-t6.tif(7)
where α(β) represent the x, y, z coordinates, Ĵα and [small pi, Greek, circumflex]α are the components of the total and vibrational angular momenta respectively, μαβ is the inverse of effective moment of inertia, and V(Q) is the potential. The number of normal modes is denoted by F, and for non-linear molecules F equals 3N − 6.

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 H3O+. 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 POLYMODE59 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 by52

image file: c8cp04990d-t7.tif(8)
where V(1)i(Qi) is the one-mode potential, i.e., the 1D cut through the full-dimensional PES in each mode, one-by-one, V(2)ij(Qi,Qj) is the intrinsic 2-mode potential among all pairs of modes, etc. Here, intrinsic means that the any n-mode term is zero if any of the arguments is zero. Also, each term in the representation is in principle of infinite order in the sense of a Taylor series expansion. So, returning to Fig. 1, V(1)(Q) is the full Morse potential.

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[thin space (1/6-em)]:[thin space (1/6-em)]1 Fermi resonances, it cannot describe 1 + 1[thin space (1/6-em)]:[thin space (1/6-em)]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 H3O+ 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 MULTIMODE68 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 reviews4–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

4 IR spectra of the formic acid dimer and isotopologs

The formic acid dimer (FAD) provides an excellent case study where theory meets experiment. This is a major challenge for theory, if all (24) vibrational modes are considered. The IR and Raman spectra of FAD have been reported and analyzed for many years.69–81 The dimer is bound symmetrically by two hydrogen bonds that strongly perturb some vibrational bands in the IR and Raman spectra. These bonds can re-arrange symmetrically by a double-proton transfer and this results in a small tunneling splitting, which has been determined by high-resolution IR spectroscopy.76,81

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 13[thin space (1/6-em)]475 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.

image file: c8cp04990d-f2.tif
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.

image file: c8cp04990d-f3.tif
Fig. 3 Double harmonic spectra of formic acid dimer and indicated isotopologs. Reproduced from ref. 35 with permission from the Royal Society of Chemistry.

image file: c8cp04990d-f4.tif
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.

image file: c8cp04990d-f5.tif
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.

image file: c8cp04990d-f6.tif
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 reported33,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.

image file: c8cp04990d-f7.tif
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.
Table 1 Comparison of the harmonic, VSCF, VCI, AIMD, and experimental frequencies (in cm−1) of the fundamental transitions of the formic acid dimer
Mode Harmonic VSCF VCI AIMDb 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 137d
3 170 250 209 168 168
4 209 204 213 202 194e
5 254 303 273 248 230d
6 275 277 286 266 268, 264
7 693 688 700 686 677d
8 716 713 726 708 698
9 956 981 945 933
10 970 1007 975 950 922
11 1084 1095 1079 1075 1060d
12 1100 1112 1091 1089 1050
13 1255 1246 1252 1236 1214d
14 1258 1248 1252 1243 1218
15 1406 1381 1388 1387 1364
16 1408 1384 1391 1397 1375d
17 1448 1421 1424 1428 1454
18 1481 1458 1457 1465 1415d
19 1715 1681 1671 1660 1670d
20 1780 1751 1761 1768 1746
21 3095 2969 2990 3050–3100 2949d
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) one97 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 H7+,104 H7O3+ and H9O4+.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 CH4 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 CH4 and the 2-body CH4–H2O 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

5 Anharmonic local monomer approach for clusters and towards the condensed phase

5.1 Local normal modes

The local monomer model begins with local normal modes. For an adsorbate the idea of treating just the modes of the adsorbate is easily motivated. For a molecular cluster consisting of many monomers, treating each monomer's vibrational modes independently seems problematic, at first glance. Similarly, treating the vibrational modes of a subset of atoms in a molecule might also seem problematic. Nevertheless both approaches have been used in the past, the latter by the molecular tailoring method111 and the former used for small water clusters110,112 and the CH4 clathrate hydrate already noted.109

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

image file: c8cp04990d-t8.tif(9)
such that the mass-scaled Cartesian displacements are arranged so that the diagonal blocks of F, i.e., F11, F22,…, FMM, are force constant matrices for each monomer. For a cluster of water monomers in this example, each block Fij is of size 9 × 9.

F can be written as a block-diagonal matrix FDIAG plus an off-diagonal one FOFFD as follows:

image file: c8cp04990d-t9.tif(10)
image file: c8cp04990d-t10.tif(11)
By discarding the off-diagonal blocks, the full mass-scaled force constant matrix is approximated by the block diagonal matrix and the normal mode analysis proceeds simply by diagonalizing each block. The eigenvectors are the local monomer normal modes and the eigenvalues are the corresponding square of the local monomer frequencies.

It is important to note that eigenvectors, denoted CDIAG are those of FDIAG and not just the diagonal blocks of FDIAG. Thus, these eigenvectors are full dimensional, but with only non-zero elements corresponding to the blocks in FDIAG. The eigenvector matrix is of the form

image file: c8cp04990d-t11.tif(12)
and we note that these eigenvectors form an orthonormal basis that spans the space of all local normal modes.

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 approach27,118,119 aims to avoid this, using a method somewhat like the local monomer one, but for large molecules instead of non-covalent clusters.

5.2 Anharmonic local monomer model

We used local normal modes extensively in numerous applications to water clusters36–40,42,125 with VSCF/VCI calculations as implemented in the code MULTIMODE.19 For each monomer the Schroedinger equation
[[T with combining circumflex]m + Vm(Qm) − E(m)i]ϕ(m)i(Qm) = 0(13)
is solved, where Qm denotes the set of local normal modes, and Vm(Qm) is the full potential of the monomer perturbed by the other monomers. For water clusters, solid water and liquid water each monomer has 9 local normal modes and in most of our calculations we have focused on the three intramolecular ones. We have included intermolecular modes in calculations of the IR spectra of water126 and heavy water.127 We illustrate this approach for local monomer calculations of CH4 clathrate hydrates.45

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

5.3 Application to methane clathrate hydrates

Clathrate hydrates are inclusion compounds where small gas molecules are trapped inside cages that are formed by hydrogen-bonded water molecules. Some typical cages in clathrates are pentagonal dodecahedral ones with 12 pentagonal faces (denoted as 512), tetrakaidecahedral (51262), hexakaidecahedral (51264), and irregular dodecahedral ones (435663). In simple methane hydrates, the methane molecules occupy both 512 and 51262 cages,129 and they can occupy other cages in binary clathrate hydrates (with another guest gas).

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 512 and 51262 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

image file: c8cp04990d-t12.tif(14)
where “M” represents methane and “W” represents water, V(1), V(2) and V(3) are one-body, intrinsic two-body, and intrinsic three-body energies, respectively. As noted in an earlier theoretical calculation of the binding energy of CH4(H2O)20,149 the contributions from higher-order terms are less than 10% of the binding energy, and thus these terms are neglected in our many-body potential. The monomer potentials of methane150 and water151 used in this expansion are spectroscopically accurate ones. The water two-body and three-body are included in the WHBB water potential.113 The intrinsic methane–water and methane–water–water potentials were reported recently.147,148 The full-dimensional two-body interaction is a fit to 36[thin space (1/6-em)]930 CCSD(T)-F12b/haTZ (aVTZ for C, O and VTZ for H) electronic energies, and the three-body is a fit to 23[thin space (1/6-em)]132 MP2-F12/haTZ energies. The details of the potential can be found in our publications.45,147,148

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 FOFFD is expected to be small. This can be verified by comparing the result of local monomer normal-mode analysis for CH4(H2O)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.

image file: c8cp04990d-f8.tif
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, Qm 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 (51262 and 51264) 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.

Table 2 Anharmonic fundamentals (cm−1) of methane symmetric stretch in four cages, calculated by MULTIMODE, and comparison with Raman spectra and CP-MD simulation. The gas-phase anharmonic frequency of CH4 calculated with MULTIMODE is also listed. Adapted with permission from ref. 45. Copyright 2016 American Chemical Society
a Ref. 135. b Ref. 146.
512 2902 2914 2876, 2879
51262 2893 2901 2865
435663 2901 2910 2875
51264 2888 2902
Gas phase 2917

The frequencies of asymmetric stretches of methane in 512 and 51262 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.

Table 3 Anharmonic fundamentals (cm−1) of asymmetric stretches of methane in four cages, calculated by MULTIMODE, and comparison with FTIR and CP-MD simulation. The anharmonic fundamental frequency of the gas-phase methane is also listed. Adapted with permission from ref. 45. Copyright 2016 American Chemical Society
a Ref. 137 and 138. b Ref. 146.
512 3010, 3021, 3027 ∼3016 2979
51262 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 ∼109 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 CH4(H2O)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 (51262 and 51264) are more delocalized than those in the small cages (512 and 435663).

image file: c8cp04990d-f9.tif
Fig. 9 Ground state nuclear wave function of methane in 512, 51262, 435663 and 51264 clathrate cages. Reprinted with permission from ref. 45. Copyright 2016 American Chemical Society.

5.4 Remarks about the condensed phase

Finally, we “dip our toes” into extending quantum approaches to the condensed phase. This is truly a taunting challenge. At present the AIMD approach is the workhorse for such simulations.155 However, semi-classical,156,157 semi-quantum ring-polymer molecular dynamics (RPMD),158,159 centroid molecular dynamics (CMD),160–162 and the local monomer39,41 methods have also been applied to the condensed phase, notably to the IR spectrum of liquid water. The CMD and RPMD methods do capture the significant effects of anharmonicity of the O–H stretch in water,39,41 and in this respect are more accurate than pure classical MD simulations of the IR spectrum. However, these alternatives to the MD approach are significantly more computationally intensive except at high temperature where they merge with MD simulations. Also, since these methods are somewhat heuristic, it is not clear how they can be made systematically more accurate.

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.

6 Summary and conclusions

In this Perspective, we began with a pedagogical review of the most widely used quantum approaches to describe anharmonicity and coupling in molecules. The vibrational configuration interaction approach was introduced as the most general and rigorous method. Since it is an expansion of the exact vibrational wavefunction in term of a (zero-order) basis, which for pedagogical purposes we took as the familiar harmonic basis from a normal mode analysis, the information content is contained in the expansion coefficients. If a single or a small number of coefficients dominate the expansion, then the interpretation is the familiar one, i.e., the wavefunction is a harmonic normal-mode function in the former case or a sum of several resonantly (Fermi, Darling–Dennison) coupled such functions in the latter case. However, if there is not a dominant set of coefficients, this presents a challenge for interpretation. This situation usually signifies a complex broad experimental band (assuming the calculation is trustworthy) which generally calls for in depth examination of the strong mode coupling.

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 CH4–H2O and 3-body CH4–H2O–H2O potentials, the vibrational energies of the caged CH4 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.

Conflicts of interest

There are no conflicts to declare.


We thank the reviewers for constructive comments. Financial support is from NASA grant number NNX16AF09G.


  1. M. Hochlaf, Phys. Chem. Chem. Phys., 2013, 15, 9967–9969 RSC .
  2. M. L. Senent, M. Hochlaf and M. Carvajal, J. Phys. Chem. A, 2016, 120, 475–476 CrossRef CAS .
  3. E. B. Wilson, J. C. Decius and P. C. Cross, Molecular vibrations: the theory of infrared and Raman vibrational spectra, Dover Publ., New York, NY, 1980, pp. 197–200 Search PubMed .
  4. J. M. Bowman, T. Carrington and H.-D. Meyer, Mol. Phys., 2008, 106, 2145–2182 CrossRef CAS .
  5. T. K. Roy and R. B. Gerber, Phys. Chem. Chem. Phys., 2013, 15, 9468–9492 RSC .
  6. O. Christiansen, Phys. Chem. Chem. Phys., 2012, 14, 6672–6687 RSC .
  7. D. Oschetzki and G. Rauhut, Phys. Chem. Chem. Phys., 2014, 16, 16426–16435 RSC .
  8. A. G. Csaszar, C. Fabri, T. Szidarovszky, E. Matyus, T. Furtenbacher and G. Czako, Phys. Chem. Chem. Phys., 2012, 14, 1085–1106 RSC .
  9. J. Tennyson, J. Chem. Phys., 2016, 145, 120901 CrossRef PubMed .
  10. T. Carrington, J. Chem. Phys., 2017, 146, 120902 CrossRef .
  11. H. H. Nielsen, Rev. Mod. Phys., 1951, 23, 90–136 CrossRef CAS .
  12. J. F. Gaw, A. Willetts, W. H. Green and N. C. Handy, Adv. Mol. Vib. Collision Dyn., 1991, 169 CAS .
  13. J. Bloino and V. Barone, J. Chem. Phys., 2012, 136, 124108 CrossRef .
  14. J. M. Bowman, J. Chem. Phys., 1978, 68, 608 CrossRef CAS .
  15. G. D. Carney, L. L. Sprandel and C. W. Kern, Adv. Chem. Phys., 1978, 37, 305 CAS .
  16. J. M. Bowman, Acc. Chem. Res., 1986, 19, 202–208 CrossRef CAS .
  17. R. B. Gerber and M. A. Ratner, Adv. Chem. Phys., 1998, 97–132 Search PubMed .
  18. K. Christoffel and J. Bowman, Chem. Phys. Lett., 1982, 85, 220 CrossRef CAS .
  19. S. Carter, J. M. Bowman and N. C. Handy, Theor. Chem. Acc., 1998, 100, 191–198 Search PubMed .
  20. S. Carter and N. Handy, Comput. Phys. Commun., 1988, 51, 49–58 CrossRef CAS .
  21. E. Matyus, G. Czako and A. G. Csaszar, J. Chem. Phys., 2009, 130, 134112 CrossRef .
  22. S. N. Yurchenko, R. J. Barber, A. Yachmenev, W. Thiel, P. Jensen and J. Tennyson, J. Phys. Chem. A, 2009, 113, 11845–11855 CrossRef CAS .
  23. D. Lauvergnat and A. Nauts, Spectrochim. Acta, Part A, 2014, 119, 18–25 CrossRef CAS .
  24. T. Halverson and B. Poirier, J. Phys. Chem. A, 2015, 119, 12417–12433 CrossRef CAS .
  25. T. Petrenko and G. Rauhut, J. Chem. Theory Comput., 2017, 13, 5515–5527 CrossRef CAS .
  26. L. S. Norris, M. A. Ratner, A. E. Roitberg and R. B. Gerber, J. Chem. Phys., 1996, 105, 11261–11267 CrossRef CAS .
  27. E. L. Klinting, C. König and O. Christiansen, J. Phys. Chem. A, 2015, 119, 11007–11021 CrossRef CAS .
  28. H.-D. Meyer, U. Manthe and L. S. Cederbaum, Chem. Phys. Lett., 1990, 165, 73–78 CrossRef CAS .
  29. U. Manthe, H.-D. Meyer and L. S. Cederbaum, J. Chem. Phys., 1992, 97, 3199–3213 CrossRef CAS .
  30. G. A. Worth, M. H. Beck, A. Jäckle, O. Vendrell and H.-D. Meyer, The MCTDH Package, Version 8.2, (2000). H.-D. Meyer, Version 8.3 (2002), Version 8.4 (2007). O. Vendrell and H.-D. Meyer Version 8.5 (2013). Version 8.5 contains the ML-MCTDH algorithm. Current versions: 8.4.12 and 8.5.5 (2016). See http://mctdh.uni-hd.de/.
  31. C. Qu and J. M. Bowman, Phys. Chem. Chem. Phys., 2016, 18, 24835 RSC .
  32. C. Qu and J. M. Bowman, J. Chem. Phys., 2018, 148, 241713 CrossRef .
  33. C. Qu and J. M. Bowman, J. Phys. Chem. Lett., 2018, 9, 2604–2610 CrossRef CAS PubMed .
  34. P. Houston, B. L. Van Hoozen, C. Qu, Q. Yu and J. M. Bowman, Faraday Discuss., 2018 10.1039/c8fd00075a .
  35. C. Qu and J. M. Bowman, Faraday Discuss., 2018 10.1039/c8fd00077h .
  36. Y. Wang and J. M. Bowman, J. Chem. Phys., 2011, 134, 154510 CrossRef .
  37. Y. Wang and J. M. Bowman, J. Chem. Phys., 2012, 136, 144113 CrossRef PubMed .
  38. Y. Wang and J. M. Bowman, J. Phys. Chem. Lett., 2013, 4, 1104–1108 CrossRef CAS .
  39. H. Liu, Y. Wang and J. M. Bowman, J. Phys. Chem. Lett., 2012, 3, 3671–3676 CrossRef CAS .
  40. H. Liu, Y. Wang and J. M. Bowman, J. Phys. Chem. B, 2014, 118, 14124–14131 CrossRef CAS .
  41. M. Rossi, H. Liu, F. Paesani, J. Bowman and M. Ceriotti, J. Chem. Phys., 2014, 141, 181101 CrossRef .
  42. H. Liu, Y. Wang and J. M. Bowman, J. Chem. Phys., 2015, 142, 194502 CrossRef .
  43. J. M. Bowman, Y. Wang, H. Liu and J. S. Mancini, J. Phys. Chem. Lett., 2015, 6, 366–373 CrossRef CAS .
  44. Z. Homayoon, R. Conte, C. Qu and J. M. Bowman, J. Chem. Phys., 2015, 143, 084302 CrossRef .
  45. C. Qu and J. M. Bowman, J. Phys. Chem. C, 2016, 120, 3167–3175 CrossRef CAS .
  46. Q. Wang and J. M. Bowman, J. Chem. Phys., 2017, 147, 161714 CrossRef .
  47. T. Q. Bui, P. B. Changala, B. J. Bjork, Q. Yu, Y. Wang, J. F. Stanton, J. Bowman and J. Ye, Mol. Phys., 2018, 116, 3710–3717 CrossRef CAS .
  48. K. Yagi and H. Otaki, J. Chem. Phys., 2014, 140, 084113 CrossRef .
  49. D. A. Matthews, J. Vázquez and J. F. Stanton, Mol. Phys., 2007, 105, 2659–2666 CrossRef CAS .
  50. R. C. Fortenberry, Q. Yu, J. S. Mancini, J. M. Bowman, T. J. Lee, T. D. Crawford, W. F. Klemperer and J. S. Francisco, J. Chem. Phys., 2015, 143, 071102 CrossRef .
  51. Q. Yu, J. M. Bowman, R. C. Fortenberry, J. S. Mancini, T. J. Lee, T. D. Crawford, W. Klemperer and J. S. Francisco, J. Phys. Chem. A, 2015, 119, 11623–11631 CrossRef CAS .
  52. S. Carter, S. J. Culik and J. M. Bowman, J. Chem. Phys., 1997, 107, 10458–10469 CrossRef CAS .
  53. J. M. Bowman, S. Carter and X. Huang, Int. Rev. Phys. Chem., 2003, 22, 533 Search PubMed .
  54. S. Carter, A. R. Sharma and J. M. Bowman, J. Chem. Phys., 2012, 137, 154301 CrossRef .
  55. A. B. McCoy, Int. Rev. Phys. Chem., 2006, 25, 77–107 Search PubMed .
  56. J. K. G. Watson, Mol. Phys., 1968, 15, 479 CrossRef CAS .
  57. X. Huang, S. Carter and J. Bowman, J. Chem. Phys., 2003, 118, 5431–5441 CrossRef CAS .
  58. E. Kamarchik, Y. Wang and J. Bowman, J. Phys. Chem. A, 2009, 113, 7556–7562 CrossRef CAS PubMed .
  59. H. Romanowski, J. M. Bowman and L. B. Harding, J. Chem. Phys., 1985, 82, 4155–4165 CrossRef CAS .
  60. K. M. Dunn, J. E. Boggs and P. Pulay, J. Chem. Phys., 1986, 85, 5838–5846 CrossRef CAS .
  61. H. Rabitz and O. Aliş, J. Math. Chem., 1999, 25, 197–233 CrossRef CAS .
  62. L. Ostrowski, B. Ziegler and G. Rauhut, J. Chem. Phys., 2016, 145, 104103 CrossRef .
  63. B. Ziegler and G. Rauhut, J. Chem. Phys., 2016, 144, 114114 CrossRef .
  64. C. König and O. Christiansen, J. Chem. Phys., 2015, 142, 144115 CrossRef .
  65. J. O. Jung and R. B. Gerber, J. Chem. Phys., 1996, 105, 10332–10348 CrossRef CAS .
  66. S. Carter, J. M. Bowman and N. C. Handy, Mol. Phys., 2012, 110, 775–781 CrossRef CAS .
  67. Q. Yu and J. M. Bowman, J. Chem. Theory Comput., 2016, 12, 5284 CrossRef CAS .
  68. J. M. Bowman, X. Huang, N. C. Handy and S. Carter, J. Phys. Chem. A, 2007, 111, 7317–7321 CrossRef CAS .
  69. F. Ito and T. Nakanaga, Chem. Phys. Lett., 2000, 318, 571–577 CrossRef CAS .
  70. F. Ito and T. Nakanaga, Chem. Phys., 2002, 277, 163–169 CrossRef CAS .
  71. G. M. Florio, T. S. Zwier, E. M. Myshakin, K. D. Jordan and E. L. Sibert, J. Chem. Phys., 2003, 118, 1735–1746 CrossRef CAS .
  72. R. Georges, M. Freytes, D. Hurtmans, I. Kleiner, J. Vander Auwera and M. Herman, Chem. Phys., 2004, 305, 187–196 CrossRef CAS .
  73. P. Zielke and M. A. Suhm, Phys. Chem. Chem. Phys., 2007, 9, 4528–4534 RSC .
  74. G. L. Barnes and E. L. Sibert, J. Mol. Spectrosc., 2008, 249, 78–85 CrossRef CAS .
  75. Y. H. Yoon, M. L. Hause, A. S. Case and F. F. Crim, J. Chem. Phys., 2008, 128, 084305 CrossRef PubMed .
  76. Ö. Birer and M. Havenith, Annu. Rev. Phys. Chem., 2009, 60, 263–275 CrossRef .
  77. Z. Xue and M. A. Suhm, J. Chem. Phys., 2009, 131, 054301 CrossRef CAS .
  78. M. W. Nydegger, W. Rock and C. M. Cheatum, Phys. Chem. Chem. Phys., 2011, 13, 6098–6104 RSC .
  79. F. Kollipost, R. W. Larsen, A. V. Domanskaya, M. Nörenberg and M. A. Suhm, J. Chem. Phys., 2012, 136, 151101 CrossRef CAS PubMed .
  80. K. Mackeprang, Z. H. Xu, Z. Maroun, M. Meuwly and H. G. Kjaergaard, Phys. Chem. Chem. Phys., 2016, 18, 24654–24662 RSC .
  81. Y. Zhang, W. Li, W. Luo, Y. Zhu and C. Duan, J. Chem. Phys., 2017, 146, 244306 CrossRef .
  82. N. Shida, P. F. Barbara and J. E. Almlöf, J. Chem. Phys., 1991, 94, 3633 CrossRef CAS .
  83. Y. Kim, J. Am. Chem. Soc., 1996, 118, 1522–1528 CrossRef CAS .
  84. T. Loerting and K. R. Liedl, J. Am. Chem. Soc., 1998, 120, 12595–12600 CrossRef CAS .
  85. S. Miura, M. E. Tuckerman and M. L. Klein, J. Chem. Phys., 1998, 109, 5290–5299 CrossRef CAS .
  86. M. V. Vener, O. Kühn and J. M. Bowman, Chem. Phys. Lett., 2001, 349, 562–570 CrossRef CAS .
  87. H. Ushiyama and K. Takatsuka, J. Chem. Phys., 2001, 115, 5903–5912 CrossRef CAS .
  88. C. S. Tautermann, A. F. Voegele and K. R. Liedl, J. Chem. Phys., 2004, 120, 631 CrossRef CAS .
  89. Z. Smedarchina, A. Fernández-Ramos and W. Siebrand, Chem. Phys. Lett., 2004, 395, 339–345 CrossRef CAS .
  90. P. R. L. Markwick, N. L. Doltsinis and D. Marx, J. Chem. Phys., 2005, 122, 054112 CrossRef .
  91. F. Fillaux, Chem. Phys. Lett., 2005, 408, 302–306 CrossRef CAS .
  92. G. V. Mil'nikov, O. Kühn and H. Nakamura, J. Chem. Phys., 2005, 123, 074308 CrossRef .
  93. D. Luckhaus, J. Phys. Chem. A, 2006, 110, 3151–3158 CrossRef CAS .
  94. C. Burisch, P. R. L. Markwick, N. L. Doltsinis and J. Schlitter, J. Chem. Theory Comput., 2008, 4, 164–172 CrossRef CAS .
  95. I. Matanović, N. Došlić and B. R. Johnson, J. Chem. Phys., 2008, 128, 084103 CrossRef .
  96. D. Luckhaus, Phys. Chem. Chem. Phys., 2010, 12, 8357–8361 RSC .
  97. S. D. Ivanov, I. M. Grant and D. Marx, J. Chem. Phys., 2015, 143, 124304 CrossRef PubMed .
  98. A. Jain and E. L. Sibert, J. Chem. Phys., 2015, 142, 084115 CrossRef PubMed .
  99. J. O. Richardson, Phys. Chem. Chem. Phys., 2017, 19, 966–970 RSC .
  100. N.-T. Van-Oanh, C. Falvo, F. Calvo, D. Lauvergnat, M. Basire, M.-P. Gaigeot and P. Parneix, Phys. Chem. Chem. Phys., 2012, 14, 2381–2390 RSC .
  101. T. Esser, H. Knorke, K. R. Asmis, W. Schöllkopf, Q. Yu, C. Qu, J. M. Bowman and M. Kaledin, J. Phys. Chem. Lett., 2018, 9, 798–803 CrossRef CAS .
  102. J. E. Bertie and K. H. Michaelian, J. Chem. Phys., 1982, 76, 886–894 CrossRef CAS .
  103. J. E. Bertie, K. H. Michaelian, H. H. Eysel and D. Hager, J. Chem. Phys., 1986, 85, 4779 CrossRef CAS .
  104. C. Qu, R. Prosmiti and J. M. Bowman, Theor. Chem. Acc., 2013, 132, 1413 Search PubMed .
  105. Q. Yu and J. M. Bowman, J. Chem. Phys., 2017, 146, 121102 CrossRef .
  106. Q. Yu and J. M. Bowman, J. Am. Chem. Soc., 2017, 139, 10984–10987 CrossRef CAS .
  107. C. H. Duong, O. Gorlova, N. Yang, P. J. Kelleher, M. A. Johnson, A. B. McCoy, Q. Yu and J. M. Bowman, J. Phys. Chem. Lett., 2017, 8, 3782–3789 CrossRef CAS .
  108. S. Yoshioki, J. Mol. Graphics Modell., 2007, 25, 856–869 CrossRef CAS PubMed .
  109. S. Yoshioki, J. Mol. Graphics Modell., 2008, 26, 1353–1364 CrossRef CAS PubMed .
  110. J. Reimers and R. Watts, Chem. Phys., 1984, 85, 83–112 CrossRef CAS .
  111. N. Sahu and S. R. Gadre, J. Chem. Phys., 2015, 142, 014107 CrossRef .
  112. T. Salmi, E. Sälli and L. Halonen, J. Phys. Chem. A, 2012, 116, 5368–5374 CrossRef CAS .
  113. Y. Wang, X. Huang, B. C. Shepler, B. J. Braams and J. M. Bowman, J. Chem. Phys., 2011, 134, 094509 CrossRef .
  114. J. S. Mancini and J. M. Bowman, J. Phys. Chem. Lett., 2014, 5, 2247–2453 CrossRef CAS PubMed .
  115. C. R. Jacob and M. Reiher, J. Chem. Phys., 2009, 130, 084106 CrossRef PubMed .
  116. P. T. Panek and C. R. Jacob, ChemPhysChem, 2014, 15, 3365–3377 CrossRef CAS PubMed .
  117. P. T. Panek and C. R. Jacob, J. Chem. Phys., 2016, 144, 164111 CrossRef PubMed .
  118. C. König, M. B. Hansen, I. H. Godtliebsen and O. Christiansen, J. Chem. Phys., 2016, 144, 074108 CrossRef PubMed .
  119. D. Madsen, O. Christiansen and C. König, Phys. Chem. Chem. Phys., 2018, 20, 3445–3456 RSC .
  120. X. Cheng and R. P. Steele, J. Chem. Phys., 2014, 141, 104105 CrossRef PubMed .
  121. S. E. Brown, A. W. Götz, X. Cheng, R. P. Steele, V. A. Mandelshtam and F. Paesani, J. Am. Chem. Soc., 2017, 139, 7082–7088 CrossRef CAS PubMed .
  122. P. M. Zimmerman and P. Smereka, J. Chem. Theory Comput., 2016, 12, 1883–1891 CrossRef CAS PubMed .
  123. A. Molina, P. Smereka and P. M. Zimmerman, J. Chem. Phys., 2016, 144, 124111 CrossRef PubMed .
  124. M. W. D. Hanson-Heine, J. Chem. Phys., 2015, 143, 164104 CrossRef .
  125. H. Liu, Y. Wang and J. M. Bowman, J. Am. Chem. Soc., 2014, 136, 5888–5891 CrossRef CAS PubMed .
  126. H. Liu, Y. Wang and J. M. Bowman, J. Phys. Chem. B, 2016, 120, 1735–1742 CrossRef CAS PubMed .
  127. H. Liu, Y. Wang and J. M. Bowman, J. Phys. Chem. B, 2016, 120, 2824–2828 CrossRef CAS PubMed .
  128. Q. Yu and J. M. Bowman, Mol. Phys., 2015, 113, 3964–3971 CrossRef CAS .
  129. E. D. Sloan and C. A. Koh, Clathrate Hydrates of Natural Gases, CRC Press, Taylor & Francis Group, 2008, pp. 1–44 Search PubMed .
  130. A. K. Sum, R. C. Burruss and E. D. Sloan, J. Phys. Chem. B, 1997, 101, 7371–7377 CrossRef CAS .
  131. C. A. Tulk, J. A. Ripmeester and D. D. Klug, Ann. N. Y. Acad. Sci., 2000, 912, 859–872 CrossRef CAS .
  132. S. Subramanian and E. D. Sloan, J. Phys. Chem. B, 2002, 106, 4348–4355 CrossRef CAS .
  133. K. C. Hester, R. M. Dunk, S. N. White, P. G. Brewer, E. T. Peltzer and E. D. Sloan, Geochim. Cosmochim. Acta, 2007, 71, 2947–2959 CrossRef CAS .
  134. B. Chazallon, C. Focsa, J.-L. Charlou, C. Bourry and J.-P. Donval, Chem. Geol., 2007, 244, 175–185 CrossRef CAS .
  135. H. Ohno, M. Kida, T. Sakurai, Y. Iizuka, T. Hondoh, H. Narita and J. Nagao, ChemPhysChem, 2010, 11, 3070–3073 CrossRef CAS .
  136. J. Qin and W. F. Kuhs, AIChE J., 2013, 59, 2155–2167 CrossRef CAS .
  137. E. Dartois and D. Deboffle, Astron. Astrophys., 2008, 490, L19–L22 CrossRef CAS .
  138. E. Dartois, D. Deboffle and M. Bouzit, Astron. Astrophys., 2010, 514, A49 CrossRef .
  139. R. Martos-Villa, M. Francisco-Márquez, M. P. Mata and C. I. Sainz-Díaz, J. Mol. Graphics Modell., 2013, 44, 253–265 CrossRef CAS .
  140. X. Cao, Y. Su, Y. Liu, J. Zhao and C. Liu, J. Phys. Chem. A, 2014, 118, 215–222 CrossRef CAS .
  141. Y. Liu and L. Ojamäe, J. Phys. Chem. C, 2015, 119, 17084–17091 CrossRef CAS .
  142. H. Itoh and K. Kawamura, Ann. N. Y. Acad. Sci., 2000, 912, 693–701 CrossRef CAS .
  143. J. S. Tse, J. Supramol. Chem., 2002, 2, 429–433 CrossRef CAS .
  144. J. A. Greathouse, R. T. Cygan and B. A. Simmons, J. Phys. Chem. B, 2006, 110, 6428–6431 CrossRef CAS .
  145. M. Hiratsuka, R. Ohmura, A. K. Sum and K. Yasuoka, J. Chem. Phys., 2012, 136, 044508 CrossRef .
  146. M. Hiratsuka, R. Ohmura, A. K. Sum and K. Yasuoka, J. Chem. Phys., 2012, 137, 144306 CrossRef PubMed .
  147. C. Qu, R. Conte, P. L. Houston and J. M. Bowman, Phys. Chem. Chem. Phys., 2015, 17, 8172 RSC .
  148. R. Conte, C. Qu and J. M. Bowman, J. Chem. Theory Comput., 2015, 11, 1631–1638 CrossRef CAS .
  149. M. J. Deible, O. Tuguldur and K. D. Jordan, J. Phys. Chem. B, 2014, 118, 8257–8263 CrossRef CAS .
  150. S. N. Yurchenko, J. Tennyson, R. J. Barber and W. Thiel, J. Mol. Spectrosc., 2013, 291, 69–76 CrossRef CAS .
  151. H. Partridge and D. W. Schwenke, J. Chem. Phys., 1997, 106, 4618–4639 CrossRef CAS .
  152. J. B. Anderson, J. Chem. Phys., 1975, 63, 1499–1503 CrossRef CAS .
  153. J. B. Anderson, J. Chem. Phys., 1976, 65, 4121–4127 CrossRef CAS .
  154. I. Kosztin, B. Faber and K. Schulten, Am. J. Phys., 1996, 64, 633–644 CrossRef .
  155. M. Thomas, M. Brehm, R. Fligg, P. Vohringer and B. Kirchner, Phys. Chem. Chem. Phys., 2013, 15, 6608–6622 RSC .
  156. W. H. Miller, J. Phys. Chem. A, 2001, 105, 2942–2955 CrossRef CAS .
  157. J. Liu, W. H. Miller, F. Paesani, W. Zhang and D. A. Case, J. Chem. Phys., 2009, 131, 164509 CrossRef PubMed .
  158. S. Habershon, G. S. Fanourgakis and D. E. Manolopoulos, J. Chem. Phys., 2008, 129, 074501 CrossRef PubMed .
  159. M. Rossi, M. Ceriotti and D. E. Manolopoulos, J. Chem. Phys., 2014, 140, 234116 CrossRef .
  160. F. Paesani and G. A. Voth, J. Chem. Phys., 2010, 132, 014105 CrossRef PubMed .
  161. G. R. Medders and F. Paesani, J. Chem. Theory Comput., 2015, 11, 1145–1154 CrossRef CAS PubMed .
  162. G. R. Medders and F. Paesani, J. Chem. Phys., 2015, 142, 212411 CrossRef PubMed .

This journal is © the Owner Societies 2019