- Email: [email protected]
Analysis of librational modes of ice XI studied by Car–Parrinello molecular dynamics
Chemical Physics 459 (2015) 102–111
Contents lists available at ScienceDirect
Chemical Physics journal homepage: www.elsevier.com/locate/chemphys
Analysis of librational modes of ice XI studied by Car–Parrinello molecular dynamics Maciej Gług, Marek Boczar ⇑, Łukasz Boda, Marek J. Wójcik ⇑ Faculty of Chemistry, Jagiellonian University, 30-060 Kraków, Ingardena 3, Poland
a r t i c l e
i n f o
Article history: Received 22 May 2015 In final form 9 August 2015 Available online 14 August 2015 Keywords: Ice XI IR spectra Car–Parrinello molecular dynamics
a b s t r a c t Theoretical infrared spectra of ice XI at 4 and 40 K were obtained by Car–Parrinello molecular dynamics calculations and evaluating the dipole moment using maximally-localized Wannier functions. In order to improve the agreement between theoretical and experimental spectra, the Kramers–Kronig relations were applied. Based on analysis of power and polarization spectra the assignment of modes was carried out. The results were compared with experimental and theoretical spectra of ice XI, performed by other authors. Such simulations can be very useful while investigating the existence of ice XI in space. Ó 2015 Elsevier B.V. All rights reserved.
1. Introduction Ice XI, discovered in 1972 by Kawada [1] and Tajima et al. [2,3], is proton-ordered, low-temperature phase of ice Ih (ordinary hexagonal ice). Few years later it was found, using neutron diffraction, that ice XI belongs to the orthorhombic Cmc21 symmetry space group with eight water molecules per unit cell [4]. The ferroelectric character of this structure was subsequently confirmed by thermal stimulated depolarization experiments [5]. Ice XI and its ferroelectricity have attracted scientific interest mainly after discovering the presence of crystalline ice in outer solar system [6]. Long-range electrostatic forces, caused by the ferroelectricity, might be an important factor for planet formation [7–9]. Fukazawa et al. [10,11] presented hypothesis that ice XI can exist in interstellar space. Electrostatic force caused by the ferroelectric properties might affect agglomeration of ice particles in the space [7,12,13]. Such presence of ice XI might be very significant for planetary formation processes. Infrared telescope measurements do not provide clear information whether ice XI exists in space or not. This is mainly due to lack of reference data. There is a strong need to obtain good-quality IR spectra of ice XI, either experimental or theoretical. Due to the high absorptivities and difficulties in preparation of pure samples it is difficult to obtain experimental infrared spectra of ice XI. For this reason, theoretical study of IR spectra of ice XI, presented in this paper, is very important to provide reference spectra for astrophysical measurements. ⇑ Corresponding authors. Fax: +48 12 634 05 15. E-mail addresses: [email protected] (M. Boczar), [email protected] (M.J. Wójcik). http://dx.doi.org/10.1016/j.chemphys.2015.08.004 0301-0104/Ó 2015 Elsevier B.V. All rights reserved.
Ice XI is the low temperature equilibrium structure of hexagonal ice (ice Ih) – it has a triple point with hexagonal ice and gaseous water at 72 K and 0 Pa. The arrangement of water molecules in hexagonal ice, in which they are surrounded by four semi-randomly directed hydrogen bonds, should change to a regular arrangement of hydrogen bonds at low temperatures. This process becomes more effective with an increase in pressure [14]. In nature the transformation is very slow without the presence of catalysts, even though ice XI is thought to have a more stable conformation than ice Ih, mainly due to relatively small energy difference between the two structures [15]. Fukazawa et al. [16] claimed the presence of ice XI in the Antarctic ice. They found that the phase transformation of ice XI to ice Ih occurs at temperature 237 K, which is higher than the temperature of the expected triple point (72 K, 0 Pa). Ice XI can be obtained under laboratory conditions from a dilute (10 mM) KOH solution kept just below 72 K at ambient pressure for about a week [17,18]. The mechanism of the influence of hydroxide and K+ ions on this process was discussed by Suga [19], however many aspects are still not exactly understood. The ice Ih M XI transition was previously reported only for the doped, not for pure water, however ice XI was also found experimentally using pure water at very low temperature (10 K) at low pressure [20]. Such conditions are thought to be present in the upper atmosphere. Ice XI was subject of numerous spectroscopic studies. The vibrational spectra, which reveal the crystal structure and atomic motions, are investigated using mainly infrared, Raman and neutron spectroscopy. In the ice crystal water molecules are connected by a network of hydrogen bonds. In the crystal there are
103
M. Gług et al. / Chemical Physics 459 (2015) 102–111
intermolecular translational lattice and librational vibrations and intramolecular OAH stretching vibrations. Fukazawa et al. [18,21] have postulated that the largest difference in the spectroscopic data between ice XI and ice Ih occurs in the region of librational modes. Infrared absorption spectra of KOH-doped ice in the librational region for different temperatures was measured by Arakawa et al. [22]. Inelastic incoherent neutron scattering of ice XI was reported by Li et al. [23,24], who also studied the ice Ih M ice XI transition using Raman spectroscopy [23]. Phase transition and the differences of the vibrational modes between ice Ih and XI in the translational and librational region was also studied by Raman spectroscopy by Abe et al. [25]. Recently experimental studies of the Raman spectra of hydrogenordered ice have been reported by Abe and Shigenari [26,27]. In our recent paper [28] we presented comparison between theoretical IR spectra of ice Ih and ice XI obtained by Car–Parrinello molecular dynamics (CPMD) calculations [29]. In this article we present the details of the crystal structure of ice XI and analysis of power and polarization spectra of ice XI. Different forms of ice have been studied in the literature. Infrared and Raman spectra of ice Ih, cubic and amorphous ice have been previously calculated using different potentials and theoretical approaches by Rice et al. [30,31], Wójcik et al. [32,33] and Skinner et al. [34,35]. CPMD calculations of vibrational frequencies, phase transitions and dielectric properties in different phases of ice were reported by Benoit et al. [36], Bernasconi et al. [37], Knight et al. [38] and Schönherr et al. [39]. In our study spectra are calculated from the velocity/dipole moment autocorrelation function (it is described in Section 2.2). In this way they are not state resolved. One of approaches for quantization of nuclear motion was used by Brela et al. [40] and Stare et al. [41]. These authors quantize selected nuclear degrees of freedom. This takes into account the influence of the quantum nature of nuclear motion on the resultant spectroscopic properties. We plan to use similar method in the future.
Born–Oppenheimer molecular dynamics, but requires shorter time step [43]. Hence, the choice of the l value should be the compromise between the accuracy and the time of the computation. 2.2. Theoretical infrared spectrum IR absorption coefficient a(x) is commonly calculated using the formula [44]
NðxÞ ¼ nðxÞ þ ikðxÞ;
All of the calculations reported in this paper were carried out within the Car–Parrinello ab initio molecular dynamics method [29]. In this approach the electronic degrees of freedom as (fictitious) dynamical variables are introduced using Car–Parrinello Lagrangian (embedded in Kohn–Sham formalism [42]):
X 1X _2 X _ _ ¼ M I RI þ lh/i j/i i hW0 jHKS Kij ðh/i j/j i dij Þ: e jW0 i þ 2 I i i;j ð1Þ
This leads to the system of coupled equations of motion for both, ions and electrons:
€ I ¼ rI hW0 jHKS jW0 i; MI R e
ð2aÞ
X Kij /j :
ð2bÞ
l/€ i ¼ HKS e /i þ
ð4Þ
where the real part of the refractive index n indicates the phase velocity, whereas the imaginary part k is responsible for the amount of absorption loss, when the electromagnetic wave propagates through the material, and k is related to the absorption coefficient by the relation
aðxÞc : x
ð5Þ
Electric susceptibility v is connected with N index by the formula
2.1. Car–Parrinello molecular dynamics
LCP
ð3Þ
where T is temperature, V is volume of cell, b is 1/(kT), c is speed of light in vacuum, n(x) is refractive index and M is total dipole moment. The brackets indicate autocorrelation function. M(t) function is the sum of electronic and ionic dipole moments computed in every step of the simulation. The ionic contribution is calculated straightforward whereas electrons dipole moment is obtained using maximally-localized Wannier functions [45]. Due to the fact, that unknown frequency dependence of n cannot be neglected, this equation provides only the product of n(x) and a(x). However, it is possible to extract absorption coefficient a(x) by using Kramers–Kronig theorem [46,47] giving the relations between the real and imaginary part of the dielectric function. The relations are of a general nature and are based on the properties of a complex, analytical response function, such as the complex refractive index in our case, defined as
kðxÞ ¼ 2. Computational details
Z 1 2bpx2 expðixtÞhMð0ÞMðtÞidt; 3VcnðxÞ 1
aðxÞ ¼
N2 ðxÞ ¼ 1 þ 4pvðxÞ
Comparing Eqs. (4) and (6) we obtain the following expression for v 2
vðxÞ ¼ v1 ðxÞ þ iv2 ðxÞ ¼
Consequently, it is not necessary to solve electronic Schrödinger Equation at each time step of simulation (opposite to Born–Oppenheimer molecular dynamics). In Eq. (1) the first component is the kinetic energy of the ions, followed by the kinetic energy of the electrons, and the finally potential energy and orbitals orthonormality constrains, respectively. The parameter l is called the fictitious mass of the electron. Reducing of the value of this parameter yields an increased accuracy of electronic energy calculations, with respect to
n2 ðxÞ k ðxÞ 1 nðxÞkðxÞ þ i 4p 2p
ð7Þ
Both real and imaginary parts of electric susceptibility are mutually connected by Kramers–Kronig relations [48]. This leads to the equation 2
v1 ðxÞ ¼
n2 ðxÞ k ðxÞ 1 2 ¼ Pr 4p p
Z
1
0
x0 nðx0 Þkðx0 Þ=2p 0 dx ; x02 x2
ð8Þ
which can be represented using absorption coefficient 2
j
ð6Þ
v1 ðxÞ ¼
2
n2 ðxÞ a ðxx2Þc 1 c ¼ 2 Pr p 4p
Z 0
1
nðx0 Þaðx0 Þ dx0 ; x02 x2
ð9Þ
where Pr denotes the Cauchy principal value of the integral, which is the limit for d ? 0 of the sum of the integrals over 1 < x0 < x d and x + d < x0 < 1. Basing on above formula it can be concluded that knowing the product of absorption coefficient and refractive index for all frequencies allows to calculate real part of susceptibility. Denoting this product as I(x) the system of two equations can be written 2
2
n2 ðxÞ a ðxx2Þc 1 ¼ v1 ðxÞ 4p
ð10Þ
104
M. Gług et al. / Chemical Physics 459 (2015) 102–111
nðxÞaðxÞ ¼ IðxÞ
ð11Þ
The right side of each equation is given by Eqs. (9) and (3) respectively. The solution with respect to n(x) and a(x) is
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u 2 2 u1 þ 4pv1 ðxÞ ð1 þ 4pv1 ðxÞÞ þ 4I2 ðxÞ xc 2 aðxÞ ¼ t 2 2 xc 2
ð12Þ
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 0v qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u 2 1 þ 4pv1 ðxÞ ð1 þ 4pv1 ðxÞÞ2 þ 4I2 ðxÞ xc 2 C Bu t C nðxÞ ¼ IðxÞB 2 A @ 2 xc 2 ð13Þ The finally formula on absorption coefficient is
aðxÞ ¼
x c
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 Z 1 u 1 2c Z 1 Iðx0 Þ 1 4c IðxÞ c2 t 0þ 0 d x P d x þ 4I2 ðxÞ 2 : 1 þ P 2 2 p 0 x02x 2 p 0 x02 x2 x
thermostat [51]. To calculate the cell dipole moment, maximally localized Wannier functions [45] were used. A fictitious mass of 300 a.u. was set for the electronic degrees of freedom. The equations of motion for nuclei and molecular orbitals were solved using the velocity Verlet algorithm with 1.5 a.u. timestep. Electrons for oxygen and hydrogen atoms were treated within the DFT framework with BLYP functional and the norm-conserving pseudopotential of Goedecker type [52–54]. The electronic wave functions were expanded in plane waves. The energy cutoff was adjusted to 230 Ry. Brillouin zone sampling was restricted to the Gamma point. The total length of simulation path after equilibrating was approx. 12 ps (330,000 steps). Periodic boundary conditions were employed for the solid state calculations with the system which theoretically extends to infinity. Power spectrum and atomic power spectra were calculated using the Fourier transformation (FT) of the atomic velocity autocorrelation function. IR spectrum was obtained from FT of the dipole moment autocorrelation function and the application of Kramers–Kronig relations. All FT calculations were performed with CPMD additional script (FOURIER.X).
ð14Þ 2.3. Simulation details Orthorhombic cell, shown in Fig. 1, containing eight water molecules with the parameters [49] included in Table 1 was selected for simulation. All calculations were executed using CPMD package [50]. Car–Parrinello molecular dynamics simulation was performed in thermodynamic equilibrium using Nosé–Hoover chain
Table 1 Positions of atoms for system from Fig. 1. The length of all OAH bond is 0.9572 Å and the value of all angles HAOAH is 104.5°. The cell is orthorhombic with a = 4.49225 Å, b = 7.78080 Å, c = 7.33581 Å [59]. Atom
Fractional coordinates
O1 O2 H1 H2 H3
0.0000 0.5000 0.0000 0.6685 0.0000
0.3333 0.1667 0.4510 0.2282 0.3386
0.5625 0.4375 0.5244 0.4836 0.6929
Fig. 1. Ice XI orthorhombic cell chosen for simulations. The space group is Cmc21. The I and II assignment is explained in the text. The crystallographic parameters of this system are listed in Table 1.
105
M. Gług et al. / Chemical Physics 459 (2015) 102–111
3. Results and discussion 3.1. Equilibrium and average geometry Prior to dynamical simulations the geometry optimization procedure was applied. Obtained parameters for the ice XI are presented in Table 2. Bond lengths and angles slightly exceed their initial values (about 0.03 A and 2–3° respectively). Analyzing the data from Table 2 one concludes that the cell divides into two sets containing molecules with similar geometry (set I: 1, 2, 5, 6 and set II: 3, 4, 7, 8; see Fig. 1). The differences between sets are caused by various molecular neighborhoods which are responsible for the presence of various hydrogen bonds. There are four different conformations of hydrogen bonds between two water molecules in ordinary ice crystal. These conformations are shown in Fig. 2 where the different H-bond types are tagged a, b, c, and d together with Wong and Whalley notation [55]. Each water molecule participates in creating four hydrogen bonds – two as a proton donor, two as a proton acceptor. The molecules from set I create two d-type hydrogen bonds as an acceptor, one d-type and one a-type as a donor. The molecules from set II create two d-type hydrogen bonds as a donor and one d-type and one a-type as an acceptor. Analysis the symmetry elements of ice XI space group can leads to conclusion that molecules enclosing in both sets are symmetry equivalent. On the base of the nuclei trajectories obtained from Car–Parrinello molecular dynamics simulations the average bond lengths and bond angles were calculated. The results are shown in Tables 3 and 4. Comparing results in Tables 3 and 4 with those presented in Table 2 one can conclude that the equilibrium distances are almost equal to the average. It is also apparent that some of the distances for 4 K are slightly shorter than these at 40 K. In the case of the simulation at 4 K the angles are larger with respect to the optimized geometry of about 0.5 for the O1–H1–H3 group of molecules (set I) and by about 0.2 for the H2–O2–H2 group (set II). Raising the temperature induces even larger changes, i.e. the respective bond angles deviate upwards as 0.8 and 0.2 from average. It can be concluded that while heating the system the two sets of molecules tend to getting similar. We conclude that at higher temperatures the neighborhood is less important for geometry of water molecules. 3.2. Theoretical spectra From the dipole moment autocorrelation function, the product of absorption coefficient and refractive index was obtained, shown in Fig. 3. Using Kramers–Kronig theorem we calculated the dependence of the absorption coefficient on wavenumber, which is shown in Fig. 4. In order to analyze the spectra, especially assignment of the peaks to the particular set of water molecules, the power spectra of atoms belonging to particular water molecules
Table 2 Optimized geometry of water molecules in crystal of ice XI. The labeling according to Fig. 1. Molecule
Set I bond lengths [Å]
1 2 3 4 5 6 7 8
H1–O1 H1–O1 H2–O2 H2–O2 H1–O1 H1–O1 H2–O2 H2–O2
Set II bond lengths [Å] 0.995 0.995 0.999 0.999 0.995 0.995 0.999 0.999
H3–O1 H3–O1 H2–O2 H2–O2 H3–O1 H3–O1 H2–O2 H2–O2
Angles [deg] 0.996 0.996 0.999 0.999 0.996 0.995 0.999 0.999
H1–O1–H3 H1–O1–H3 H2–O2–H2 H2–O2–H2 H1–O1–H3 H1–O1–H3 H2–O2–H2 H2–O2–H2
106.33 106.33 107.75 107.75 106.33 106.33 107.75 107.75
Fig. 2. H-bonded dimer configurations in ice Ih, A (in Wong and Whalley notation: h-cis), B (h-trans), C (c-cis) and D (c-trans).
Table 3 Average geometry of water molecules in crystal of ice XI for simulation performed in 4 K. The numbering according to Fig. 1. Molecule
Set I bond lengths [Å]
1 2 3 4 5 6 7 8
H1–O1 H1–O1 H2–O2 H2–O2 H1–O1 H1–O1 H2–O2 H2–O2
Set II bond lengths [Å] 0.995 0.995 0.999 0.999 0.995 0.995 0.999 0.999
H3–O1 H3–O1 H2–O2 H2–O2 H3–O1 H3–O1 H2–O2 H2–O2
Angles [deg] 0.995 0.995 0.999 0.999 0.995 0.995 0.999 0.999
H1–O1–H3 H1–O1–H3 H2–O2–H2 H2–O2–H2 H1–O1–H3 H1–O1–H3 H2–O2–H2 H2–O2–H2
106.83 106.80 107.93 107.93 106.81 106.81 107.94 107.91
Table 4 Average geometry of water molecules in crystal of ice XI for simulation performed in 40 K. The numbering according to Fig. 1. Molecule
Set I bond lengths [Å]
1 2 3 4 5 6 7 8
H1–O1 H1–O1 H2–O2 H2–O2 H1–O1 H1–O1 H2–O2 H2–O2
Set II bond lengths [Å] 0.995 0.997 1.000 0.999 0.996 0.995 0.999 1.000
H3–O1 H3–O1 H2–O2 H2–O2 H3–O1 H3–O1 H2–O2 H2–O2
Angles [deg] 0.996 0.996 1.000 0.999 0.995 0.996 1.000 0.999
H1–O1–H3 H1–O1–H3 H2–O2–H2 H2–O2–H2 H1–O1–H3 H1–O1–H3 H2–O2–H2 H2–O2–H2
107.10 107.04 107.94 107.94 107.10 107.07 107.96 107.93
were calculated. It should be emphasized that the power spectrum does not always accurately reproduces the form of IR spectrum (especially intensities) because it does not take into account changes of the dipole moment during calculations. Nevertheless, it is assumed that influences of vibrations of atoms belonging to I or II sets on the IR spectra are proportional to intensities presented in the I or II power spectrum. Furthermore, using components of dipole moment along the directions of crystallographic axes, the polarized IR spectra were computed. Comparing Figs. 3 and 4 one can conclude that the representations obtained using Kramers–Kronig relations, presented in Fig. 4, better reproduce experimental spectra, and may be considered as theoretical infrared spectra of ice XI.
106
M. Gług et al. / Chemical Physics 459 (2015) 102–111
Fig. 3. The product of absorption coefficient and refractive index vs. wavenumber for ice XI at 4 K (a), and 40 K (b).
Fig. 4. The absorption coefficient vs. wavenumber of ice XI at 4 K (a) and 40 K (b) (reproduced from Ref. [41]).
Three prominent bands are present in the IR spectra at both temperatures 4 and 40 K: in the range 600–1100 cm1 from librational mode (mL), at 1500–1750 cm1 from the HOH bending mode (mB), and at 2700–3300 from the OH stretching mode (mS). Additionally, at 40 K there appear weak bands in the range 1800–2400 cm1. They are either the overtone of the librational mode (3mL) or bands originating from combined bending and librational mode (mL + mB) [56]. Temperature increase does not significantly affects band broadening. The peak position of mS is shifted
to lower frequency at lower temperature. This phenomenon occurs in the system with strong hydrogen bonds [57]. These outcomes are consistent with experimental studies of infrared spectra of ice XI [22]. In order to study the effect of hydrogen bonding on the form of vibrational spectra we performed simulations for isolated water molecule using the same parameters as for ice. In case of the OH stretching modes we observe pronounced red shifts, whereas for the HOH bending modes very small blue shifts are present (about
M. Gług et al. / Chemical Physics 459 (2015) 102–111
2–3 cm1). Comparing intensities of bending modes we can state that for the HOH bending band, the intensity is significantly greater for single molecule than for the crystal. It is consistent with result of Brubach et al. [58]. They concluded that bending mode band decreases in intensity upon cooling and almost vanishes at the crystallization. Figs. 5 and 6(a) present an enlarged view of the low frequency (librational) region of Fig. 4 compared with power spectra of set I and II. Figs. 5 and 6(b) show comparison of theoretical IR spectra with polarized theoretical spectra in the librational region. The simplest approach to analyze these charts is to determine how each librational mode changes components of the dipole moment DP. In the case that the change is pronounced, it is expected that this vibration should be present in the polarized spectrum along the corresponding direction. Fig. 7 illustrates values and directions of DP for librational modes of water in the crystal structure of ice XI. The twist librations were not taken into account because there is no change of dipole moment in these modes (the axes of these modes are parallel to the dipole moment of water molecules). Fig. 7(a) shows subsystems in crystal structure of ice XI, Fig. 7 (b) three librational modes of single water molecule, whereas Fig. 7(c) shows arbitrary selected positive directions of wagging and rocking liberations for subsystem ABCD. The displacement of modes were obtained using the projection operator for the symmetry operations of Cmc21, following Ref. [27]. Motions of oxygen atoms were ignored because of large oxygen to hydrogen mass ratio. Table 5 presents changes of dipole moment related to proper modes presented in Fig. 7. The content of Table 5 refers to the modes of subsystem A, B, C, D or A0 , B0 , C0 , D0 . It is worth mentioning that the change of dipole moment in b direction is smaller than in c direction. This is because the angle between axis of water molecule and c axis is about 54.7°. The couplings between subsystems were not taken into account. In such case, the following modes are expected to appear in the IR spectrum: W(A)-W⁄(C), W(B)-W(D), W(B)-W⁄(D), R(A)-
107
R(C), R(A)-R⁄(C), R(B)-R(D). ‘‘W” and ‘‘R” refer to ‘‘wagging” and ‘‘rocking” respectively, symbol ‘‘⁄” indicates that direction of motion is negative (see. Fig. 7b), whereas the letter in parentheses represents the molecules, which are in motion. Some modes of subsystem A, B, C, D can couple with modes of subsystem A0 , B0 , C0 , D0 by hydrogen bonds between A and D0 , B and C0 , C and B0 or D and A0 molecules. All of these H-bonds are of a-type (see Fig. 2) with A, B0 , C and D0 as proton acceptor and A0 , B, C0 and D as proton donor. The condition for the occurrence of coupling is the same symmetry of modes. The following modes can couple: W(A)-W(C) with R(B)-R⁄(D), W(A)-W⁄(C) with R(B)-R(D), W(B)-W(D) with R(A)-R⁄(C) and W(B)-W⁄(D) with R(A)-R(C). Table 6 presents the symmetry and behavior of the dipole moment for modes formed by the coupling. In the case of absence of the couplings, the particular peak in the IR spectrum is associated with the power spectrum from only one set (set I and II contain molecules labeled A, C, A0 , C0 and B, D, B0 , D0 respectively). Otherwise, if the coupling is present, the IR peak is associated with peaks from power spectra for two sets. This remark is very useful in following analysis which was performed based on the spectra at 4 K. Comparing power spectra with theoretical IR spectrum (Fig. 5a) we conclude that peaks 1, 2, 4, 5 and 6 represent vibrations of molecules from both set I and set II, whereas peak 3 results from vibrations of molecules from only set I. Additional way to perform the assignments is projecting a set of velocities of particular atoms onto eigenvector corresponding to the particular vibrational mode. For this we have chosen single coordinate for each mode and performed FT of the autocorrelation function calculated for velocity connected with this coordinate. The calculated spectrum can be identified as a contribution of considered mode to the total DOS spectrum. In order to construct normal coordinate we selected first a coordinate describing the motion of each molecule (A, B, C, D) within rocking and wagging librations. In both cases we calculated (for each step of simulation) the sum of unit vectors connected with OAH direction. For rocking
Fig. 5. Comparison of theoretical IR spectrum with power spectra from set I and II (a) and with theoretical polarized IR spectra (b) for ice XI at 4 K. The selected peaks are labeled by numbers for further discussion.
108
M. Gług et al. / Chemical Physics 459 (2015) 102–111
Fig. 6. Comparison of theoretical IR spectrum with power spectra from set I and II (a) and with theoretical polarized IR spectra (b) for ice XI at 40 K.
Fig. 7. (a) Subsystems ABCD and A0 B0 C0 D0 in crystal structure of ice XI, (b) three librational modes of water molecule, (c) positive directions of wagging and rocking modes of water molecules for subsystem ABCD.
Table 5 Symmetry of rocking and wagging librational modes of four molecules subsystems (A, B, C, D or A0 , B0 , C0 , D0 ) with the changes of dipole moment. The symbol ‘‘⁄” indicates that phase is reversed (direction of motion is opposite to that from Fig. 7c). The direction of change of the molecular dipole moment is given in // symbols. The DP means change of the total dipole moment. ‘‘W” and ‘‘R” refer to ‘‘wagging” and ‘‘rocking” respectively. Molecule
A B C D Direction of DP
Wagging
Rocking
A2
B1
B2
A1
A1
B2
B1
A2
W/a/ – W/a/ – 0
W/a/ – W⁄/a/ – a
– W/b + c/ – W/b c/ b
– W/b + c/ – W⁄/b + c/ c
R/b + c/ – R/b + c/ – c
R/b + c/ – R⁄/b c/ – b
– R/a/ – R/a/ a
– R/a/ – R⁄/a/ 0
libration, the resulting vector is projected onto equilibrium plane of molecule. Subsequently, the angle between tangential component and equilibrium sum of unit vectors was calculated. For the wagging case, the procedure was similar, with one
exception – projection was performed onto equilibrium symmetry plane of water molecule. The signs of the angles were consistent with the directions presented in Fig 7. We used OAH unit vectors to ensure that the angles were not varied as an effect of bending
109
M. Gług et al. / Chemical Physics 459 (2015) 102–111
Table 6 Symmetry of coupled rocking/wagging librational modes in ice XI unit cell with the change of dipole moment. The symbol ‘‘⁄” indicates that phase is reversed (direction of motion is opposite to Fig. 7c). The direction of change of the molecular dipole moment is given in // symbols. The DP means change of the total dipole moment. ‘‘W” and ‘‘R” refer to ‘‘wagging” and ‘‘rocking” respectively. Molecule
A B C D Direction of DP
Rocking/wagging A2 in phase
A2 out-of-phase
B1 in phase
B1 out-of-phase
B2 in phase
B2 out-of-phase
A1 in phase
A1 out-of-phase
W/a/ R/a/ W/a/ R⁄/a/ 0
W⁄/a/ R/a/ W⁄/a/ R⁄/a/ 0
W/a/ R/a/ W⁄/a/ R/a/ a
W⁄/a/ R/a/ W/a/ R/a/ a
R/b + c/ W/b + c/ R⁄/bc/ W/bc/ b
R⁄/bc/ W/b + c/ R/b + c/ W/bc/ b
R/b + c/ W/b + c/ R/b + c/ W⁄/b + c/ c
R/b + c/ W⁄/bc/ R/b + c/ W/bc/ c
and stretching modes. The appropriate projections eradicate influence of rocking motion on wagging coordinate and vice versa. The coordinates of simple modes (Table 5) and coupled modes (Table 6) are considered as the sum of angles describing particular librations for given molecules. In case of opposite direction (indicated by ‘‘⁄” symbol) the ‘‘ + ” is changed to ‘‘”. In Figs. 8–11 we
present the power spectra of the modes described above grouped on the base of their symmetry. Peak 1 (at about 620 cm1) overlaps with the peak in the a-polarized spectrum. This indicates that this peak is derived from the mode which changes a-component of the dipole moment. Based on the Tables 5 and 6 we conclude that the origin of this band is in the mode W(A)-W⁄(C) or R(B)-R(D) or coupled mode W(A)-R(B)-W⁄ (C)-R(D) or W⁄(A)-R(B)-W(C)-R(D). As previously stated (using the power spectra), this peak is the result of coupling of vibrations of molecules from both sets. For this reason peak 1 can be assigned to one of the coupled modes. In the case of peak 5 (945 cm1), the conclusion is analogous. Assigning W(A)-R(B)-W⁄(C)-R(D) mode to peak 5 is consistent with experimental Raman spectroscopic study [27] and model calculations [59]. Then the W⁄(A)-R(B)-W(C)-R(D) mode can be assigned to peak 1. Fig. 8 confirms assignments for peaks 1 and 5. Spectra originating from W(A)-W⁄(C) and R(B)-R(D) contain two shared peaks – at 620 and 945 cm1. This fact indicates that these modes can be coupled. In order to determine properly in phase/out-of-phase mode coupling, two top spectra were drawn. Based on positions of main peaks, we conclude that the peak at 620 cm1 originates from the W⁄(A)-R(B)-W(C)-R(D) mode, whereas the peak at 945 cm1 is connected with the W(A)-R(B)-W⁄(C)-R(D) mode. Additionally, there is a peak at 860 cm1 in the W(A)-W⁄(C) spectrum and at 760 cm1 in the R(B)-R(D) spectrum. These peaks have not corresponding equivalent in another spectrum. It is highly probable that these peaks are the remains of simple modes W(A)-W⁄(C) and R(B)-R(D). This conjecture complies with the fact that coupling of two vibrations leads to one mode characterized by lower frequency and the other by higher frequency.
Fig. 9. DOS functions for B2 modes. Black line – R(A)-R*(C), red line – W(B)-W(D), green line – R(A)-W(B)-R*(C)-W(D), blue line – R*(A)-W(B)-R(C)-W(D). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 10. DOS functions for A1 modes. Black – R(A)-R(C), red – W(B)-W*(D), green – R (A)-W(B)-R(C)-W*(D), blue – R(A)-W*(B)-R(C)-W(D). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 8. DOS functions for B1 modes. Black – W(A)-W*(C), red – R(B)-R(D), green – W (A)-R(B)-W*(C)-R(D), blue – W*(A)-R(B)-W(C)-R(D). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
110
M. Gług et al. / Chemical Physics 459 (2015) 102–111
Fig. 11. DOS functions for A2 modes. Black – W(A)-W(C), red – R(B)-R*(D), green – W(A)-R(B)-W(C)-R*(D), blue – W*(A)-R(B)-W*(C)-R*(D). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Peak 3 (at about 710 cm1) is associated with the change of b-component of the dipole moment. It may result from R(A)R⁄(C) or W(B)-W(D) or coupled R(A)-W(B)-R⁄(C)-W(D) or R⁄(A)W(B)-R(C)-W(D) vibrations. Based on power spectra, one can conclude, this peak originates from set I. This facts leads to conclusion that peak 3 can be assigned to the R(A)-R⁄(C) mode. Peak 6 (at 1015 cm1) is related to the change of the dipole moment along b direction. Power spectra of both sets indicate that the origin of this peak can be a coupled mode. Based on the analysis of Fig. 9, similar to analysis of Fig. 8, we can state that peak 6 is associated with the coupled mode R⁄(A)-W(B)-R(C)-W(D) and confirm that peak 3 is related to R(A)-R⁄(C). Furthermore, the R(A)-W(B)-R⁄(C)-W(D) mode is localized at about 620 cm1. This means that peak 1 originates from two modes W⁄(A)-R(B)-W(C)-R(D) and R(A)-W(B)R⁄(C)-W(D). Peak at about 900 cm1 corresponds to the W(B)-W (D) mode (as for peak 3, it is a remain of the coupling). This peak could be present in IR spectra because of change of the dipole moment for the W(B)-W(D) mode. However, the intensity is insufficient to distinguish it. Peaks 2 and 4 (at about 700 and 925 cm1) are related to the change of the dipole moment in c-direction. It may result from W(B)-W⁄(D), R(A)-R(C), R(A)-W(B)-R(C)-W⁄(D) or R(A)-W⁄(B)-R (C)-W(D) vibrations. Based on power spectra, they origin from set I and II. This facts leads to the conclusion that these peaks are associated with coupled mode. Based on Fig. 10, we can conclude that peak 2 corresponds to R(A)-W⁄(B)-R(C)-W(D) and peak 4 to R(A)-W(B)-R(C)-W⁄(D). In addition to previously discussed bands, the peaks at 800–875 cm1 are present in the power spectra. They have no counterparts in the IR spectra. These peaks can arise from the modes which do not change the dipole moment, i.e. W(A)-W(C), R(B)-R⁄(D) and the coupled modes: W(A)-R(B)-W(C)-R⁄(D) or W⁄(A)-R(B)-W⁄(C)-R⁄(D) (symmetry A2). This conclusion is partly confirmed by spectra in Fig. 11. Another possible origin of these peak can be twist modes. For spectra calculated at 40 K the analysis is more difficult because of greater noises. However the results of analysis are very similar. The main difference is slight shift to lower wavenumbers. In addition, the peak at 1100 cm1 is present in the IR spectrum, in polarized spectrum (along b axis) and in power spectra for both sets. This can be the result of combination of a librational and a translational modes [60,27].
Fig. 12. Comparison of the DOS function obtained from CPMD calculations, IINS spectrum (taken from Ref. [24]), DOS function obtained by model calculation (taken from Ref. [59]), and DOS function obtained from classical MD calculations (taken from Ref. [61]) in the librational region of H2O in the ice XI phase.
3.3. Comparison with other theoretical and experimental data In order to verify the accuracy of our investigation, the entire power spectrum of ice XI at 4 K (which is the sum of all atomic power spectra) was compared, in Fig. 12, with density of states (DOS) observed in the inelastic incoherent neutron scattering (IINS) of H2O in ice XI [24], with DOS obtained from model calculations [59] and with DOS from classical molecular dynamics simulation [61]. In order to facilitate the comparison, our calculated CPMD DOS was presented as the smooth bandshape obtained by convoluting the CPMD resulting data with gaussian functions. There is very good agreement in the positions and relative intensities of our CPMD calculated DOS and IINS spectrum, except the peak at about 800 cm1, which it is the most intensive one in the CPMD DOS. The agreement between DOS obtained by Iwano et al. [59] with our CPMD DOS, as well as with IINS DOS, is good in the higher-energy part, whereas the peak at about 600 cm1 seems to be doubly splitted in the DOS obtained from model calculations. The DOS obtained by Itoh et al. [61] using classical molecular dynamics calculations reveal some discrepancies in peak frequencies, compared with three DOS functions discussed above, however qualitative similarities are obviously seen. The relative peak positions are reproduced well and the entire DOS function is shifted towards higher wavenumbers. The difference in the peak frequencies was explained by larger elastic constants in the classical MD calculations. Analyzing Fig. 12, we conclude that our power spectrum reproduces most features of the experimental IINS spectrum, including both peak positions and relative intensities, and is also in qualitative agreement with other calculated DOS functions. Based on measurement of infrared spectra of KOH doped water at 4–285 K, Arakawa et. al. [22] suggested that ordered ice in space maybe detectable using infrared telescopes. They stated, there are two peaks in librational region of infrared spectrum of ice XI – major at 850 cm1 and minor at 590 cm1 (in the IR spectrum of ice Ih only one broad peak is observed in this region). In order to verify the above thesis we performed CPMD calculations of the IR spectra of ice XI at 4 and 40 K. These spectra show two major peaks in the librational region, in agreement with Arakawa predictions, with similar intensities similar to the experimental ones.
M. Gług et al. / Chemical Physics 459 (2015) 102–111
4. Conclusions Theoretical infrared spectra and DOS spectra of ordered phase ferroelectric ice XI have been calculated at 4 and 40 K using Car–Parrinello molecular dynamics and Kramers–Kronig relations. The detailed analysis of vibrations in the librational region (at 550–1150 cm1) was performed. The peaks present in the IR spectra were assigned to particular modes based on power spectra and polarized IR spectra. Our study shows that not all peaks result from the coupling between librational modes. Comparison presented in Fig. 8 shows that our results are in good agreement with available experimental and theoretical spectra of ice XI. Conflict of interest There is no conflict of interest. Acknowledgments This research was supported in part by PL-Grid Infrastructure. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18]
S. Kawada, J. Phys. Soc. Jpn. 32 (1972) 1442. Y. Tajima, T. Matsuo, H. Suga, J. Phys. Chem. Solids 45 (1984) 1135. T. Matsuo, Y. Tajima, H. Suga, J. Phys. Chem. Solids 47 (1986) 165. Y. Tajima, T. Matsuo, H. Suga, J. Chem. Phys. 82 (1985) 424. M. Jackson, W. Whithworth, J. Chem. Phys. 103 (1995) 7647. M.E. Brown, W.M. Calvin, Science 287 (2000) 107. M.J. Iedema, M.J. Dresser, D.L. Doering, J.B. Rowland, W.P. Hess, A.A. Tsekouras, J.P. Cowin, J. Phys. Chem. B 102 (1998) 9203. H. Wang, R.C. Bell, M.J. Iedema, A.A. Tsekouras, J.P. Cowin, Astrophys. J. 620 (2005) 1027. H.X. Zhao, X.J. Kong, H. Li, Y.C. Jin, L.S. Long, X.C. Zeng, R.B. Huang, L.S. Zheng, Proc. Natl. Acad. Sci. 108 (2011) 3481. H. Fukazawa, A. Hoshikawa, Y. Ishii, B.C. Chakoumakos, J.A. Fernandez-Baca, ApJ 652 (2006) L57. H. Fukazawa, A. Hoshikawa, B.C. Chakoumakos, J.A. Fernandez-Baca, Nucl. Instrum. Methods Phys. Res. A 600 (2009) 279. X.C. Su, L. Lianos, Y.R. Shen, G.A. Somorjai, Phys. Rev. Lett. 80 (1998) 1533. H. Wang, R.C. Bell, M.J. Iedema, A.A. Tsekouras, J.P. Cowin, ApJ 620 (2005) 1027. A.H. Castro Neto, P. Pujol, E. Fradkin, Phys. Rev. B 74 (2006) 024302. X. Fan, D. Bing, J. Zhang, Z. Shen, J.-L. Kuo, Comput. Mater. Sci. 49 (2010) S170. H. Fukazawa, S. Mae, S. Ikeda, O. Watanabe, Chem. Phys. Lett. 294 (1998) 554. S. Kawada, J. Phys. Soc. Jpn. 58 (1989) 295. H. Fukazawa, S. Ikeda, S. Mae, Chem. Phys. Lett. 282 (1998) 215.
111
[19] H. Suga, Thermochim. Acta 300 (1997) 117. [20] K. Furic´, V. Volovšek, J. Mol. Struct. 976 (2010) 174. [21] H. Fukazawa, S. Ikeda, M. Oguro, S.M. Bennington, S. Mae, J. Chem. Phys. 118 (2003) 1577. [22] M. Arakawa, H. Kagi, H. Fukazawa, Astrophys. J. Suppl. 184 (2009) 361. [23] J.C. Li, V.M. Nield, S.M. Jackson, Chem. Phys. Lett. 241 (1995) 290. [24] J.C. Li, J. Chem. Phys. 105 (1996) 6733. [25] K. Abe, T. Miasa, Y. Ohtake, K. Nakano, M. Nakajima, H. Yamamoto, T. Shigenari, J. Korean Phys. Soc. 46 (2005) 300. [26] K. Abe, T. Shigenari, J. Chem. Phys. 134 (2011) 104506. [27] T. Shigenari, K. Abe, J. Chem. Phys. 136 (2012) 174504. [28] M.J. Wójcik, M. Gług, M. Boczar, Ł. Boda, Chem. Phys. Lett. 612 (2014) 162. [29] R. Car, M. Parrinello, Phys. Rev. Lett. 55 (1985) 2471. [30] M.S. Bergren, S.A. Rice, J. Chem. Phys. 77 (1982) 583. [31] S.A. Rice, M.S. Bergren, A.C. Belch, G. Nielson, J. Phys. Chem. 87 (1983) 4295. [32] M.J. Wójcik, V. Buch, J.P. Devlin, J. Chem. Phys. 99 (1993) 2332. [33] M.J. Wójcik, K. Szczeponek, S. Ikeda, J. Chem. Phys. 117 (2002) 9850. [34] F. Li, J.L. Skinner, J. Chem. Phys. 132 (2010) 204505. [35] F. Li, J.L. Skinner, J. Chem. Phys. 133 (2010) 244504. [36] M. Benoit, M. Bernasconi, P. Focher, M. Parrinello, Phys. Rev. Lett. 76 (1996) 2934. [37] M. Bernasconi, P.L. Silvestrelli, M. Parrinello, Phys. Rev. Lett. 81 (1998) 1235. [38] C. Knight, S.J. Singer, J.-L. Kuo, T.K. Hirsch, L. Ojamäe, M.L. Klein, Phys. Rev. E 73 (2006) 056113. [39] M. Schönherr, B. Slater, J. Hutter, J.V. Vondele, J. Phys. Chem. B 118 (2014) 590. [40] M. Brela, J. Stare, G. Pirc, M. Sollner-Dolenc, M. Boczar, M.J. Wójcik, J. Mavri, J. Phys. Chem. B 116 (2012) 4510. [41] J. Stare, J. Panek, J. Eckert, J. Grdadolnik, J. Mavri, D. Hadzi, J. Phys. Chem. A 112 (2008) 1576. [42] W. Kohn, J. Sham, Phys. Rev. 140 (1965) A1133. [43] F.A. Bornemann, C. Schütte, Numer. Math. 78 (1998) 359. [44] R. Ramírez, T. López-Ciudad, P.P. Kumar, D. Marx, J. Chem. Phys. 121 (2004) 3973. [45] N. Marzari, D. Vanderbilt, Phys. Rev. B 56 (1997) 12847. [46] R. de L. Kronig, J. Opt. Soc. Am. 12 (1926) 547. [47] H.A. Kramers, Atti Cong. Intern. Fisici, (Transactions of Volta Centenary Congress) Como 2, 1927, 545. [48] F. Wooten, Optical Properties of Solids, Academic Press, New York, 1972. [49] T.K. Hirsch, L. Ojamäe, J. Phys. Chem. B 108 (2004) 15856. [50] CPMD,
Contents lists available at ScienceDirect
Chemical Physics journal homepage: www.elsevier.com/locate/chemphys
Analysis of librational modes of ice XI studied by Car–Parrinello molecular dynamics Maciej Gług, Marek Boczar ⇑, Łukasz Boda, Marek J. Wójcik ⇑ Faculty of Chemistry, Jagiellonian University, 30-060 Kraków, Ingardena 3, Poland
a r t i c l e
i n f o
Article history: Received 22 May 2015 In final form 9 August 2015 Available online 14 August 2015 Keywords: Ice XI IR spectra Car–Parrinello molecular dynamics
a b s t r a c t Theoretical infrared spectra of ice XI at 4 and 40 K were obtained by Car–Parrinello molecular dynamics calculations and evaluating the dipole moment using maximally-localized Wannier functions. In order to improve the agreement between theoretical and experimental spectra, the Kramers–Kronig relations were applied. Based on analysis of power and polarization spectra the assignment of modes was carried out. The results were compared with experimental and theoretical spectra of ice XI, performed by other authors. Such simulations can be very useful while investigating the existence of ice XI in space. Ó 2015 Elsevier B.V. All rights reserved.
1. Introduction Ice XI, discovered in 1972 by Kawada [1] and Tajima et al. [2,3], is proton-ordered, low-temperature phase of ice Ih (ordinary hexagonal ice). Few years later it was found, using neutron diffraction, that ice XI belongs to the orthorhombic Cmc21 symmetry space group with eight water molecules per unit cell [4]. The ferroelectric character of this structure was subsequently confirmed by thermal stimulated depolarization experiments [5]. Ice XI and its ferroelectricity have attracted scientific interest mainly after discovering the presence of crystalline ice in outer solar system [6]. Long-range electrostatic forces, caused by the ferroelectricity, might be an important factor for planet formation [7–9]. Fukazawa et al. [10,11] presented hypothesis that ice XI can exist in interstellar space. Electrostatic force caused by the ferroelectric properties might affect agglomeration of ice particles in the space [7,12,13]. Such presence of ice XI might be very significant for planetary formation processes. Infrared telescope measurements do not provide clear information whether ice XI exists in space or not. This is mainly due to lack of reference data. There is a strong need to obtain good-quality IR spectra of ice XI, either experimental or theoretical. Due to the high absorptivities and difficulties in preparation of pure samples it is difficult to obtain experimental infrared spectra of ice XI. For this reason, theoretical study of IR spectra of ice XI, presented in this paper, is very important to provide reference spectra for astrophysical measurements. ⇑ Corresponding authors. Fax: +48 12 634 05 15. E-mail addresses: [email protected] (M. Boczar), [email protected] (M.J. Wójcik). http://dx.doi.org/10.1016/j.chemphys.2015.08.004 0301-0104/Ó 2015 Elsevier B.V. All rights reserved.
Ice XI is the low temperature equilibrium structure of hexagonal ice (ice Ih) – it has a triple point with hexagonal ice and gaseous water at 72 K and 0 Pa. The arrangement of water molecules in hexagonal ice, in which they are surrounded by four semi-randomly directed hydrogen bonds, should change to a regular arrangement of hydrogen bonds at low temperatures. This process becomes more effective with an increase in pressure [14]. In nature the transformation is very slow without the presence of catalysts, even though ice XI is thought to have a more stable conformation than ice Ih, mainly due to relatively small energy difference between the two structures [15]. Fukazawa et al. [16] claimed the presence of ice XI in the Antarctic ice. They found that the phase transformation of ice XI to ice Ih occurs at temperature 237 K, which is higher than the temperature of the expected triple point (72 K, 0 Pa). Ice XI can be obtained under laboratory conditions from a dilute (10 mM) KOH solution kept just below 72 K at ambient pressure for about a week [17,18]. The mechanism of the influence of hydroxide and K+ ions on this process was discussed by Suga [19], however many aspects are still not exactly understood. The ice Ih M XI transition was previously reported only for the doped, not for pure water, however ice XI was also found experimentally using pure water at very low temperature (10 K) at low pressure [20]. Such conditions are thought to be present in the upper atmosphere. Ice XI was subject of numerous spectroscopic studies. The vibrational spectra, which reveal the crystal structure and atomic motions, are investigated using mainly infrared, Raman and neutron spectroscopy. In the ice crystal water molecules are connected by a network of hydrogen bonds. In the crystal there are
103
M. Gług et al. / Chemical Physics 459 (2015) 102–111
intermolecular translational lattice and librational vibrations and intramolecular OAH stretching vibrations. Fukazawa et al. [18,21] have postulated that the largest difference in the spectroscopic data between ice XI and ice Ih occurs in the region of librational modes. Infrared absorption spectra of KOH-doped ice in the librational region for different temperatures was measured by Arakawa et al. [22]. Inelastic incoherent neutron scattering of ice XI was reported by Li et al. [23,24], who also studied the ice Ih M ice XI transition using Raman spectroscopy [23]. Phase transition and the differences of the vibrational modes between ice Ih and XI in the translational and librational region was also studied by Raman spectroscopy by Abe et al. [25]. Recently experimental studies of the Raman spectra of hydrogenordered ice have been reported by Abe and Shigenari [26,27]. In our recent paper [28] we presented comparison between theoretical IR spectra of ice Ih and ice XI obtained by Car–Parrinello molecular dynamics (CPMD) calculations [29]. In this article we present the details of the crystal structure of ice XI and analysis of power and polarization spectra of ice XI. Different forms of ice have been studied in the literature. Infrared and Raman spectra of ice Ih, cubic and amorphous ice have been previously calculated using different potentials and theoretical approaches by Rice et al. [30,31], Wójcik et al. [32,33] and Skinner et al. [34,35]. CPMD calculations of vibrational frequencies, phase transitions and dielectric properties in different phases of ice were reported by Benoit et al. [36], Bernasconi et al. [37], Knight et al. [38] and Schönherr et al. [39]. In our study spectra are calculated from the velocity/dipole moment autocorrelation function (it is described in Section 2.2). In this way they are not state resolved. One of approaches for quantization of nuclear motion was used by Brela et al. [40] and Stare et al. [41]. These authors quantize selected nuclear degrees of freedom. This takes into account the influence of the quantum nature of nuclear motion on the resultant spectroscopic properties. We plan to use similar method in the future.
Born–Oppenheimer molecular dynamics, but requires shorter time step [43]. Hence, the choice of the l value should be the compromise between the accuracy and the time of the computation. 2.2. Theoretical infrared spectrum IR absorption coefficient a(x) is commonly calculated using the formula [44]
NðxÞ ¼ nðxÞ þ ikðxÞ;
All of the calculations reported in this paper were carried out within the Car–Parrinello ab initio molecular dynamics method [29]. In this approach the electronic degrees of freedom as (fictitious) dynamical variables are introduced using Car–Parrinello Lagrangian (embedded in Kohn–Sham formalism [42]):
X 1X _2 X _ _ ¼ M I RI þ lh/i j/i i hW0 jHKS Kij ðh/i j/j i dij Þ: e jW0 i þ 2 I i i;j ð1Þ
This leads to the system of coupled equations of motion for both, ions and electrons:
€ I ¼ rI hW0 jHKS jW0 i; MI R e
ð2aÞ
X Kij /j :
ð2bÞ
l/€ i ¼ HKS e /i þ
ð4Þ
where the real part of the refractive index n indicates the phase velocity, whereas the imaginary part k is responsible for the amount of absorption loss, when the electromagnetic wave propagates through the material, and k is related to the absorption coefficient by the relation
aðxÞc : x
ð5Þ
Electric susceptibility v is connected with N index by the formula
2.1. Car–Parrinello molecular dynamics
LCP
ð3Þ
where T is temperature, V is volume of cell, b is 1/(kT), c is speed of light in vacuum, n(x) is refractive index and M is total dipole moment. The brackets indicate autocorrelation function. M(t) function is the sum of electronic and ionic dipole moments computed in every step of the simulation. The ionic contribution is calculated straightforward whereas electrons dipole moment is obtained using maximally-localized Wannier functions [45]. Due to the fact, that unknown frequency dependence of n cannot be neglected, this equation provides only the product of n(x) and a(x). However, it is possible to extract absorption coefficient a(x) by using Kramers–Kronig theorem [46,47] giving the relations between the real and imaginary part of the dielectric function. The relations are of a general nature and are based on the properties of a complex, analytical response function, such as the complex refractive index in our case, defined as
kðxÞ ¼ 2. Computational details
Z 1 2bpx2 expðixtÞhMð0ÞMðtÞidt; 3VcnðxÞ 1
aðxÞ ¼
N2 ðxÞ ¼ 1 þ 4pvðxÞ
Comparing Eqs. (4) and (6) we obtain the following expression for v 2
vðxÞ ¼ v1 ðxÞ þ iv2 ðxÞ ¼
Consequently, it is not necessary to solve electronic Schrödinger Equation at each time step of simulation (opposite to Born–Oppenheimer molecular dynamics). In Eq. (1) the first component is the kinetic energy of the ions, followed by the kinetic energy of the electrons, and the finally potential energy and orbitals orthonormality constrains, respectively. The parameter l is called the fictitious mass of the electron. Reducing of the value of this parameter yields an increased accuracy of electronic energy calculations, with respect to
n2 ðxÞ k ðxÞ 1 nðxÞkðxÞ þ i 4p 2p
ð7Þ
Both real and imaginary parts of electric susceptibility are mutually connected by Kramers–Kronig relations [48]. This leads to the equation 2
v1 ðxÞ ¼
n2 ðxÞ k ðxÞ 1 2 ¼ Pr 4p p
Z
1
0
x0 nðx0 Þkðx0 Þ=2p 0 dx ; x02 x2
ð8Þ
which can be represented using absorption coefficient 2
j
ð6Þ
v1 ðxÞ ¼
2
n2 ðxÞ a ðxx2Þc 1 c ¼ 2 Pr p 4p
Z 0
1
nðx0 Þaðx0 Þ dx0 ; x02 x2
ð9Þ
where Pr denotes the Cauchy principal value of the integral, which is the limit for d ? 0 of the sum of the integrals over 1 < x0 < x d and x + d < x0 < 1. Basing on above formula it can be concluded that knowing the product of absorption coefficient and refractive index for all frequencies allows to calculate real part of susceptibility. Denoting this product as I(x) the system of two equations can be written 2
2
n2 ðxÞ a ðxx2Þc 1 ¼ v1 ðxÞ 4p
ð10Þ
104
M. Gług et al. / Chemical Physics 459 (2015) 102–111
nðxÞaðxÞ ¼ IðxÞ
ð11Þ
The right side of each equation is given by Eqs. (9) and (3) respectively. The solution with respect to n(x) and a(x) is
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u 2 2 u1 þ 4pv1 ðxÞ ð1 þ 4pv1 ðxÞÞ þ 4I2 ðxÞ xc 2 aðxÞ ¼ t 2 2 xc 2
ð12Þ
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi11 0v qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u 2 1 þ 4pv1 ðxÞ ð1 þ 4pv1 ðxÞÞ2 þ 4I2 ðxÞ xc 2 C Bu t C nðxÞ ¼ IðxÞB 2 A @ 2 xc 2 ð13Þ The finally formula on absorption coefficient is
aðxÞ ¼
x c
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 Z 1 u 1 2c Z 1 Iðx0 Þ 1 4c IðxÞ c2 t 0þ 0 d x P d x þ 4I2 ðxÞ 2 : 1 þ P 2 2 p 0 x02x 2 p 0 x02 x2 x
thermostat [51]. To calculate the cell dipole moment, maximally localized Wannier functions [45] were used. A fictitious mass of 300 a.u. was set for the electronic degrees of freedom. The equations of motion for nuclei and molecular orbitals were solved using the velocity Verlet algorithm with 1.5 a.u. timestep. Electrons for oxygen and hydrogen atoms were treated within the DFT framework with BLYP functional and the norm-conserving pseudopotential of Goedecker type [52–54]. The electronic wave functions were expanded in plane waves. The energy cutoff was adjusted to 230 Ry. Brillouin zone sampling was restricted to the Gamma point. The total length of simulation path after equilibrating was approx. 12 ps (330,000 steps). Periodic boundary conditions were employed for the solid state calculations with the system which theoretically extends to infinity. Power spectrum and atomic power spectra were calculated using the Fourier transformation (FT) of the atomic velocity autocorrelation function. IR spectrum was obtained from FT of the dipole moment autocorrelation function and the application of Kramers–Kronig relations. All FT calculations were performed with CPMD additional script (FOURIER.X).
ð14Þ 2.3. Simulation details Orthorhombic cell, shown in Fig. 1, containing eight water molecules with the parameters [49] included in Table 1 was selected for simulation. All calculations were executed using CPMD package [50]. Car–Parrinello molecular dynamics simulation was performed in thermodynamic equilibrium using Nosé–Hoover chain
Table 1 Positions of atoms for system from Fig. 1. The length of all OAH bond is 0.9572 Å and the value of all angles HAOAH is 104.5°. The cell is orthorhombic with a = 4.49225 Å, b = 7.78080 Å, c = 7.33581 Å [59]. Atom
Fractional coordinates
O1 O2 H1 H2 H3
0.0000 0.5000 0.0000 0.6685 0.0000
0.3333 0.1667 0.4510 0.2282 0.3386
0.5625 0.4375 0.5244 0.4836 0.6929
Fig. 1. Ice XI orthorhombic cell chosen for simulations. The space group is Cmc21. The I and II assignment is explained in the text. The crystallographic parameters of this system are listed in Table 1.
105
M. Gług et al. / Chemical Physics 459 (2015) 102–111
3. Results and discussion 3.1. Equilibrium and average geometry Prior to dynamical simulations the geometry optimization procedure was applied. Obtained parameters for the ice XI are presented in Table 2. Bond lengths and angles slightly exceed their initial values (about 0.03 A and 2–3° respectively). Analyzing the data from Table 2 one concludes that the cell divides into two sets containing molecules with similar geometry (set I: 1, 2, 5, 6 and set II: 3, 4, 7, 8; see Fig. 1). The differences between sets are caused by various molecular neighborhoods which are responsible for the presence of various hydrogen bonds. There are four different conformations of hydrogen bonds between two water molecules in ordinary ice crystal. These conformations are shown in Fig. 2 where the different H-bond types are tagged a, b, c, and d together with Wong and Whalley notation [55]. Each water molecule participates in creating four hydrogen bonds – two as a proton donor, two as a proton acceptor. The molecules from set I create two d-type hydrogen bonds as an acceptor, one d-type and one a-type as a donor. The molecules from set II create two d-type hydrogen bonds as a donor and one d-type and one a-type as an acceptor. Analysis the symmetry elements of ice XI space group can leads to conclusion that molecules enclosing in both sets are symmetry equivalent. On the base of the nuclei trajectories obtained from Car–Parrinello molecular dynamics simulations the average bond lengths and bond angles were calculated. The results are shown in Tables 3 and 4. Comparing results in Tables 3 and 4 with those presented in Table 2 one can conclude that the equilibrium distances are almost equal to the average. It is also apparent that some of the distances for 4 K are slightly shorter than these at 40 K. In the case of the simulation at 4 K the angles are larger with respect to the optimized geometry of about 0.5 for the O1–H1–H3 group of molecules (set I) and by about 0.2 for the H2–O2–H2 group (set II). Raising the temperature induces even larger changes, i.e. the respective bond angles deviate upwards as 0.8 and 0.2 from average. It can be concluded that while heating the system the two sets of molecules tend to getting similar. We conclude that at higher temperatures the neighborhood is less important for geometry of water molecules. 3.2. Theoretical spectra From the dipole moment autocorrelation function, the product of absorption coefficient and refractive index was obtained, shown in Fig. 3. Using Kramers–Kronig theorem we calculated the dependence of the absorption coefficient on wavenumber, which is shown in Fig. 4. In order to analyze the spectra, especially assignment of the peaks to the particular set of water molecules, the power spectra of atoms belonging to particular water molecules
Table 2 Optimized geometry of water molecules in crystal of ice XI. The labeling according to Fig. 1. Molecule
Set I bond lengths [Å]
1 2 3 4 5 6 7 8
H1–O1 H1–O1 H2–O2 H2–O2 H1–O1 H1–O1 H2–O2 H2–O2
Set II bond lengths [Å] 0.995 0.995 0.999 0.999 0.995 0.995 0.999 0.999
H3–O1 H3–O1 H2–O2 H2–O2 H3–O1 H3–O1 H2–O2 H2–O2
Angles [deg] 0.996 0.996 0.999 0.999 0.996 0.995 0.999 0.999
H1–O1–H3 H1–O1–H3 H2–O2–H2 H2–O2–H2 H1–O1–H3 H1–O1–H3 H2–O2–H2 H2–O2–H2
106.33 106.33 107.75 107.75 106.33 106.33 107.75 107.75
Fig. 2. H-bonded dimer configurations in ice Ih, A (in Wong and Whalley notation: h-cis), B (h-trans), C (c-cis) and D (c-trans).
Table 3 Average geometry of water molecules in crystal of ice XI for simulation performed in 4 K. The numbering according to Fig. 1. Molecule
Set I bond lengths [Å]
1 2 3 4 5 6 7 8
H1–O1 H1–O1 H2–O2 H2–O2 H1–O1 H1–O1 H2–O2 H2–O2
Set II bond lengths [Å] 0.995 0.995 0.999 0.999 0.995 0.995 0.999 0.999
H3–O1 H3–O1 H2–O2 H2–O2 H3–O1 H3–O1 H2–O2 H2–O2
Angles [deg] 0.995 0.995 0.999 0.999 0.995 0.995 0.999 0.999
H1–O1–H3 H1–O1–H3 H2–O2–H2 H2–O2–H2 H1–O1–H3 H1–O1–H3 H2–O2–H2 H2–O2–H2
106.83 106.80 107.93 107.93 106.81 106.81 107.94 107.91
Table 4 Average geometry of water molecules in crystal of ice XI for simulation performed in 40 K. The numbering according to Fig. 1. Molecule
Set I bond lengths [Å]
1 2 3 4 5 6 7 8
H1–O1 H1–O1 H2–O2 H2–O2 H1–O1 H1–O1 H2–O2 H2–O2
Set II bond lengths [Å] 0.995 0.997 1.000 0.999 0.996 0.995 0.999 1.000
H3–O1 H3–O1 H2–O2 H2–O2 H3–O1 H3–O1 H2–O2 H2–O2
Angles [deg] 0.996 0.996 1.000 0.999 0.995 0.996 1.000 0.999
H1–O1–H3 H1–O1–H3 H2–O2–H2 H2–O2–H2 H1–O1–H3 H1–O1–H3 H2–O2–H2 H2–O2–H2
107.10 107.04 107.94 107.94 107.10 107.07 107.96 107.93
were calculated. It should be emphasized that the power spectrum does not always accurately reproduces the form of IR spectrum (especially intensities) because it does not take into account changes of the dipole moment during calculations. Nevertheless, it is assumed that influences of vibrations of atoms belonging to I or II sets on the IR spectra are proportional to intensities presented in the I or II power spectrum. Furthermore, using components of dipole moment along the directions of crystallographic axes, the polarized IR spectra were computed. Comparing Figs. 3 and 4 one can conclude that the representations obtained using Kramers–Kronig relations, presented in Fig. 4, better reproduce experimental spectra, and may be considered as theoretical infrared spectra of ice XI.
106
M. Gług et al. / Chemical Physics 459 (2015) 102–111
Fig. 3. The product of absorption coefficient and refractive index vs. wavenumber for ice XI at 4 K (a), and 40 K (b).
Fig. 4. The absorption coefficient vs. wavenumber of ice XI at 4 K (a) and 40 K (b) (reproduced from Ref. [41]).
Three prominent bands are present in the IR spectra at both temperatures 4 and 40 K: in the range 600–1100 cm1 from librational mode (mL), at 1500–1750 cm1 from the HOH bending mode (mB), and at 2700–3300 from the OH stretching mode (mS). Additionally, at 40 K there appear weak bands in the range 1800–2400 cm1. They are either the overtone of the librational mode (3mL) or bands originating from combined bending and librational mode (mL + mB) [56]. Temperature increase does not significantly affects band broadening. The peak position of mS is shifted
to lower frequency at lower temperature. This phenomenon occurs in the system with strong hydrogen bonds [57]. These outcomes are consistent with experimental studies of infrared spectra of ice XI [22]. In order to study the effect of hydrogen bonding on the form of vibrational spectra we performed simulations for isolated water molecule using the same parameters as for ice. In case of the OH stretching modes we observe pronounced red shifts, whereas for the HOH bending modes very small blue shifts are present (about
M. Gług et al. / Chemical Physics 459 (2015) 102–111
2–3 cm1). Comparing intensities of bending modes we can state that for the HOH bending band, the intensity is significantly greater for single molecule than for the crystal. It is consistent with result of Brubach et al. [58]. They concluded that bending mode band decreases in intensity upon cooling and almost vanishes at the crystallization. Figs. 5 and 6(a) present an enlarged view of the low frequency (librational) region of Fig. 4 compared with power spectra of set I and II. Figs. 5 and 6(b) show comparison of theoretical IR spectra with polarized theoretical spectra in the librational region. The simplest approach to analyze these charts is to determine how each librational mode changes components of the dipole moment DP. In the case that the change is pronounced, it is expected that this vibration should be present in the polarized spectrum along the corresponding direction. Fig. 7 illustrates values and directions of DP for librational modes of water in the crystal structure of ice XI. The twist librations were not taken into account because there is no change of dipole moment in these modes (the axes of these modes are parallel to the dipole moment of water molecules). Fig. 7(a) shows subsystems in crystal structure of ice XI, Fig. 7 (b) three librational modes of single water molecule, whereas Fig. 7(c) shows arbitrary selected positive directions of wagging and rocking liberations for subsystem ABCD. The displacement of modes were obtained using the projection operator for the symmetry operations of Cmc21, following Ref. [27]. Motions of oxygen atoms were ignored because of large oxygen to hydrogen mass ratio. Table 5 presents changes of dipole moment related to proper modes presented in Fig. 7. The content of Table 5 refers to the modes of subsystem A, B, C, D or A0 , B0 , C0 , D0 . It is worth mentioning that the change of dipole moment in b direction is smaller than in c direction. This is because the angle between axis of water molecule and c axis is about 54.7°. The couplings between subsystems were not taken into account. In such case, the following modes are expected to appear in the IR spectrum: W(A)-W⁄(C), W(B)-W(D), W(B)-W⁄(D), R(A)-
107
R(C), R(A)-R⁄(C), R(B)-R(D). ‘‘W” and ‘‘R” refer to ‘‘wagging” and ‘‘rocking” respectively, symbol ‘‘⁄” indicates that direction of motion is negative (see. Fig. 7b), whereas the letter in parentheses represents the molecules, which are in motion. Some modes of subsystem A, B, C, D can couple with modes of subsystem A0 , B0 , C0 , D0 by hydrogen bonds between A and D0 , B and C0 , C and B0 or D and A0 molecules. All of these H-bonds are of a-type (see Fig. 2) with A, B0 , C and D0 as proton acceptor and A0 , B, C0 and D as proton donor. The condition for the occurrence of coupling is the same symmetry of modes. The following modes can couple: W(A)-W(C) with R(B)-R⁄(D), W(A)-W⁄(C) with R(B)-R(D), W(B)-W(D) with R(A)-R⁄(C) and W(B)-W⁄(D) with R(A)-R(C). Table 6 presents the symmetry and behavior of the dipole moment for modes formed by the coupling. In the case of absence of the couplings, the particular peak in the IR spectrum is associated with the power spectrum from only one set (set I and II contain molecules labeled A, C, A0 , C0 and B, D, B0 , D0 respectively). Otherwise, if the coupling is present, the IR peak is associated with peaks from power spectra for two sets. This remark is very useful in following analysis which was performed based on the spectra at 4 K. Comparing power spectra with theoretical IR spectrum (Fig. 5a) we conclude that peaks 1, 2, 4, 5 and 6 represent vibrations of molecules from both set I and set II, whereas peak 3 results from vibrations of molecules from only set I. Additional way to perform the assignments is projecting a set of velocities of particular atoms onto eigenvector corresponding to the particular vibrational mode. For this we have chosen single coordinate for each mode and performed FT of the autocorrelation function calculated for velocity connected with this coordinate. The calculated spectrum can be identified as a contribution of considered mode to the total DOS spectrum. In order to construct normal coordinate we selected first a coordinate describing the motion of each molecule (A, B, C, D) within rocking and wagging librations. In both cases we calculated (for each step of simulation) the sum of unit vectors connected with OAH direction. For rocking
Fig. 5. Comparison of theoretical IR spectrum with power spectra from set I and II (a) and with theoretical polarized IR spectra (b) for ice XI at 4 K. The selected peaks are labeled by numbers for further discussion.
108
M. Gług et al. / Chemical Physics 459 (2015) 102–111
Fig. 6. Comparison of theoretical IR spectrum with power spectra from set I and II (a) and with theoretical polarized IR spectra (b) for ice XI at 40 K.
Fig. 7. (a) Subsystems ABCD and A0 B0 C0 D0 in crystal structure of ice XI, (b) three librational modes of water molecule, (c) positive directions of wagging and rocking modes of water molecules for subsystem ABCD.
Table 5 Symmetry of rocking and wagging librational modes of four molecules subsystems (A, B, C, D or A0 , B0 , C0 , D0 ) with the changes of dipole moment. The symbol ‘‘⁄” indicates that phase is reversed (direction of motion is opposite to that from Fig. 7c). The direction of change of the molecular dipole moment is given in // symbols. The DP means change of the total dipole moment. ‘‘W” and ‘‘R” refer to ‘‘wagging” and ‘‘rocking” respectively. Molecule
A B C D Direction of DP
Wagging
Rocking
A2
B1
B2
A1
A1
B2
B1
A2
W/a/ – W/a/ – 0
W/a/ – W⁄/a/ – a
– W/b + c/ – W/b c/ b
– W/b + c/ – W⁄/b + c/ c
R/b + c/ – R/b + c/ – c
R/b + c/ – R⁄/b c/ – b
– R/a/ – R/a/ a
– R/a/ – R⁄/a/ 0
libration, the resulting vector is projected onto equilibrium plane of molecule. Subsequently, the angle between tangential component and equilibrium sum of unit vectors was calculated. For the wagging case, the procedure was similar, with one
exception – projection was performed onto equilibrium symmetry plane of water molecule. The signs of the angles were consistent with the directions presented in Fig 7. We used OAH unit vectors to ensure that the angles were not varied as an effect of bending
109
M. Gług et al. / Chemical Physics 459 (2015) 102–111
Table 6 Symmetry of coupled rocking/wagging librational modes in ice XI unit cell with the change of dipole moment. The symbol ‘‘⁄” indicates that phase is reversed (direction of motion is opposite to Fig. 7c). The direction of change of the molecular dipole moment is given in // symbols. The DP means change of the total dipole moment. ‘‘W” and ‘‘R” refer to ‘‘wagging” and ‘‘rocking” respectively. Molecule
A B C D Direction of DP
Rocking/wagging A2 in phase
A2 out-of-phase
B1 in phase
B1 out-of-phase
B2 in phase
B2 out-of-phase
A1 in phase
A1 out-of-phase
W/a/ R/a/ W/a/ R⁄/a/ 0
W⁄/a/ R/a/ W⁄/a/ R⁄/a/ 0
W/a/ R/a/ W⁄/a/ R/a/ a
W⁄/a/ R/a/ W/a/ R/a/ a
R/b + c/ W/b + c/ R⁄/bc/ W/bc/ b
R⁄/bc/ W/b + c/ R/b + c/ W/bc/ b
R/b + c/ W/b + c/ R/b + c/ W⁄/b + c/ c
R/b + c/ W⁄/bc/ R/b + c/ W/bc/ c
and stretching modes. The appropriate projections eradicate influence of rocking motion on wagging coordinate and vice versa. The coordinates of simple modes (Table 5) and coupled modes (Table 6) are considered as the sum of angles describing particular librations for given molecules. In case of opposite direction (indicated by ‘‘⁄” symbol) the ‘‘ + ” is changed to ‘‘”. In Figs. 8–11 we
present the power spectra of the modes described above grouped on the base of their symmetry. Peak 1 (at about 620 cm1) overlaps with the peak in the a-polarized spectrum. This indicates that this peak is derived from the mode which changes a-component of the dipole moment. Based on the Tables 5 and 6 we conclude that the origin of this band is in the mode W(A)-W⁄(C) or R(B)-R(D) or coupled mode W(A)-R(B)-W⁄ (C)-R(D) or W⁄(A)-R(B)-W(C)-R(D). As previously stated (using the power spectra), this peak is the result of coupling of vibrations of molecules from both sets. For this reason peak 1 can be assigned to one of the coupled modes. In the case of peak 5 (945 cm1), the conclusion is analogous. Assigning W(A)-R(B)-W⁄(C)-R(D) mode to peak 5 is consistent with experimental Raman spectroscopic study [27] and model calculations [59]. Then the W⁄(A)-R(B)-W(C)-R(D) mode can be assigned to peak 1. Fig. 8 confirms assignments for peaks 1 and 5. Spectra originating from W(A)-W⁄(C) and R(B)-R(D) contain two shared peaks – at 620 and 945 cm1. This fact indicates that these modes can be coupled. In order to determine properly in phase/out-of-phase mode coupling, two top spectra were drawn. Based on positions of main peaks, we conclude that the peak at 620 cm1 originates from the W⁄(A)-R(B)-W(C)-R(D) mode, whereas the peak at 945 cm1 is connected with the W(A)-R(B)-W⁄(C)-R(D) mode. Additionally, there is a peak at 860 cm1 in the W(A)-W⁄(C) spectrum and at 760 cm1 in the R(B)-R(D) spectrum. These peaks have not corresponding equivalent in another spectrum. It is highly probable that these peaks are the remains of simple modes W(A)-W⁄(C) and R(B)-R(D). This conjecture complies with the fact that coupling of two vibrations leads to one mode characterized by lower frequency and the other by higher frequency.
Fig. 9. DOS functions for B2 modes. Black line – R(A)-R*(C), red line – W(B)-W(D), green line – R(A)-W(B)-R*(C)-W(D), blue line – R*(A)-W(B)-R(C)-W(D). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 10. DOS functions for A1 modes. Black – R(A)-R(C), red – W(B)-W*(D), green – R (A)-W(B)-R(C)-W*(D), blue – R(A)-W*(B)-R(C)-W(D). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 8. DOS functions for B1 modes. Black – W(A)-W*(C), red – R(B)-R(D), green – W (A)-R(B)-W*(C)-R(D), blue – W*(A)-R(B)-W(C)-R(D). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
110
M. Gług et al. / Chemical Physics 459 (2015) 102–111
Fig. 11. DOS functions for A2 modes. Black – W(A)-W(C), red – R(B)-R*(D), green – W(A)-R(B)-W(C)-R*(D), blue – W*(A)-R(B)-W*(C)-R*(D). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Peak 3 (at about 710 cm1) is associated with the change of b-component of the dipole moment. It may result from R(A)R⁄(C) or W(B)-W(D) or coupled R(A)-W(B)-R⁄(C)-W(D) or R⁄(A)W(B)-R(C)-W(D) vibrations. Based on power spectra, one can conclude, this peak originates from set I. This facts leads to conclusion that peak 3 can be assigned to the R(A)-R⁄(C) mode. Peak 6 (at 1015 cm1) is related to the change of the dipole moment along b direction. Power spectra of both sets indicate that the origin of this peak can be a coupled mode. Based on the analysis of Fig. 9, similar to analysis of Fig. 8, we can state that peak 6 is associated with the coupled mode R⁄(A)-W(B)-R(C)-W(D) and confirm that peak 3 is related to R(A)-R⁄(C). Furthermore, the R(A)-W(B)-R⁄(C)-W(D) mode is localized at about 620 cm1. This means that peak 1 originates from two modes W⁄(A)-R(B)-W(C)-R(D) and R(A)-W(B)R⁄(C)-W(D). Peak at about 900 cm1 corresponds to the W(B)-W (D) mode (as for peak 3, it is a remain of the coupling). This peak could be present in IR spectra because of change of the dipole moment for the W(B)-W(D) mode. However, the intensity is insufficient to distinguish it. Peaks 2 and 4 (at about 700 and 925 cm1) are related to the change of the dipole moment in c-direction. It may result from W(B)-W⁄(D), R(A)-R(C), R(A)-W(B)-R(C)-W⁄(D) or R(A)-W⁄(B)-R (C)-W(D) vibrations. Based on power spectra, they origin from set I and II. This facts leads to the conclusion that these peaks are associated with coupled mode. Based on Fig. 10, we can conclude that peak 2 corresponds to R(A)-W⁄(B)-R(C)-W(D) and peak 4 to R(A)-W(B)-R(C)-W⁄(D). In addition to previously discussed bands, the peaks at 800–875 cm1 are present in the power spectra. They have no counterparts in the IR spectra. These peaks can arise from the modes which do not change the dipole moment, i.e. W(A)-W(C), R(B)-R⁄(D) and the coupled modes: W(A)-R(B)-W(C)-R⁄(D) or W⁄(A)-R(B)-W⁄(C)-R⁄(D) (symmetry A2). This conclusion is partly confirmed by spectra in Fig. 11. Another possible origin of these peak can be twist modes. For spectra calculated at 40 K the analysis is more difficult because of greater noises. However the results of analysis are very similar. The main difference is slight shift to lower wavenumbers. In addition, the peak at 1100 cm1 is present in the IR spectrum, in polarized spectrum (along b axis) and in power spectra for both sets. This can be the result of combination of a librational and a translational modes [60,27].
Fig. 12. Comparison of the DOS function obtained from CPMD calculations, IINS spectrum (taken from Ref. [24]), DOS function obtained by model calculation (taken from Ref. [59]), and DOS function obtained from classical MD calculations (taken from Ref. [61]) in the librational region of H2O in the ice XI phase.
3.3. Comparison with other theoretical and experimental data In order to verify the accuracy of our investigation, the entire power spectrum of ice XI at 4 K (which is the sum of all atomic power spectra) was compared, in Fig. 12, with density of states (DOS) observed in the inelastic incoherent neutron scattering (IINS) of H2O in ice XI [24], with DOS obtained from model calculations [59] and with DOS from classical molecular dynamics simulation [61]. In order to facilitate the comparison, our calculated CPMD DOS was presented as the smooth bandshape obtained by convoluting the CPMD resulting data with gaussian functions. There is very good agreement in the positions and relative intensities of our CPMD calculated DOS and IINS spectrum, except the peak at about 800 cm1, which it is the most intensive one in the CPMD DOS. The agreement between DOS obtained by Iwano et al. [59] with our CPMD DOS, as well as with IINS DOS, is good in the higher-energy part, whereas the peak at about 600 cm1 seems to be doubly splitted in the DOS obtained from model calculations. The DOS obtained by Itoh et al. [61] using classical molecular dynamics calculations reveal some discrepancies in peak frequencies, compared with three DOS functions discussed above, however qualitative similarities are obviously seen. The relative peak positions are reproduced well and the entire DOS function is shifted towards higher wavenumbers. The difference in the peak frequencies was explained by larger elastic constants in the classical MD calculations. Analyzing Fig. 12, we conclude that our power spectrum reproduces most features of the experimental IINS spectrum, including both peak positions and relative intensities, and is also in qualitative agreement with other calculated DOS functions. Based on measurement of infrared spectra of KOH doped water at 4–285 K, Arakawa et. al. [22] suggested that ordered ice in space maybe detectable using infrared telescopes. They stated, there are two peaks in librational region of infrared spectrum of ice XI – major at 850 cm1 and minor at 590 cm1 (in the IR spectrum of ice Ih only one broad peak is observed in this region). In order to verify the above thesis we performed CPMD calculations of the IR spectra of ice XI at 4 and 40 K. These spectra show two major peaks in the librational region, in agreement with Arakawa predictions, with similar intensities similar to the experimental ones.
M. Gług et al. / Chemical Physics 459 (2015) 102–111
4. Conclusions Theoretical infrared spectra and DOS spectra of ordered phase ferroelectric ice XI have been calculated at 4 and 40 K using Car–Parrinello molecular dynamics and Kramers–Kronig relations. The detailed analysis of vibrations in the librational region (at 550–1150 cm1) was performed. The peaks present in the IR spectra were assigned to particular modes based on power spectra and polarized IR spectra. Our study shows that not all peaks result from the coupling between librational modes. Comparison presented in Fig. 8 shows that our results are in good agreement with available experimental and theoretical spectra of ice XI. Conflict of interest There is no conflict of interest. Acknowledgments This research was supported in part by PL-Grid Infrastructure. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18]
S. Kawada, J. Phys. Soc. Jpn. 32 (1972) 1442. Y. Tajima, T. Matsuo, H. Suga, J. Phys. Chem. Solids 45 (1984) 1135. T. Matsuo, Y. Tajima, H. Suga, J. Phys. Chem. Solids 47 (1986) 165. Y. Tajima, T. Matsuo, H. Suga, J. Chem. Phys. 82 (1985) 424. M. Jackson, W. Whithworth, J. Chem. Phys. 103 (1995) 7647. M.E. Brown, W.M. Calvin, Science 287 (2000) 107. M.J. Iedema, M.J. Dresser, D.L. Doering, J.B. Rowland, W.P. Hess, A.A. Tsekouras, J.P. Cowin, J. Phys. Chem. B 102 (1998) 9203. H. Wang, R.C. Bell, M.J. Iedema, A.A. Tsekouras, J.P. Cowin, Astrophys. J. 620 (2005) 1027. H.X. Zhao, X.J. Kong, H. Li, Y.C. Jin, L.S. Long, X.C. Zeng, R.B. Huang, L.S. Zheng, Proc. Natl. Acad. Sci. 108 (2011) 3481. H. Fukazawa, A. Hoshikawa, Y. Ishii, B.C. Chakoumakos, J.A. Fernandez-Baca, ApJ 652 (2006) L57. H. Fukazawa, A. Hoshikawa, B.C. Chakoumakos, J.A. Fernandez-Baca, Nucl. Instrum. Methods Phys. Res. A 600 (2009) 279. X.C. Su, L. Lianos, Y.R. Shen, G.A. Somorjai, Phys. Rev. Lett. 80 (1998) 1533. H. Wang, R.C. Bell, M.J. Iedema, A.A. Tsekouras, J.P. Cowin, ApJ 620 (2005) 1027. A.H. Castro Neto, P. Pujol, E. Fradkin, Phys. Rev. B 74 (2006) 024302. X. Fan, D. Bing, J. Zhang, Z. Shen, J.-L. Kuo, Comput. Mater. Sci. 49 (2010) S170. H. Fukazawa, S. Mae, S. Ikeda, O. Watanabe, Chem. Phys. Lett. 294 (1998) 554. S. Kawada, J. Phys. Soc. Jpn. 58 (1989) 295. H. Fukazawa, S. Ikeda, S. Mae, Chem. Phys. Lett. 282 (1998) 215.
111
[19] H. Suga, Thermochim. Acta 300 (1997) 117. [20] K. Furic´, V. Volovšek, J. Mol. Struct. 976 (2010) 174. [21] H. Fukazawa, S. Ikeda, M. Oguro, S.M. Bennington, S. Mae, J. Chem. Phys. 118 (2003) 1577. [22] M. Arakawa, H. Kagi, H. Fukazawa, Astrophys. J. Suppl. 184 (2009) 361. [23] J.C. Li, V.M. Nield, S.M. Jackson, Chem. Phys. Lett. 241 (1995) 290. [24] J.C. Li, J. Chem. Phys. 105 (1996) 6733. [25] K. Abe, T. Miasa, Y. Ohtake, K. Nakano, M. Nakajima, H. Yamamoto, T. Shigenari, J. Korean Phys. Soc. 46 (2005) 300. [26] K. Abe, T. Shigenari, J. Chem. Phys. 134 (2011) 104506. [27] T. Shigenari, K. Abe, J. Chem. Phys. 136 (2012) 174504. [28] M.J. Wójcik, M. Gług, M. Boczar, Ł. Boda, Chem. Phys. Lett. 612 (2014) 162. [29] R. Car, M. Parrinello, Phys. Rev. Lett. 55 (1985) 2471. [30] M.S. Bergren, S.A. Rice, J. Chem. Phys. 77 (1982) 583. [31] S.A. Rice, M.S. Bergren, A.C. Belch, G. Nielson, J. Phys. Chem. 87 (1983) 4295. [32] M.J. Wójcik, V. Buch, J.P. Devlin, J. Chem. Phys. 99 (1993) 2332. [33] M.J. Wójcik, K. Szczeponek, S. Ikeda, J. Chem. Phys. 117 (2002) 9850. [34] F. Li, J.L. Skinner, J. Chem. Phys. 132 (2010) 204505. [35] F. Li, J.L. Skinner, J. Chem. Phys. 133 (2010) 244504. [36] M. Benoit, M. Bernasconi, P. Focher, M. Parrinello, Phys. Rev. Lett. 76 (1996) 2934. [37] M. Bernasconi, P.L. Silvestrelli, M. Parrinello, Phys. Rev. Lett. 81 (1998) 1235. [38] C. Knight, S.J. Singer, J.-L. Kuo, T.K. Hirsch, L. Ojamäe, M.L. Klein, Phys. Rev. E 73 (2006) 056113. [39] M. Schönherr, B. Slater, J. Hutter, J.V. Vondele, J. Phys. Chem. B 118 (2014) 590. [40] M. Brela, J. Stare, G. Pirc, M. Sollner-Dolenc, M. Boczar, M.J. Wójcik, J. Mavri, J. Phys. Chem. B 116 (2012) 4510. [41] J. Stare, J. Panek, J. Eckert, J. Grdadolnik, J. Mavri, D. Hadzi, J. Phys. Chem. A 112 (2008) 1576. [42] W. Kohn, J. Sham, Phys. Rev. 140 (1965) A1133. [43] F.A. Bornemann, C. Schütte, Numer. Math. 78 (1998) 359. [44] R. Ramírez, T. López-Ciudad, P.P. Kumar, D. Marx, J. Chem. Phys. 121 (2004) 3973. [45] N. Marzari, D. Vanderbilt, Phys. Rev. B 56 (1997) 12847. [46] R. de L. Kronig, J. Opt. Soc. Am. 12 (1926) 547. [47] H.A. Kramers, Atti Cong. Intern. Fisici, (Transactions of Volta Centenary Congress) Como 2, 1927, 545. [48] F. Wooten, Optical Properties of Solids, Academic Press, New York, 1972. [49] T.K. Hirsch, L. Ojamäe, J. Phys. Chem. B 108 (2004) 15856. [50] CPMD,