CHEMPHYSCHEM ARTICLES DOI: 10.1002/cphc.201402371

Ab Initio H2O in Realistic Hydrophilic Confinement Christoph Allolio,[a] Felix Klameth,[b] Michael Vogel,[b] and Daniel Sebastiani*[a] A protocol for the ab initio construction of a realistic cylindrical pore in amorphous silica, serving as a geometric nanoscale confinement for liquids and solutions, is presented. Upon filling the pore with liquid water at different densities, the structure and dynamics of the liquid inside the confinement can be characterized. At high density, the pore introduces long-range oscillations into the water density profile, which makes the

water structure unlike that of the bulk across the entire pore. The tetrahedral structure of water is also affected up to the second solvation shell of the pore wall. Furthermore, the effects of the confinement on hydrogen bonding and diffusion, resulting in a weakening and distortion of the water structure at the pore walls and a slowdown in diffusion, are characterized.

1. Introduction The behavior of water under confinement has been the focus of attention at the interface of physics and chemistry for many years.[1–3] Confinements induce drastic changes in the properties of water compared with those of the bulk. This has theoretically been explored, especially for hydrophobic confinements, such as carbon nanotubes.[4–9] The practical realization of experiments in geometrically well-defined hydrophobic confinements is, however, rather more difficult. For this reason, hydrophilic confinements, especially amorphous silicas, such as MCM-41 or SBA-15, have been a focus of experimental attention.[10–12] Confined water touches on surprisingly fundamental questions; hydration water is often subdiffusive,[13–15] undermining the standard assumptions of white noise and Brownian motion, and leads to revised Langevin and Fokker–Planck approaches.[16] The application of phenomenological bulk thermodynamics-based approaches for microscopic systems has also led to counterintuitive results. At the same time, suppressed crystallization that is observed in confinement has been exploited to study metastable states, such as the longsought-after fragile-to-strong transition of supercooled water.[10, 17–19] However, confined water is far from exotic—it is ubiquitous. Water in a biological cell is heavily confined between proteins and other macromolecules, and confined solvent is an important factor in protein folding and biochemical interfaces.[20] A number of computational studies have been dedicated to study water in amorphous silicates.[13, 18, 21] These studies have mainly used conventional force-field (FF)-based

[a] Dr. C. Allolio, Prof. D. Sebastiani Department of Chemistry, Martin-Luther Universitt Halle-Wittenberg von-Danckelmann-Platz 4 06120 Halle/Saale (Germany) Fax: (+ 49) 0345-5527157 E-mail: [email protected] [b] F. Klameth, Prof. M. Vogel Institut fr Festkçrperphysik, Technische Universitt Darmstadt Hochschulstr. 6-8, 64289 Darmstadt (Germany) Supporting Information for this article is available on the WWW under http://dx.doi.org/10.1002/cphc.201402371.

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

methods for water that are well established and well validated, but actually only based on bulk water reference calculations. We present herein a consistent treatment of confined water at the ab initio level of theory, including a quantum chemical design of the actual silica pore itself. In this context, it is especially relevant to gauge the response to changes in temperature and pressure, to test transferability, and to determine to what extent comparisons to the bulk phase are justified.

Computational Details Herein, three different computational approaches were employed: ab initio molecular dynamics (MD), reactive FF MD, and standard FF MD. For ab-initio MD we used Kohn-Sham density functional theory (DFT), with the Becke-Lee-Yang-Parr[22,23] (BLYP) exchangecorrelation functional and dispersion corrections by Grimme[24] . We employed the Gaussian and plane waves (GPW) method[25,26] , with Goedecker-Teter-Hutter pseudopotentials[27,28] a cutoff of 280 Ry and the TZV2P basis sets for all atom species, including the MOLOPT[29] basis set used for silicon. The self-consistent field (SCF) convergence was set to 6  107 a.u. and an orbital transformation with iterative refinement (OT/IR)[30] algorithm was used to converge the Kohn–Sham (KS) orbitals. The system was thermostatically controlled to 350 K by using a CSVR thermostat[31] with a time constant of 100 fs. The increased temperature with respect to standard conditions was meant to compensate for the overstructuring of water.[32–34] We used the deuteron mass for the hydrogen atoms and a time step of 1 fs. All ab initio simulations were performed using the CP2K code[35] .Despite the use of a stochastical thermostat in the ab initio trajectory, diffusion coefficients and power spectra were not expected to be strongly affected by thermostatic control.[31] For the reactive FF simulations, we employed the ReaxFF FF and parametrization by van Duin and others.[36–38] In particular, we used the SERIALREAX implementation and a time step of 0.95 fs for SiO2 and 0.65 fs for the filled pore. The adjustment of the time step was due to fast OH vibrations introduced when inserting the water molecules. The simulation software package GROMACS was utilized to perform standard FF MD simulations.[39] For H2O, we applied the extended simple point charge (SPC/E) model;[40] interactions between the silica pore and H2O were treated within the framework of Brdka and Zerda.[41] A time step of ChemPhysChem 2014, 15, 3955 – 3962

