Insights in quantum dynamical effects in the infrared spectroscopy of liquid water from a semiclassical study with an ab initio-based flexible and polarizable force field

The dynamical properties of liquid water play an important role in many processes in nature. In this paper, we focus on the infrared (IR) absorption spectrum of liquid water based on the linearized semiclassical initial value representation (LSC-IVR) with the local Gaussian approximation (LGA) [J. Liu and W. H. Miller, J. Chem. Phys. 131, 074113 (2009)] and an ab initio based, flexible, polarizable Thole-type model (TTM3-F) [G. S. Fanourgakis and S. S. Xantheas, J. Chem. Phys. 128, 074506 (2008)]. Although the LSC-IVR (LGA) gives the exact result for the isolated three-dimensional shifted harmonic stretching model, it yields a blueshifted peak position for the more realistic anharmonic stretching potential. By using the short-time information of the LSC-IVR correlation function; however, it is shown how one can obtain more accurate results for the position of the stretching peak. Due to the physical decay in the condensed phase system, the LSC-IVR (LGA) is a good and practical approximate quantum approach for the IR spectrum of liquid water. The present results offer valuable insight into future attempts to improve the accuracy of the TTM3-F potential or other ab initio-based models in reproducing the IR spectrum of liquid water.     

Using the thermal Gaussian approximation for the Boltzmann operator in semiclassical initial value time correlation functions

The thermal Gaussian approximation (TGA) recently developed by Frantsuzov et al. [Chem. Phys. Lett. 381, 117 (2003)] has been demonstrated to be a practical way for approximating the Boltzmann operator exp(-betaH) for multidimensional systems. In this paper the TGA is combined with semiclassical (SC) initial value representations (IVRs) for thermal time correlation functions. Specifically, it is used with the linearized SC-IVR (LSC-IVR, equivalent to the classical Wigner model), and the "forward-backward semiclassical dynamics" approximation developed by Shao and Makri [J. Phys. Chem. A 103, 7753 (1999); 103, 9749 (1999)]. Use of the TGA with both of these approximate SC-IVRs allows the oscillatory part of the IVR to be integrated out explicitly, providing an extremely simple result that is readily applicable to large molecular systems. Calculation of the force-force autocorrelation for a strongly anharmonic oscillator demonstrates its accuracy, and calculation of the velocity autocorrelation function (and thus the diffusion coefficient) of liquid neon demonstrates its applicability.     

Real time correlation function in a single phase space integral beyond the linearized semiclassical initial value representation

It is shown how quantum mechanical time correlation functions [defined, e.g., in Eq. (1.1)] can be expressed, without approximation, in the same form as the linearized approximation of the semiclassical initial value representation (LSC-IVR), or classical Wigner model, for the correlation function [cf. Eq. (2.1)], i.e., as a phase space average (over initial conditions for trajectories) of the Wigner functions corresponding to the two operators. The difference is that the trajectories involved in the LSC-IVR evolve classically, i.e., according to the classical equations of motion, while in the exact theory they evolve according to generalized equations of motion that are derived here. Approximations to the exact equations of motion are then introduced to achieve practical methods that are applicable to complex (i.e., large) molecular systems. Four such methods are proposed in the paper--the full Wigner dynamics (full WD) and the second order WD based on "Wigner trajectories" [H. W. Lee and M. D. Scully, J. Chem. Phys. 77, 4604 (1982)] and the full Donoso-Martens dynamics (full DMD) and the second order DMD based on "Donoso-Martens trajectories" [A. Donoso and C. C. Martens, Phys. Rev. Lett. 8722, 223202 (2001)]--all of which can be viewed as generalizations of the original LSC-IVR method. Numerical tests of the four versions of this new approach are made for two anharmonic model problems, and for each the momentum autocorrelation function (i.e., operators linear in coordinate or momentum operators) and the force autocorrelation function (nonlinear operators) have been calculated. These four new approximate treatments are indeed seen to be significant improvements to the original LSC-IVR approximation.     

Complex autocorrelation function and energy spectrum by classical trajectory calculations

