Article pubs.acs.org/JPCB
Vibrational Analysis of an Ice Ih Model from 0 to 4000 cm−1 Using the Ab Initio WHBB Potential Energy Surface Hanchao Liu, Yimin Wang, and Joel M. Bowman* Cherry L. Emerson Center for Scientific Computation and Department of Chemistry, Emory University, Atlanta, Georgia 30322, United States ABSTRACT: We present an analysis of the vibrational modes of a model of hexagonal ice, ice Ih, comprised of 192 monomers with a core region of 105 monomers, using the ab initio WHBB potential energy surface [Wang, Y.; Shepler, B.; Braams, B.; Bowman, J. M. J. Chem. Phys. 2011, 134, 094509]. A standard normal-mode analysis and a local-monomer normal-mode analysis of 105 core monomers are performed to obtain harmonic frequencies and state densities of the “pseudo-translation” (0−400 cm−1), “libration” (500−1100 cm−1), monomer bend fundamental (∼1600 cm), and O−H stretch (∼3000−3700 cm−1) bands. In addition, the coupled local-monomer model is used to obtain the vibrational density of states in the bend fundamental and O−H stretch regions. The harmonic and local-monomer vibrational density of states obtained from core monomers are in good agreement with those of inelastic neutron scattering spectra, especially the latter, which accounts for anharmonic coupling of monomer modes. Full deuteration is also considered, and the vibrational density of states is again compared to experiment, where good agreement is found.
1. INTRODUCTION Understanding the vibrational dynamics of ice is one of the central goals in ice research. The complexity in structures and molecular interactions, especially through hydrogen-bond networks, makes this study challenging for both experimentalists and theorists. Experimentally, infrared and Raman spectroscopy are the two most popular techniques to probe the vibrational properties of ice.1−5 These spectroscopies are especially useful for studying the intramolecular vibrations of ice.4,5 However, both techniques are limited in the lowfrequency region.1−3 This is due to IR and Raman selection rules specifically for ice,6 and in general, the intensities are governed by the well-known matrix elements and selection rules, and therefore, they do not “report” the entire vibrational density of states (DOS). Inelastic neutron scattering (INS) spectroscopy can, in principle, directly provide the complete vibrational DOS and thus is a more direct tool for obtaining this important quantity. The INS spectra of several phases of ice have been measured,6−15 including hexagonal ice (ice Ih),6−10 which is the focus of this paper. Previous interpretations of the INS spectra were based on calculations using empirical model water potentials, mostly using rigid monomers, from the 1970s to the 2000s.16−22 Recently, density functional theory calculations (with several functionals) were undertaken by Li and co-workers to obtain the vibrational DOS of ice Ih in the range of 0−3500 cm−1.23 These calculations used a 64-monomer supercell and, in part, were aimed at addressing (and contradicting) the conjecture of two types of hydrogen bonding in ice proposed in the 1990s.9,10 These calculations “overestimate” the vibrational © 2013 American Chemical Society
frequencies, that is, they are systematically blue-shifted by roughly 100 cm−1 relative to experiment. In addition, the calculations show a doublet structure in the monomer bend and OH stretch regions, roughly 1600 and 2800−3300 cm−1, respectively, in qualitative disagreement with experiment for those bands. Shortly after those calculations were reported, Hirata, Xantheas, and co-workers reported MP2/aug-cc-pVDZ calculations of the harmonic vibrational DOS (and IR and Raman spectra).24 They also employed a supercell of 64 monomers, which was subjected to an external self-consistent charged environment using the binary interaction method (BIM).25,26 In this method, the electronic energy of water monomers is obtained by one- and two-body interactions in a periodic lattice, perturbed by the electric field of surrounding charges. These calculations obtained the geometry of ice Ih in very good agreement with the value inferred from experiment. In addition, the calculated vibrational DOS is overall in good agreement with experiment, with the exception of not getting the sharp peak in the experimental spectrum at roughly 60 cm−1. As expected with the harmonic approximation, the peaks in the monomer bend and OH-stretching regions are substantially blue-shifted relative to those in experiment. As these authors noted, there are several sources of error in these state-of-the-art, direct calculations. To quote, these are the “inherent error in MP2 theory, a moderate basis size, embedded-fragment scheme and quasi-periodicity of the proton Received: June 13, 2013 Revised: August 1, 2013 Published: August 7, 2013 10046
dx.doi.org/10.1021/jp405865c | J. Phys. Chem. B 2013, 117, 10046−10052
The Journal of Physical Chemistry B
Article
WHBB PES. A core region of 105 monomers was also considered, as described briefly below and in detail elsewhere.27 This core region does not contain any free OH stretch modes, whereas the surface of the full 192-monomer cluster does. Thus, both models are of interest; however, here, the focus is on comparisons with experiment for bulk ice, and therefore, the core model is the relevant one. Vibrational DOS. With the ab initio potential energy surface in hand, performing a normal-mode analysis is a straightforward, almost trivial task compared to a direct ab initio calculation of the force constant matrix (and subsequent trivial diagonalization), for example, the recent one carried out by He et al.24 A standard finite-difference approximation was used to evaluate the mass-weighted Hessian of the 192-mer and 105-mer clusters. For both calculations, every monomer is included in the Hessian calculation; therefore, the matrix dimension is 1728 for the former and 945 for the latter cluster. Diagonalization of the mass-weighted Hessians generates 1728 and 945 harmonic frequencies and corresponding massweighted normal-mode eigenvectors, respectively. The fulldimensional normal-mode analysis took roughly 20 h on 16 processors of a single-node workstation for the 192-mer and 6 h for the 105-mer. Most of the CPU time is in the calculation of the potential. The LMon calculation starts with a local normal-mode (“lormal”) analysis for each monomer.31 This analysis is done for a given monomer by fixing all other monomers at their equilibrium positions. Because the LMon analysis is done for the high-frequency bend and O−H stretch regions, the core 105-monomer model was used,as this does not contain any free O−H stretch modes. The lormal-mode Hessian for each of the monomers is nine-dimensional, and diagonalization of each Hessian gives six “frustrated” translational and rotational modes and three intramolecular vibrational modes. This calculation is, as expected, much faster than diagonalizing the full Hessian. For a single monomer, the CPU time is 15 s on 16 CPUs on a single node, and the calculation is obviously trivially parallel. Again, the CPU time is dominated by the evaluation of the potential. While the results of the harmonic lormal-mode analysis are of potential interest, they are used here for subsequent coupledanharmonic analysis. This postharmonic vibrational analysis was performed by solving the three-intramolecular-mode coupled Schroedinger equation for each monomer.31 This was done using the code MULTIMODE,33−36 with a slight modification to permit parallel computation. As described in detail elsewhere,33−36 the coupled Schroedinger equation is solved using a “configuration interaction” approach, in which the basis functions are the eigenstates states of the zero-point vibrational self-consistent field (VSCF) Hamiltonian for each mode. A small basis, that is, 92 basis functions, is sufficient to obtain the energies of the zero-point, the fundamentals, and bend overtone to within a few wavenumbers or less. The expense of a calculation for one monomer is roughly 90 min on a single CPU, where, again, the major effort is in the evaluation of the many-monomer potential. The evaluations were done on moderate quadrature grids, which result in roughly 1000 potential evaluations per monomer.
positions”. In addition, the authors noted the lack of anharmonicity in their calculations. Other important approximations of the model include (1) using a modest size cluster of 64 water monomers and (2) not explicitly including a threebody potential. This group also focused on the two now“infamous” experimental peaks in the INS at 229 and 306 cm−1 and also concluded that they are not signatures of two types of hydrogen bonding. Inspired by the steps taken by Hirata, Xantheas, and coworkers, we were motivated to extend our recent study of the IR spectra of ice Ih and amorphous solid water in the monomer bend and “O−H stretch” regions27 to a calculation of the vibrational DOS of a model for ice Ih. As with the earlier calculation, we used a 192-monomer model with the ab initio WHBB potential energy surface,28 which is briefly reviewed below. A group of 105 core monomers was selected27 and is also used here for the vibrational analysis. The results from this smaller core cluster are considered under the assumption that they are the ones directly comparable to the INS experimental data. Specifically, we calculated the vibrational DOS in the range of 0−4000 cm−1 using standard and local harmonic normal-mode analysis and also using the coupled-anharmonic local-monomer (LMon) model, also reviewed briefly below, in the range of 1500−4000 cm−1. In addition, we calculated the vibrational DOS for pure D2O and compared both spectra with experimental INS spectra. The character of the vibrational states in different frequency ranges was analyzed from the normal-mode eigenvectors. A short summary and conclusions are given in the last section. The paper is organized as follows. The Computational Details are given in the next section. This consists of a brief review of the WHBB potential and the 192-model of ice Ih using it. Then details of the harmonic and coupled-anharmonic calculations of the vibrational DOS are given. The results and comparisons with experiment are given in section 3. An analysis of the normal-mode eigenvectors is also given in that section, with the aim of assigning the nature of the modes in four regions of the vibrational DOS spectrum.
2. COMPUTATIONAL DETAILS Potential Energy Surface. The WHBB potential for N water monomers is given by the following one-, two-, and three-body (monomer) expansion28 N
V=
N
N
∑ V1‐body(i)+ ∑ V2‐body(i , j) + ∑ i=1
i