3955

CHEMPHYSCHEM ARTICLES

www.chemphyschem.org

1 fs was employed. The temperature was controlled by means of a Nos–Hoover thermostat.[42] The system was equilibrated for about 1 ns at temperatures higher than 240 K; at lower temperatures we used up to 10 ns. Data production runs extended up to 100 ns, depending on temperature. The classical MD simulations were performed with a fourfold replication of the silica backbone in the z direction.

2. Construction of the Silica Pore Model The construction of the silica pore started with a high-temperature ReaxFF simulation of a silica melt: A slice of trigonal aquartz crystal structure was transformed into an orthorhombic cell of 39.31  39.31  10.27 . The projection reduced the density to 2.41 g cm3 and 0.024 silicon atoms per 3. This density was somewhat higher than the usual values used for fused silica (2.2 g cm3) or cristobalite (2.33 g cm3 as in P. Gallo’s MCM-41 model),[13] but lower than quartz (2.6 g cm3). Molten quartz was tempered down to 600 K from 5655 K within 23 ps and then to 200 K within another 10 ps. To model experimental templated amorphous silica pore geometries such as MCM-41 or SBA-15,[11] we then cut out a hexagonal piece of the pore. To allow a diameter of the final pore of 2.2 nm, we chose a height of 22 . The 195 H2O molecules (corresponding to a density of ca. 1.5 g cm3) were then inserted into the pore. The water molecules were left to expand and react with the pore wall for 70 ps in another ReaxFF simulation. During this time, the system was quenched from 2100 K, through 300 K, to 100 K to create a low-energy structure for the following ab initio run and avoid excess hydrolysis. In the final structure, we found 174 H2O molecules, one H3O + ion, and two HO ions. We protonated the hydroxyl ions because the silica surface was acidic and the OH ions were due to the reaction of naked SiO with water. In addition, one H2 molecule formed and one SiH bond was substituted during pore generation. We then switched to ab initio MD. The simulation was continued in the ab initio molecular dynamics (AIMD) setup for 24 ps. After the initial few picoseconds of relaxation, we found 177 H2O molecules, of which one was protonated and a corresponding wall SiOgroup was present. A snapshot is shown in Figure 1. After ab initio relaxation, the pore was re-equilibrated by using classical MD at HD and also at a density of 1 g cm3. To run the classical FF MD, we further sanitized the pore, reprotonated SiOH, and normalized a hexacoordinated SiOH2 complex and SiH to tetrahedral silanols. Using the structure file obtained from ab initio simulations and employing the potential and charges of Brdka and Zerda results in a positive net charge of the silica matrix. Therefore, the charge of the silicon atoms decreased by about 1 % and the charges of bonded oxygen atoms increased by about 1 % to obtain an uncharged system. Apart from this, the silica matrix geometry remained frozen. The classically equilibrated geometry at 1 g cm3 density was again run for 19.5 ps by using AIMD.  2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

Figure 1. A snapshot of the high-density (HD) silica pore after structural relaxation at T = 350 K.

3. Results and Discussion 3.1. Characterization of the Pore Features The amorphous silica pore obtained by virtual synthesis was compared with the experimental features of real amorphous silicates.[11] We computed the solvent accessible surface using the tool included in the computer program vmd, with the probe radius corresponding to the Brunauer–Emmett–Teller (BET) area of the N2 molecule. This allowed us to determine the surface SiOH density of the pore. It was found to be 3 nm2

. This is in good agreement with experimental data for MCM41.[12] We also examined the connectivity of the silicon atoms in the interface. A SiO2 silicon atom without silanol groups is called a Q4-type atom; for each silanol group attached to silicon, the index is decreased by one. The ratio of Q3/Q2 is therefore a measure of the grade of hydrolysis and roughness of the surface. We found a Q3/Q2 ratio of seven, which corresponded more closely to the more rugged SBA-15 type of amorphous silica (ratio of 11) because MCM-41 contained very few Q2-type silica groups.[11] We counted exclusively silanol groups at the interface. This more rugged surface structure also exhibits hydrogen bonding between silanol groups on different silicon atoms, as found in SBA-15, but not prevalent in MCM-41.[11] In conclusion, pore surface features are largely within the reasonable expectation values for amorphous silica, but closer to SBA-15 than MCM-41. Q3/Q4 ratios are less important, because we are not modeling a pore unit cell, but simply the interaction with the pore surface. In addition to these pore features, we found a number of defects, including under-/overcoordinated silicon atoms buried in the pore walls. For further details, see the Supporting Information.