A quasiclassical method which enables evaluation of complex autocorrelation function from classical trajectory calculations is proposed. The method is applied for two highly excited nonlinearly coupled harmonic oscillators in regimes prevailed either by regular or chaotic classical motions. A good agreement of classical and quantum autocorrelation functions is found within short (Ehrnfest) time limit. Fourier transforms of the autocorrelation functions provide moderate resolved energy spectra, where classical and quantum results nearly coincide. The actual energy levels are obtained from approximate short-time autocorrelation functions with the help of filter diagonalization. This study is a follow up to our previous work [P. Zdanska and N. Moiseyev, J. Chem. Phys. 115, 10608 (2001)], where the complex autocorrelation has been obtained up to overall phase factors of recurrences.     

Unimolecular rovibrational bound and resonance states for large angular momentum: J=20 calculations for HO2

We explore the calculation of unimolecular bound states and resonances for deep-well species at large angular momentum using a Chebychev filter diagonalization scheme incorporating doubling of the autocorrelation function as presented recently by Neumaier and Mandelshtam [Phys. Rev. Lett. 86, 5031 (2001)]. The method has been employed to compute the challenging J=20 bound and resonance states for the HO2 system. The methodology has firstly been tested for J=2 in comparison with previous calculations, and then extended to J=20 using a parallel computing strategy. The quantum J-specific unimolecular dissociation rates for HO2-->H+O2 in the energy range from 2.114 to 2.596 eV have been reported for the first time, and comparisons with the results of Troe and co-workers [J. Chem. Phys. 113, 11019 (2000) Phys. Chem. Chem. Phys. 2, 631 (2000)] from statistical adiabatic channel method/classical trajectory calculations have been made. For most of the energies, the reported statistical adiabatic channel method/classical trajectory rate constants agree well with the average of the fluctuating quantum-mechanical rates. Near the dissociation threshold, quantum rates fluctuate more severely, but their average is still in agreement with the statistical adiabatic channel method/classical trajectory results.     

Test of the consistency of various linearized semiclassical initial value time correlation functions in application to inelastic neutron scattering from liquid para-hydrogen

The linearized approximation to the semiclassical initial value representation (LSC-IVR) is used to calculate time correlation functions relevant to the incoherent dynamic structure factor for inelastic neutron scattering from liquid para-hydrogen at 14 K. Various time correlations functions were used which, if evaluated exactly, would give identical results, but they do not because the LSC-IVR is approximate. Some of the correlation functions involve only linear operators, and others involve nonlinear operators. The consistency of the results obtained with the various time correlation functions thus provides a useful test of the accuracy of the LSC-IVR approximation and its ability to treat correlation functions involving both linear and nonlinear operators in realistic anharmonic systems. The good agreement of the results obtained from different correlation functions, their excellent behavior in the spectral moment tests based on the exact moment constraints, and their semiquantitative agreement with the inelastic neutron scattering experimental data all suggest that the LSC-IVR is indeed a good short-time approximation for quantum mechanical correlation functions.     

Critical comparison of classical and quantum mechanical treatments of the phase equilibria of water

The Gibbs ensemble Monte Carlo simulation technique was used to compare the phase equilibria of the rigid TIP4P water model [Jorgensen et al., J. Chem. Phys. 79, 926 (1983)] utilizing classical and quantum statistical mechanics. The quantum statistical mechanical treatment generally resulted in lower liquid densities and higher vapor densities, narrowing the phase envelope. As a result, the calculated critical temperatures and normal boiling points were lower from the quantum simulations than the classical by 22 and 17 K, respectively, but the critical densities were equal within the estimated uncertainties. When the phase diagram from the quantum statistical mechanical treatment was increased by 22 K, it agreed with the classical results quite well throughout the entire simulated temperature range. A semiclassical treatment, involving a low order expansion in Planck's constant, resulted in good agreement with the path integral results for second virial coefficients, but gave densities and vapor pressures that fluctuated between the values for the classical and quantum statistical mechanics values, with no definite agreement with either.     

An approach for generating trajectory-based dynamics which conserves the canonical distribution in the phase space formulation of quantum mechanics. I. Theories

We have reformulated and generalized our recent work [J. Liu and W. H. Miller, J. Chem. Phys. 126, 234110 (2007)] into an approach for generating a family of trajectory-based dynamics methods in the phase space formulation of quantum mechanics. The approach (equilibrium Liouville dynamics) is in the spirit of Liouville's theorem in classical mechanics. The trajectory-based dynamics is able to conserve the quantum canonical distribution for the thermal equilibrium system and approaches classical dynamics in the classical (? → 0), high temperature (β → 0), and harmonic limits. Equilibrium Liouville dynamics provides the framework for the development of novel theoretical∕computational tools for studying quantum dynamical effects in large∕complex molecular systems.     