ChemPhysChem 2014, 15, 3955 – 3962

3956

CHEMPHYSCHEM ARTICLES 3.2. Water Structure We examined the water structure by using the density profile and hydrogen-bonding coordination number, and proceeded to the tetrahedral-order parameter and associated tetrahedral entropy. The focus of interest in the structural characterization is the range of surface effects and their dependence on density, temperature, and simulation approach.

3.2.1. Water Density The computed density profiles from the density functional theory (DFT) and classical trajectories are given in Figure 2. The

Figure 2. Water density as a function of the distance from the pore wall. We computed the density profile at different pressures at 350 K (top), as well as different temperatures at HD (bottom). Density profiles were obtained from ab initio MD (DFT) and SPC/E (FF) trajectories at a mean density at the pore center of 1.3 (HD) and 1.0 g cm3 [low density (LD)]. The top panel also shows the density of the silica pore from the LD-DFT trajectory, which is very similar to the HD-DFT and FF simulations.

HD simulation, is designed to have a density of about 1.3 g cm3. This corresponds to a pressure of about 2500 bar at room temperature, according to experiments in MCM-41-S by Zhang et al.[17] This pressure regime is of special interest because under these conditions an experimental observation of the fragile–strong transition of the system has been reported. Specifically, a density hysteresis in water was very pronounced  2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

www.chemphyschem.org at such a high pressure and visible even at relatively high temperatures.[17] It should be noted, however, that, according to the IAPWS95[44] equations, which extrapolate from water density/pressure data at 350 K from Bridgman,[45] the bulk-water pressure corresponding to this density should be close to 15 000 bar, bringing it near the liquid–ice VI phase transition. Soper reported the possibility of erroneous experimental densities due to the pore diameter being larger than that reported.[46] Another possible source of errors is the presence of voids and microporosity.[46] However, Zheng et al. reported that they believed the results to be unaffected, despite a possible pore diameter of 25  and microporosity.[47] Our parameters are therefore well suited to help clarify the issue. We report the density profile relative to the pore wall, because in the following we concentrate on interface effects. The pore wall is defined to be 13  away from the center of the box, although some water has diffused further into the wall. The reason for this somewhat convoluted definition is the irregular (rounded hexagonal) geometrical shape of the pore, which prevents us from using an interface atom-referenced coordinate system for volumetric data. Near the chosen origin, a density minimum is present in all of the simulations and the silica density reaches the density of fused quartz used to prepare the simulation. The LD simulation roughly corresponds to atmospheric pressure. At higher temperatures, we found good agreement between BLYP/D- and SPC/E-based density profiles, at both densities (Figure 2, top). Particularly noticeable is a strongly layered structure in the density profile at higher densities, which persists over a large temperature range (Figure 2, bottom). The layering is a long-range confinement effect on water, which is not present at lower pressures. Specifically, at normal density, the oscillations reach a plateau at 6  from the wall (see FF-based density profile at Figure 2, top). We observe that the layered density profile is nearly invariant to changes in the temperature. Furthermore, the density at the border of the pore decreases less rapidly than that observed in a comparable system by Gallo and Rovere.[48] We assign this effect to the irregular shape of the pore and to microporosity in our single pore. By this we refer to cavities in the wall, which adsorb small amounts of water molecules. These cavities are a direct consequence of the partially hydrolyzed character of the pore and the high-pressure regime during the virtual synthesis. Accordingly, we find some H2O to be located deeply in the wall, that is, at more than 13  from the center (compare the slower decay of the HD pore in Figure 2, top). In conclusion, the observation of long-range density oscillations makes the assignment of water-phase behavior based on bulk structure factors, such as that done by Zhang et al.,[17] a questionable undertaking at high pressures. The unreliability of the assumption of homogeneous density in the pore center was also pointed out by Soper.[46] By contrast, no layering effects were observed near neutral walls formed by fixed water molecules.[49] 3.2.2. Hydrogen-Bond Coordination To better understand hydrogen bonding across the pore, the average number of hydrogen bonds per H2O molecule was ChemPhysChem 2014, 15, 3955 – 3962

3957

CHEMPHYSCHEM ARTICLES computed as a function of the distance from the interface. There are several geometrical definitions for hydrogen bonds available in the literature.[50] Our choice was rOO < 3.55 , a(O2 ;O1;H1) < p/6, which is a common definition.[50, 51] We also include silica oxygen atoms and silanol hydrogen atoms into the hydrogen-bond-formation possibilities, using the same definition as that for inter-water hydrogen bonding. The distribution in Figure 3 shows that the hydrogen-bond coordination number reaches the value of 3.7 with increasing distance from the pore silica atoms. This bulk value is reached from a distance of about 4 ; we computed the hydrogen-bond coordination number from a previous trajectory of bulk water at the same temperature (see ref. [33] for details) to be 3.67. What is more striking, however, is that the hydrogen-bond coordination number is only very weakly affected by the high-pressure regime and does not reflect the oscillatory behavior of the density profile. This suggests that the hydrogen bonds of water can adapt their structure to large variations in density.

www.chemphyschem.org The return to bulk-like hydrogen bonding after the first solvation shell of the interface is in good agreement with our previous study of a model silica interface, for which we made a similar observation. Specifically, we found that the 1H NMR chemical shifts of H2O molecules approached the bulk value if the distance from the wall was greater than one solvation shell.[52] Chemical shifts are a good proxy for the strength of hydrogen bonding, and entirely non-hydrogen-bonded water exhibits a strong upfield shift.[53–55] Hydrogen bonding at the amorphous silica interface was also studied in detail in a recent publication by Cimas et al.[56] 3.2.3. Tetrahedral Structure The local water structure can be characterized by a tetrahedral order parameter, Q.[57] It is defined by the angles, yikj, between water oxygen k and its four nearest neighbors i and j. Taking into account all those angles, one obtains Equation (1): Q1

Figure 3. Hydrogen-bond coordination number across the pore in the ab ini-

 3 X 4  3X 1 2 cosyikj þ 8 i j¼iþ1 3

ð1Þ

The four nearest neighbors in this case are not only water oxygen atoms, but also contain the oxygen atoms in the pore walls. The reason for this choice is that the silanol groups in the walls may act as hydrogen-bond donors as well as acceptors. This order parameter is defined such that, in a perfectly tetrahedral environment Q = 1, since cos(109.48) = 1/3. The normalization factor of 3/8 is chosen so that random angular distributions will result in a mean value of hQi = 0.[57, 58] In contrast to the original formulation,[57] the definition in Equation (1) allows for negative values of the order parameter, which occur in very strong deviations from tetrahedral geometries. We computed the probability density distribution of the order parameter as a function of the distance from the wall silica atoms (see Figure 4).

Figure 4. Probability density of the tetrahedral order parameter, Q, as a function from the pore silicon atoms. For the HD-DFT high-pressure simulation (left) and normal-pressure LD-DFT (right), the inset on the right magnifies the order parameter distribution near the wall, as found in the LD-DFT simulation.

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim tio simulations as a function of the distance to silicon.

ChemPhysChem 2014, 15, 3955 – 3962

3958

CHEMPHYSCHEM ARTICLES In the bulk, the tetrahedral order parameter, Q, is known to have two characteristic peaks: one centered at around 0.8 and another centered at around 0.5. These peaks characterize the local water structure as ice-like and disordered, respectively, and their relative weight is known to be correlated with transport properties, such as the diffusion coefficient.[58] The relationship of the height of these peaks is temperature dependent; lower temperatures lead to an increased prevalence of the higher value peak.[57] Our calculations show a drastic deviation from the bulk behavior at the interface, leading to a new peak at 0.2. Spatially, this behavior corresponds to the massive decrease in hydrogen bonding at the interface (see Figure 4, cf. also Figure 3). As the coordination is no longer effectively tetrahedral, the order parameter distribution is shifted to lower values by the angles of the adjacent SiO2 oxygen atoms. Notably, this is not due to overly large distances to SiO2 oxygen atoms, the hydrogen-bond length criterion of 3.55  is always fulfilled for the nearest neighbors. These nearest neighbors are however not all hydrogen-bonding partners; for example, oxygen atoms connected to Q4 silica centers cannot act as hydrogen-bond donors. The peak around Q = 0.2 reflects this rupture of the hydrogen-bond network. However, apart from this particular behavior of interface water, there is also an indirect effect at an increased distance away from the interface (cf. Figure 4). Namely, up to about 4 , the disordered peak is increased with respect to the ice-like peak. It appears that the hydrogen-bond network is slightly disordered by the presence of the interface. In the higher density simulation, “cold” highly structured water coordination is suppressed throughout the pore; an effect that is very similar in effect to an increase in temperature. This behavior is well known from bulk water.[58] We evaluated the direct effect of pressure by comparing the distribution of the tetrahedral order parameter P(Q) in the center of the pore across a range of temperatures, within the classical MD simulation (Figure 5). The results are in good