Higher excitations for an exponential multireference wavefunction Ansatz and single-reference based multireference coupled cluster Ansatz: application to model systems H(4), P(4), and BeH(2)

This article reports on the convergence of the exponential multireference wavefunction Ansatz (MRexpT) [J. Chem. Phys. 123, 84102 (2005)] and the single-reference based multireference coupled cluster Ansatz [J. Chem. Phys. 94, 1229 (1991)] with respect to higher cluster excitations. The approaches are applied to the H(4), P(4), and BeH(2) model systems according to the recently published analysis by Evangelista et al. [J. Chem. Phys. 125, 154113 (2006)]. The results show both MRexpT and SRMRCC to be highly accurate although SRMRCC shows problems due to its lack of Fermi vacuum invariance (symmetry breaking).     

Calculations of molecular properties in hybrid coupled-cluster and molecular mechanics approach

We report benchmark calculations obtained with our new coupled-cluster singles and doubles (CCSD) code for calculating the first- and second-order molecular properties. This code can be easily incorporated into combined [Valiev, M.; Kowalski, K. J. Chem. Phys. 2006, 125, 211101] classical molecular mechanics (MM) and ab initio coupled-cluster (CC) calculations using NWChem, enabling us to study molecular properties in a realistic environment. To test this methodology, we discuss the results of calculations of dipole moments and static polarizabilities for the Cl2O system in the CCl4 solution using the CCSD (CC with singles and doubles) linear response approach. We also discuss the application of the asymptotic extrapolation scheme (AES) [Kowalski, K.; Valiev, M. J. Phys. Chem. A 2006, 110, 13106] in reducing the numerical cost of CCSD calculations.     

Thermodynamics and equilibrium structure of Ne38 cluster: quantum mechanics versus classical

The equilibrium properties of classical Lennard-Jones (LJ38) versus quantum Ne38 Lennard-Jones clusters are investigated. The quantum simulations use both the path-integral Monte Carlo (PIMC) and the recently developed variational-Gaussian wave packet Monte Carlo (VGW-MC) methods. The PIMC and the classical MC simulations are implemented in the parallel tempering framework. The classical heat capacity Cv(T) curve agrees well with that of Neirotti et al. [J. Chem. Phys. 112, 10340 (2000)], although a much larger confining sphere is used in the present work. The classical Cv(T) shows a peak at about 6 K, interpreted as a solid-liquid transition, and a shoulder at approximately 4 K, attributed to a solid-solid transition involving structures from the global octahedral (Oh) minimum and the main icosahedral (C5v) minimum. The VGW method is used to locate and characterize the low energy states of Ne38, which are then further refined by PIMC calculations. Unlike the classical case, the ground state of Ne38 is a liquidlike structure. Among the several liquidlike states with energies below the two symmetric states (Oh and C5v), the lowest two exhibit strong delocalization over basins associated with at least two classical local minima. Because the symmetric structures do not play an essential role in the thermodynamics of Ne38, the quantum heat capacity is a featureless curve indicative of the absence of any structural transformations. Good agreement between the two methods, VGW and PIMC, is obtained. The present results are also consistent with the predictions by Calvo et al. [J. Chem. Phys. 114, 7312 (2001)] based on the quantum superposition method within the harmonic approximation. However, because of its approximate nature, the latter method leads to an incorrect assignment of the Ne38 ground state as well as to a significant underestimation of the heat capacity.     

Closed-shell ring coupled cluster doubles theory with range separation applied on weak intermolecular interactions

We explore different variants of the random phase approximation to the correlation energy derived from closed-shell ring-diagram approximations to coupled cluster doubles theory. We implement these variants in range-separated density-functional theory, i.e., by combining the long-range random phase approximations with short-range density-functional approximations. We perform tests on the rare-gas dimers He(2), Ne(2), and Ar(2), and on the weakly interacting molecular complexes of the S22 set of Jurec?ka et al. [P. Jurec?ka, J. S?poner, J. C?erny?, and P. Hobza, Phys. Chem. Chem. Phys. 8, 1985 (2006)]. The two best variants correspond to the ones originally proposed by Szabo and Ostlund [A. Szabo and N. S. Ostlund, J. Chem. Phys. 67, 4351 (1977)]. With range separation, they reach mean absolute errors on the equilibrium interaction energies of the S22 set of about 0.4 kcal/mol, corresponding to mean absolute percentage errors of about 4%, with the aug-cc-pVDZ basis set.     