www.chemphyschem.org gous effect on P(Q) as an increase of about 40 K in temperature. From the distribution of Q, a tetrahedral entropy SQ can be obtained.[59] The tetrahedral entropy has been used to characterize the tetrahedral structure of ab initio water; this mirrors the trends of more sophisticated entropy models.[60] At a given temperature, it may be computed in analogy to the Shannon entropy [Eq. (2)]: 3 SQ  S0 þ kB 2

Z

Qmax

logð1  QÞPðQÞdQ

ð2Þ

Qmin

in which S0 is a temperature/pressure-independent constant, which arises in the derivation of SQ under the assumption of independent Qk for each molecule.[59] The results are shown in Figure 6; they reveal temperature, position, and density de-

Figure 6. Tetrahedral entropy of the DFT water as a function of the distance from the pore silica atoms (top). Evolution of the tetrahedral entropy as a function of temperature for bulk and confined systems, computed by using a FF for water more than 5  from the wall atoms (bottom).

Figure 5. Temperature dependence of the probability distribution of water at the pore center in the HD system. Averaging was performed over all water more than 5  from the wall atoms.

agreement with ab initio MD, upon comparing the room-temperature distributions. The bimodal distribution in the LD pore at a temperature of 300 K is very similar to the one found at higher density at 260 K. The increase in pressure has an analo 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

pendence of the tetrahedral entropy. The position dependence reflects the effects already visible from the distribution of Q. Recent experiments indicate an increase of water entropy inside silica pores,[61] which, from our data, can only reflect interface effects because water in the pore middle at normal density has a bulk-like tetrahedral entropy over a large temperature range (see Figure 6, bottom). The ab initio and classical MD results agree well at high pressures (compare Figure 6, top and bottom). At lower pressures, 350 K DFT water has a lower tetrahedral entropy than 300 K SPC/E water, despite the inChemPhysChem 2014, 15, 3955 – 3962

3959

CHEMPHYSCHEM ARTICLES

www.chemphyschem.org

creased temperature, which is in agreement with the slight overstructuring of DFT water.[60] Experimental observations of water in confinement have postulated a phase transition observable at high pressures, but our simulations do not exhibit any jumps in the tetrahedral entropy in the HD pore. However, in the LD pore, we find a mild change of the temperature dependence near 220 K, which may be related to a fragile-tostrong transition previously reported in simulations of confined SPC/E water at a similar temperature.[18] According to Brovchenko et al., the liquid–liquid critical point of bulk SPC/E water is at 1 1.06 g cm3 and slightly above 200 K,[62] which might explain the disappearance of the kink at high densities. The LD-DFT entropy profile hints at an effect of the wall beyond the first solvation shell, as do lower temperature entropy profiles from the FF-based simulation (see the Supporting Information). Although this indicates a finite size explanation of the observed kink, the change in temperature dependence also occurs to water molecules below the 5  threshold. We judge the final resolution of this puzzle to be beyond the scope of this work. In conclusion, DFT and FF-based molecular dynamics methods show a similar evolution of the water structure across the pore. The density profile at HD shows a confinement effect over relatively large distances. The H2O hydrogen-bond coordination number is affected only at the first water layer of the wall. The tetrahedral entropy is influenced until the second layer. At normal density and, hence, normal pressure, we find possible indication of a fragile-to-strong transition in the form of a kink in the temperature-dependent tetrahedral entropy. We find this kink to be absent at high pressure.

We consider only the one-dimensional diffusion parallel to the pore axis (z direction). The position, z(t), of a H2O molecule is taken to be identical to the z position of its oxygen atom. We first define different regions of the pore: the water in the first hydration layer of the pore wall at a distance of < 4  to the next silicon atom (@walls), and the region in the middle of the pore at a distance > 7  (@middle) to any silicon atom. Results are shown in Figure 7. We find a remarkable slowing of

Figure 7. MSD along the z direction for different layers of D2O, smoothed by window averaging and averaged over 1000 starting points after 6 ps of relaxation, for the HD pore.

diffusion at the wall, while the diffusion in the middle of the pore (in the z direction) is in good agreement with literature values for bulk D2O (see Table 1).[63] The absolute value of the

3.3. Water Dynamics The dynamics of confined water are known to show strong deviations from the bulk. Most prominent are diffusion effects. We investigate the influence of the interface on self-diffusion and velocity autocorrelation spectra. To ascertain interface effects, it is necessary to analyze water dynamics in a spatially resolved manner. Our approach is to simply take the initial positions of the atoms at various time origins, t0, as the “position” of the dynamic quantity we wish to compute. This is only valid for short timescales or slow dynamics. Because exchange between regions is less than 30 % over 10 ps, in our case, it is admissible to make qualitative statements about local dynamic quantities. In the following section, we discuss diffusion constants and frequency spectra from DFT simulations. Because the ab initio simulations use heavy hydrogen atoms, we obtain D2O-like dynamics.

3.3.1. Diffusion The diffusion coefficient is computed from the mean square displacement (MSD) by the Einstein relationship [Eq. (3)]. D¼

1d lim h½zðt0 Þ  zðt0 þ tÞ2 it0 2 dt t!1

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

ð3Þ

Table 1. Diffusion coefficients of the regions of the pore at high and low densities. Region

HD [2 ps1]

LD [2 ps1]

full pore @middle @walls exptl 298 K exptl 308 K

0.08 0.16 < 0.03

0.13 0.23 0.07 0.186[63] 0.239[63]

computed diffusion constant is roughly equal to the experimentally observed one at D2O at 308 K. The value calculated is subject to considerable uncertainty, due to the relatively short simulation time and geometrical finite size effects.[64] Diffusion coefficients of water are not very sensitive to pressure: the small differences we find have the expected trend, since the experimentally measured coefficient decreases at very high densities/pressures.[58, 64] The ratio of the diffusion coefficients for bulk water versus all water in the pore is about two (cf. Figure 1), which is in broad agreement with what has been reported for water confined in mesoporous silica in experimental and theoretical studies.[52, 65] The slowdown of water dynamics near the wall is consistent with results from MD simulations for water in silica pores[49] and neutral confinements.[13] ChemPhysChem 2014, 15, 3955 – 3962

3960

CHEMPHYSCHEM ARTICLES

www.chemphyschem.org

The water at the wall shows indications of being subdiffusive. However, to verify this, the dynamics have to be monitored at longer timescales and diffusion out of the layers has to be taken into account. We attribute some of the slowed diffusion to microporosity of the wall. In conclusion, the ab initio molecular dynamics run delivers values in agreement with expectation; this clearly shows the effects of the interface.

3.3.2. Frequency Spectrum In order to compute the vibrational frequency spectrum of the confined water, we first computed the velocity autocorrelation function (VACF), from the velocities v of the water atoms [Eq. 4]. The VACF C(t) is normalized to the time origin t0, and averaged over all water atoms  CðtÞ ¼

vðt0 Þ  vðt0 þ tÞ vðt0 Þ  vðt0 Þ

 ð4Þ H2 O

The VACF may be Fourier transformed to yield the vibrational density of states. We compute the power spectrum by means of Equation (5), where we have averaged over the time origins t0 : CðwÞ ¼
6 ps of relaxation by using 2 ps sampling windows, respectively. Intensities are normalized so as to be comparable; the spectra computed from all H2O are shifted upwards by 0.2.

When computing C(w), we have further enhanced the spectrum by zero-padding and Blackman–Harris window apodization. The resulting spectrum is given in terms of wavenumbers in Figure 8. Because the ab initio trajectory uses heavy hydrogen atoms, we obtain a D2O frequency spectrum. The high-frequency vibrational bands of water in the center of the pore are in good agreement with experimental bulk D2O vibrational frequency bands obtained from IR spectroscopy, centered at n˜  2500 cm1 for the OD stretching mode and at n˜  2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

This work has been supported by the German Research Foundation (DFG) within the framework of FOR 1583 (grants SE1008/8 and VO 905/9-1). Computing infrastructure was provided by the Northern German Supercomputing Alliance (HLRN) under grant HLRN/bec00073. We also acknowledge the ZEDAT HPC team of the FU Berlin for exclusive access to the Soroban cluster for an entire week. In addition C.A. would like to acknowledge Adri van Duin for providing access to his ReaxFF code and parameters and Tobias Watermann for the BET surface simulation. Keywords: density functional calculations · mesoporous materials · molecular dynamics · phase transitions · water chemistry

ChemPhysChem 2014, 15, 3955 – 3962

3961