Extended hydrodynamic approach to quantum-classical nonequilibrium evolution. II. Application to nonpolar solvation

The mixed quantum-classical formulation derived in our companion paper [D. Bousquet, K. H. Hughes, D. Micha, and I. Burghardt, J. Chem. Phys. 134, 064116 (2011)], which is based upon a hydrodynamic representation of the classical sector, is applied to nonequilibrium nonpolar solvation dynamics as exemplified by the solvation of the electronically excited NO molecule in a rare gas environment. Derived from a partition of the Hamiltonian into a primary (quantum) part and a secondary (classical) part the hydrodynamic equations are formulated for multi-quantum states and result in explicit equations of motion for populations and coherences. The hierarchy of hydrodynamic equations is truncated by the following approximate closure schemes: Gauss-Hermite closure, dynamical density functional theory approximation, and a generalized Maxwellian closure. A comparison of the dynamics using these three closure methods showed that the suitability of a particular closure scheme was dependent on the initial conditions and the nonequilibrium character of the dynamics.     

Classical and quantum gibbs free energies and phase behavior of water using simulation and cell theory

A method to calculate the classical and quantum free energy of a liquid from a computer simulation by using cell theory [J. Chem. Phys. 2007, 126, 064504] is tested for liquid water and ice Ih against experiment as a function of temperature. This fast and efficient method reproduces reasonably well the experimental values of entropy, enthalpy, and free energy of a liquid across the supercooled, stable, and superheated range of temperatures considered. There are small differences between classical and quantum results of water at 298 K, necessitating a small correction term to reproduce water's enthalpy of vaporisation. Only at higher temperatures is entropy underestimated by up to 9 J K(-1) mol(-1) as verified by thermodynamic integration calculations. Satisfactory agreement for ice, however, is only obtained by using the quantum formulation. Even then, at higher temperatures, the entropies exceed experiment by up to 15 J K(-1) mol(-1). Further insight into the quantum nature of water is provided by inspecting the temperature dependence of the frequencies. The harmonic approximation is further supported by the harmonic force and torque distributions and the very similar entropies obtained from the force and torque variances. All these results suggest that the single molecule harmonic oscillator approximation for water, although not exact, provides a rapid, insightful, and useful means to evaluate the thermodynamic properties of water from a computer simulation in a way that can account for the quantization of water's energy levels.     

Thermochemistry of new molecular species: SBr and HSBr

Total energies, optimized geometries, and vibrational frequencies of SBr and HSBr have been evaluated at the coupled cluster level of theory with the correlation consistent basis sets. Extrapolated to the complete basis set limit and with corrections for core-valence, scalar relativistic, and spin-orbit effects, atomization energies were computed and then combined with the experimental heats of formation of the atomic species to generate very accurate heats of formation for the species SBr and HSBr. For SBr, we predict 37.45 and 36.07 kcal/mol for DeltaHf(0 K) and DeltaHf (298.15 K), respectively, in very good agreement with the inferred experimental values of 37.98 and 36.15 kcal/mol. For HSBr, the estimate turns out to be 10.38 and 8.29 kcal/mol for DeltaHf (0 K) and DeltaHf (298.15 K), respectively. Using the more recent HBrO experimental heat of formation at 298.15 K of Lock et al., [J. Phys. Chem. 100, 7972 (1996)] the inferred experimental value for HSBr is predicted to be 8.15 kcal/mol, compared with 8.65 kcal/mol derived from the data of Ruscic and Berkowitz [J. Chem. Phys. 101, 7795 (1994)]. Considering the better agreement of the result with that predicted using the experimental value of DeltaHf(298.15 K) of Lock et al., the author also supports the suggestion made by Denis [J. Phys. Chem. A. 110, 5887 (2006)] that the result of Lock et al. should be preferred over the one of Ruscic and Berkowitz. For DeltaHf(0 K), the author found 10.38 and 10.56 kcal/mol, respectively, for the theoretical and inferred experimental estimates.     

Molecular nitrogen-N2 properties: the intermolecular potential and the equation of state

Quantum mechanical (QM) high precision calculations were used to determine N(2)-N(2) intermolecular interaction potential. Using QM numerical data the anisotropic potential energy surface was obtained for all orientations of the pair of the nitrogen molecules in the rotation invariant form. The new N(2)-N(2) potential is in reasonably good agreement with the scaled potential obtained by van der Avoird et al. using the results of Hartree-Fock calculations [J. Chem. Phys. 84, 1629 (1986)]. The molecular dynamics (MD) of the N(2) molecules has been used to determine nitrogen equation of state. The classical motion of N(2) molecules was integrated in rigid rotor approximation, i.e., it accounted only translational and rotational degrees of freedom. Fincham [Mol. Simul. 11, 79 (1993)] algorithm was shown to be superior in terms of precision and energy stability to other algorithms, including Singer [Mol. Phys. 33, 1757 (1977)], fifth order predictor-corrector, or Runge-Kutta, and was therefore used in the MD modeling of the nitrogen pressure [S. Krukowski and P. Strak, J. Chem. Phys. 124, 134501 (2006)]. Nitrogen equation of state at pressures up to 30 GPa (300 kbars) and temperatures from the room temperature to 2000 K was obtained using MD simulation results. Results of MD simulations are in very good agreement (the error below 1%) with the experimental data on nitrogen equation of state at pressures below 1 GPa (10 kbars) for temperatures below 1800 K [R. T. Jacobsen et al., J. Phys. Chem. Ref. Data 15, 735 (1986)]. For higher temperatures, the deviation is slightly larger, about 2.5% which still is a very good agreement. The slightly larger difference may be attributed to the vibrational motion not accounted explicitly by rigid rotor approximation, which may be especially important at high temperatures. These results allow to obtain reliable equation of state of nitrogen for pressures up to 30 GPa (300 kbars), i.e., close to molecular nitrogen stability limit, determined by Nellis et al. [Phys. Rev. Lett. 53, 1661 (1984)].     

Rotational excitation of water by hydrogen molecules: comparison of results from classical and quantum mechanics

Quasiclassical trajectory calculations are carried out for rotational excitation of water by hydrogen molecules. State-to-state rate coefficients are determined at 100 K and are compared to available quantum results. A good agreement between classical and quantum rates is observed for downward transitions, with an average accuracy of classical results better than a factor of 2. It is thus found that the ambiguities described by Faure and Wiesenfeld [J. Chem. Phys. 121, 6771 (2004)] can be solved in the particular case of waterlike asymmetric-top molecules.     

Symmetry broken self-consistent field approximation for coupled oscillators

A self-consistent field approximation allowing for a change of coordinates is applied to a model of two coupled oscillators considered by Davis and Heller [J. Chem. Phys. 75 , 246 (1981)]. The model accounts well for the unsymmetrical nonstationary states involved in the quantum dynamical tunneling phenomenon.     

Quasi-classical trajectories study of Ne79Br2(B) vibrational predissociation

A full-dimensional quasi-classical trajectories study on the vibrational predissociation (VP) of the Ne79Br2(B) complex is presented. Following the most recent experiments, the Br2(B) vibrational levels v'=16-29 were explored. The total angular momentum, J, was taken to be zero, and a semiclassical Franck-Condon model to compute initial conditions from quantum distributions was employed. Predissociation lifetimes were extracted from Ne79Br2 population decay by using two different exponential laws. Predicted lifetimes are in excellent agreement with the last experimental results [J. A. Cabrera, C. R. Bieler, B. C. Olbricht, W. E. van der Veer and K. C. Janda, J. Chem. Phys., 2005, 123, 054311]. The Br2 fragment ro-vibrational distributions resulting from the VP of the molecule were obtained from the statistics of classical magnitudes using the standard binning procedure. Computed rotational distributions (for the Deltav'=-1, -2 channels) are also in very good agreement with the experimental results [M. Nejad-Sattari and T. A. Stephenson, J. Chem. Phys., 1997, 106 5454]. The influence of two quantum effects-the closing of the Deltav'=-1 dissociation channel and the intramolecular vibrational relaxation (IVR) mechanism-on the agreement with experimental rotational distributions, is discussed. Due to the classical character of our calculations and the binning procedure we used, the agreement of computed vibrational distributions with experimental and quantum theoretical is qualitative. For instance, for v'=28-for which the Deltav'=-1 channel is experimentally found to be closed-the Deltav'=-2 channel becomes statistically more significant. A discussion on the viability of similar quasi-classical methods to model the VP dynamics of analogous clusters is presented.