CHEMPHYSCHEM ARTICLES [1] M. Vogel, Eur. Phys. J. Spec. Top. 2010, 189, 47 – 64. [2] M. Fayer, N. E. Levinger, Annu. Rev. Anal. Chem. 2010, 3, 89 – 107. [3] I. Brovchenko, A. Oleinikova, Interfacial and Confined Water, Elsevier, Amsterdam, 2008. [4] J. C. Rasaiah, S. Garde, G. Hummer, Annu. Rev. Phys. Chem. 2008, 59, 713 – 740. [5] J. R. Bordin, A. B. de Oliveira, A. Diehl, M. C. Barbosa, J. Chem. Phys. 2012, 137, 084504. [6] R. Zangi, J. Phys. Condens. Matter 2004, 16, S5371. [7] D. J. Bonthuis, K. F. Rinne, K. Falk, C. N. Kaplan, D. Horinek, A. N. Berker, L. Bocquet, R. R. Netz, J. Phys. Condens. Matter 2011, 23, 184110. [8] G. Hummer, J. C. Rasaiah, J. P. Noworyta, Nature 2001, 414, 188 – 190. [9] T. B. Sisan, S. Lichter, Phys. Rev. Lett. 2014, 112, 044501. [10] C. E. Bertrand, Y. Zhang, S.-H. Chen, Phys. Chem. Chem. Phys. 2013, 15, 721 – 745. [11] I. G. Shenderovich, D. Mauder, D. Akcakayiran, G. Buntkowsky, H.-H. Limbach, G. H. Findenegg, J. Phys. Chem. B 2007, 111, 12088 – 12096. [12] B. Grnberg, T. Emmler, E. Gedat, J. Shenderovich, G. H. Findenegg, H. H. Limbach, G. Buntkowsky, Chem. Eur. J. 2004, 10, 5689 – 5696. [13] P. Gallo, M. Rovere, S.-H. Chen, J. Phys. Condens. Matter 2010, 22, 284102. [14] M. Settles, W. Doster, Faraday Discuss. 1996, 103, 269 – 279. [15] A. Bizzarri, S. Cannistraro, J. Phys. Chem. B 2002, 106, 6617 – 6633. [16] R. Metzler, J. Klafter, Phys. Rep. 2000, 339, 1 – 77. [17] Y. Zhang, A. Faraone, W. A. Kamitakahara, K.-H. Liu, C.-Y. Mou, J. B. Leao, S. Chang, S.-H. Chen, Proc. Natl. Acad. Sci. USA 2011, 108, 12206 – 12211. [18] P. Gallo, M. Rovere, S.-H. Chen, J. Phys. Condens. Matter 2012, 24, 064109. [19] F. Mallamace, C. Corsaro, S.-H. Chen, H. E. Stanley in Liquid Polymorphism, Vol. 152 (Ed.: H. E. Stanley), Wiley, Hoboken, 2013, pp. 203 – 262. [20] D. Lucent, V. Vishal, V. S. Pande, Proc. Natl. Acad. Sci. USA 2007, 104, 10430 – 10434. [21] K. Malek, J. Sun, Mendeleev Commun. 2006, 16, 11 – 13. [22] A. D. Becke, Phys. Rev. A 1988, 38, 3098. [23] C. Lee, W. Yang, R. G. Parr, Phys. Rev. B 1988, 37, 785 – 789. [24] S. Grimme, J. Comput. Chem. 2006, 27, 1787 – 1799. [25] G. Lippert, J. Hutter, M. Parrinello, Mol. Phys. 1997, 92, 477 – 487. [26] J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing, J. Hutter, Comput. Phys. Commun. 2005, 167, 103 – 128. [27] C. Hartwigsen, S. Goedecker, J. Hutter, Phys. Rev. B 1998, 58, 3641. [28] S. Goedecker, M. Teter, J. Hutter, Phys. Rev. B 1996, 54, 1703 – 1710. [29] J. VandeVondele, J. Hutter, J. Chem. Phys. 2007, 127, 114105. [30] V. Weber, J. VandeVondele, J. Hutter, A. M. N. Niklasson, J. Chem. Phys. 2008, 128, 084113. [31] G. Bussi, D. Donadio, M. Parrinello, J. Chem. Phys. 2007, 126, 014101. [32] D. Banyai, T. Murakhtina, D. Sebastiani, Magn. Reson. Chem. 2010, 48, S56 – S60. [33] C. Allolio, N. Salas-Illanes, Y. S. Desmukh, M. R. Hansen, D. Sebastiani, J. Phys. Chem. B 2013, 117, 9939 – 9946. [34] C. Allolio, M. Sajadi, N. P. Ernsting, D. Sebastiani, Angew. Chem. 2013, 125, 1858 – 1861; Angew. Chem. Int. Ed. 2013, 52, 1813 – 1816. [35] The CP2K Developers Group, computer code CP2K, available at http:// cp2k.org.

 2014 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim

www.chemphyschem.org [36] A. C. T. van Duin, A. Strachan, S. Stewman, Q. Zhang, X. Xu, W. A. Goddard, J. Phys. Chem. A 2003, 107, 3803 – 3811. [37] J. C. Fogarty, H. M. Aktulga, A. Y. Grama, A. C. T. van Duin, S. A. Pandit, J. Chem. Phys. 2010, 132, 174704. [38] A. C. T. van Duin, S. Dasgupta, F. Lorant, W. A. Goddard, J. Phys. Chem. A 2001, 105, 9396 – 9409. [39] B. Hess, C. Kutzner, D. Van Der Spoel, E. Lindahl, J. Chem. Theory Comput. 2008, 4, 435 – 447. [40] H. Berendsen, J. Grigera, T. Straatsma, J. Phys. Chem. 1987, 91, 6269 – 6271. [41] A. Brdka, T. Zerda, J. Chem. Phys. 1996, 104, 6319. [42] S. Nos, Mol. Phys. 1984, 52, 255. [43] W. Humphrey, A. Dalke, K. Schulten, J. Mol. Graphics 1996, 14, 33 – 38. [44] W. Wagner, A. Pruß, J. Phys. Chem. Ref. Data 2002, 31, 387 – 535. [45] P. W. Bridgman, Proc. Am. Acad. Arts Sci. 1942, 74, 399 – 424. [46] A. K. Soper, Proc. Natl. Acad. Sci. USA 2011, 108, E1192. [47] Y. Zhang, A. Faraone, W. A. Kamitakahara, K.-H. Liu, C.-Y. Mou, J. B. Leao, S. Chang, S.-H. Chen, Proc. Natl. Acad. Sci. USA 2011, 108, E1193 – E1194. [48] P. Gallo, M. Rovere, J. Chem. Phys. 2012, 137, 164503. [49] F. Klameth, M. Vogel, J. Chem. Phys. 2013, 138, 134503. [50] R. Kumar, J. R. Schmidt, J. L. Skinner, J. Chem. Phys. 2007, 126, 204107. [51] A. Luzar, D. Chandler, Phys. Rev. Lett. 1996, 76, 928 – 931. [52] X. Y. Guo, T. Watermann, S. Keane, C. Allolio, D. Sebastiani, Z. Phys. Chem. 2012, 226, 1415 – 1424. [53] H. H. Limbach, Magn. Reson. Chem. 2001, 39, S1 – S2. [54] H. E. Gottlieb, V. Kotlyar, A. Nudelman, J. Org. Chem. 1997, 62, 7512 – 7515. [55] T. Watermann, H. Elgabarty, D. Sebastiani, Phys. Chem. Chem. Phys. 2014, 16, 6146 – 6152. [56] A. Cimas, F. Tielens, M. Sulpizi, M.-P. Gaigeot, D. Costa, J. Phys. Condens. Matter 2014, 26, 244106. [57] P.-L. Chau, A. J. Hardwick, Mol. Phys. 1998, 93, 511 – 518. [58] J. R. Errington, P. G. Debenedetti, Nature 2001, 409, 318 – 321. [59] P. Kumar, S. V. Buldyrev, H. E. Stanley, Proc. Natl. Acad. Sci. USA 2009, 106, 22130 – 22134. [60] C. Zhang, L. Spanu, G. Galli, J. Phys. Chem. B 2011, 115, 14190 – 14195. [61] S. Kittaka, S. Takahara, H. Matsumoto, Y. Wada, T. J. Satoh, T. Yamaguchi, J. Chem. Phys. 2013, 138, 204714. [62] I. Brovchenko, A. Geiger, A. Oleinikova, J. Chem. Phys. 2005, 123, 044515. [63] E. H. Hardy, A. Zygar, M. D. Zeidler, M. Holz, F. D. Sacher, J. Chem. Phys. 2001, 114, 3174. [64] H.-S. Lee, M. E. Tuckerman, J. Chem. Phys. 2007, 126, 164501. [65] K. R. Harris, L. A. Woolf, J. Chem. Soc. Faraday Trans. 1 1980, 76, 377 – 385. [66] K. Shirono, H. Daiguji, J. Phys. Chem. C 2007, 111, 7938 – 7946. [67] O. Bonner, J. D. Curry, Infrared Phys. 1970, 10, 91 – 94. [68] M. Sulpizi, M.-P. Gaigeot, M. Sprik, J. Chem. Theory Comput. 2012, 8, 1037 – 1047.

Received: May 26, 2014 Published online on September 11, 2014

ChemPhysChem 2014, 15, 3955 – 3962

3962

Ab initio H2O in realistic hydrophilic confinement.

A protocol for the ab initio construction of a realistic cylindrical pore in amorphous silica, serving as a geometric nanoscale confinement for liquid...
2MB Sizes 0 Downloads 10 Views