Vibrational dynamics and spectroscopy of water at porous g-C₃N₄ and C₂N surfaces

Abstract

Porous graphitic materials containing nitrogen are promising catalysts for photo(electro)chemical reactions, notably water splitting, but can also serve as “molecular sieves”. Nitrogen increases the hydrophilicity of the graphite parent material, among other effects. A deeper understanding of how water interacts with C- and N-containing layered materials, if and which differences exist between materials with different N content and pore size, and what the role of water dynamics is – a prerequsite for catalysis and sieving – is largely absent, however. Vibrational spectroscopy can answer some of these questions. In this work, the vibrational dynamics and spectroscopy of deuterated water molecules (D₂O) mimicking dense water layers at room temperature on the surfaces of two different C/N-based materials with different N content and pore size, namely graphitic C₃N₄ (g-C₃N₄) and C₂N, are studied using ab initio molecular dynamics (AIMD). In particular, time-dependent vibrational sum frequency generation (TD-vSFG) spectra of the OD modes and also time-averaged vSFG spectra and OD frequency distributions are computed. This allows us to distinguish “free” (dangling) OD bonds from OD bonds that are bound in a H-bonded water network or at the surface – with subtle differences between the two surfaces and also to a pure water/air interface. It is found that the temporal decay of OD modes is very similar on both surfaces with a correlation time near 4 ps. In contrast, TD-vSFG spectra reveal that the interconversion time from “bonded” to “free” OD bonds is about 8 ps for water on C₂N and thus twice as long as for g-C₃N₄, demonstrating a propensity of the former material to stabilize bonded OD bonds.

---

1 Introduction

Water at interfaces and on surfaces is of immense significance for it plays a role as a solvent and as a catalyst in many chemical processes.^1–3 In biology, interfacial water at lipid bilayers and proteins contributes to the overall biological activity and selectivity.⁴ Water on semiconductor surfaces offers an alternative to conventional energy sources by means of electro- or photo-catalytic water splitting. Many photo(electro)catalysts contain metals that are expensive and sometimes hazardous.^5–7 In contrast, recently synthesized organic materials like graphitic carbon nitride (g-C₃N₄)^8–10 and C₂N ^11,12 have not only demonstrated their propensity to facilitate water-splitting, but they are also easy to synthesize and non-hazardous to the environment. These quasi-two-dimensional, porous, layered materials offer a large surface area, and their properties can be tuned by the C/N ratio and the size of the pores. Nitrogen is particularly important in these materials, notably for uptake of water, as it considerably increases the hydrophilicity of the material w.r.t. graphite or graphene.^12,13 This makes these materials also interesting in a non-catalytic context, e.g., as sieves to separate liquids when acting as water-selective membranes.

Both for catalytic and non-catalytic purposes, like selective adsorption and separation, a precise understanding of the interaction of water with C/N-containing porous materials is desirable. Several theoretical studies have been performed that explored the geometric and electronic structure of these carbon-based materials and how the latter is modified by water adsorption.^{12,13,15–18} Often (but not always) single- or few-layer models were studied, and low coverages were studied using stationary models. For example, in a first-principles density functional theory (DFT) study including van der Waals corrections and employing a (3 × 3) periodic supercell model, it was found that a single water molecule preferentially adsorbs in a triangular pore of the g-C₃N₄ surface, with an adsorption energy between −0.55 and −0.62 eV depending on the functional.¹⁸ The water molecule lies almost flat in the pore of the buckled g-C₃N₄ surface, forming two hydrogen bonds with acceptor N atoms of the cavity. Also, cluster models using highly correlated wavefunction methods support these results.^13,14 (Note that with smaller (1 × 1) unit cells, the surface is falsely predicted to be flat and other adsorption geometries are more stable.¹⁶) In ref. 18, other (metastable) adsorption sites for water were also found using a potential energy surface (PES) scan, and barriers for lateral diffusion of water to these sites were determined, which were in the order of 0.2 eV. In ref. 17, in an attempt to go to higher-coverage situations, multi-layer systems and dynamical behaviour relevant for water-selective membranes and “molecular sieves”, the diffusion of water (and ethanol) through pristine and modified g-C₃N₄ sheets was studied using classical molecular dynamics simulations based on empirical forcefields. An inverse correlation between the diffusion coefficient of water and the lifetime of intermolecular hydrogen bonds (HBs) was found, i.e., short lifetimes caused larger diffusion coefficients (and vice versa). Diffusion coefficients (for water) on/between g-C₃N₄ were also found to be smaller than for bulk water, and the HB-lifetime was correspondingly larger for water at/in g-C₃N₄ compared to bulk water.¹⁷

In the case of C₂N, which consists of flat layers under ideal conditions and with hexagonal-shaped pores slightly larger than those of g-C₃N₄, both periodic DFT and cluster calculations predict a similar adsorption energy for a single H₂O molecule on a C₂N monolayer as for g-C₃N₄. In this case, however, the water resides above the (center of the) pore, but still forming two H-bonds to two adjacent, basic N atoms.^12,13 Beyond single-layer models for C₂N and after loading with larger amounts of water, according to the combined experimental–theoretical work in ref. 12, water molecules reside between the layers, forming H-bonds with N atoms of C₂N that are stronger than the H-bonds between different water molecules in bulk water. The adsorption energy decreases with increasing water load.

In order to understand the interaction with and adsorption of liquid water on 2D C/N-containing materials better, it is necessary to consider high coverages, finite temperatures, and dynamical behaviour that will naturally arise at a finite temperature. (In the studies mentioned above, this was only done within the classical MD work reported in ref. 17.) Our goal here is to unravel the differences of the C/N material/water systems compared to bulk water or water/air interfaces, and also the differences between various C/N materials. This is, in the end, a prerequisite to fully understanding the performance of N-containing graphitic materials for water splitting or as sieves, and their rational design. Another motivation for our work is the mentioned observation that the H-network, and also the transport and chemical behaviours of water, in confined spaces/at surfaces can be very different from what is observed in the bulk (see also ref. 19 and 20).

Appropriate “fingerprints” are needed by which differences between various systems and general effects of a confining surface can be accessed. In this context, we note that the spectral signature peaks for the concerted intermolecular motion of hydrogen-bonded water molecules, as well as the intramolecular fluctuations of the OH (or OD) modes of water molecules, can be seen in vibrational spectra, such as the infrared (IR) spectrum. In particular, the IR spectrum provides a fingerprint region to decipher the hydrogen-bond dynamics in aqueous solutions and at interfaces.³ Nevertheless, conventional vibrational spectroscopy cannot disseminate the contribution of interfacial molecules from the bulk. Vibrational sum frequency generation spectroscopy (vSFG) enables us to selectively study the structure, hydrogen-bond network and orientational profile of interfaces and surfaces.^21–28 Moreover, the temporal evolution of the interfacial molecules can be explored using time-dependent vSFG (TD-vSFG).^29–31

Existing theoretical studies have not explored the vibrational dynamics of water molecules on two-dimensional, N-containing graphitic materials. In the present work, we have investigated the vibrational dynamics and vSFG spectra of interfacial water on g-C₃N₄ and C₂N surfaces using ab initio molecular dynamics. We concentrate on surfaces in this work, using single-layer models for the 2D material, and leaving investigations of the effects of other layers and the bulk to future work. We will consider thin layers of D₂O molecules (rather than H₂O, for practical reasons) adsorbed on one side of the 2D layer, and analyze OD vibrations in detail. Fluctuations in the vibrational frequency of OD modes were obtained using the wavelet transform of the time-series analysis and the surface-specific velocity–velocity autocorrelation function (ssVVAF) approach³² was used to calculate the vSFG and TD-vSFG spectra of interfacial water molecules.

The paper is organized as follows. In the following Section 2, we briefly outline theoretical methods and computational models. In Section 3, we study the vibrational dynamics, time-averaged vibrational spectroscopy and time-resolved vibrational spectroscopy. Section 4 summarizes and concludes our work.

2 Computational details and models

2.1 Models

We performed AIMD simulations for D₂O molecules adsorbed either on C₂N or g-C₃N₄ monolayers, using periodic supercell models. For g-C₃N₄, a supercell of dimensions (20.4 × 20.4 × 24.0) Å with 126 atoms (C₅₄N₇₂) was employed representing a (3 × 3) 2D trigonal elementary cell to model the wave-like reconstruction found elsewhere^13,18 (see Fig. 1(a) in ref. 13 for a pictorial representation). We then adsorbed 60 D₂O molecules on one side (the surface) of the monolayer. A snapshot resulting from an AIMD simulation (see below), is shown in Fig. 1(a). As a result, about 14 Å of vacuum separate different periodic images along the z-axis (perpendicular to the surface).

Fig. 1 Snapshots of water (D₂O) molecules simulated on (a) a C₃N₄ (3 × 3) monolayer and (b) a C₂N monolayer. The employed unit cells are shown in side-views, with C: green, N: blue, O: red, D: white.

C₂N was modelled using a primitive unit cell of dimensions (8.32 × 8.32 × 20.0) Å with 18 atoms (composition C₁₂N₆), corresponding to a 2D trigonal elementary cell, (see Fig. 1(c) in ref. 13). On one side of this layer, 12 D₂O molecules were added, as shown in the form of an AIMD snapshot in Fig. 1(b). This way, a vacuum gap of around 12 Å remained between the water-covered layers.

2.2 AIMD simulations

First-principles simulations of water molecules on the g-C₃N₄ and C₂N surfaces were performed using the Vienna Ab initio Simulation Package^33,34 (VASP, version 5.4). The inner core electrons were represented by projector augmented-wave (PAW) pseudopotentials^35,36 and valence electrons using the Perdew–Becke–Ernzerhof (PBE) exchange correlation functional³⁷ together with Grimme's D3 dispersion correction.^38,39 (For a general assessment of DFT for water, see ref. 40.) The plane-wave kinetic energy cutoff was set to 400 eV for both systems. For C₂N/water, a 3 × 3 × 1 k-point grid was used for Brillouin zone sampling; for g-C₃N₄/water, only the Γ-point was included.

Simulations were performed in the NVT ensemble using the Nosé–Hoover thermostat⁴¹ at 300 K with deuterium masses for hydrogen and a time step of 1 fs for integrating the equations of motion. The vibrational spectra were sampled by using single, long NVT trajectories of 45 ps length (C₃N₄/D₂O) and 55 ps (C₂N/D₂O). D₂O instead of H₂O was used in order to allow for longer time steps and total propagation times in AIMD, hence better statistics, and also to inspire experiments that are often done with deuterated water.⁴²

2.3 Frequency distributions

A (time-dependent) distribution of vibrational frequencies of OD modes of water molecules on g-C₃N₄ and C₂N surfaces was determined using the wavelet transform of a time-series analysis. The method is well documented⁴³ and is based on the principle that a time-dependent function (f(t)) can be expressed in terms of basis functions obtained by the translations and dilations of a mother waveletwhich is represented in the so-called Morlet–Grossman form in our present study and is mathematically given asThe parameters λ and σ are assigned the values of 1 and 2 s⁻¹ in our present study. The coefficients of the wavelet expansion are given by the wavelet transform of f(t), i.e.,where a and b are both real quantities. Here, a > 0 is a scale parameter that is directly related to the instantaneous frequency content of the system over a small time-window centered around b. Accordingly, based on the fluctuations in the time series f(t), the wavelet transform L_ψf(a, b) provides the vibrational frequency for the given small time-window around t = b. Since the frequency value is proportional to the time-window narrows for high frequency (small a) and widens for low frequency (large a). The value of parameter a that maximizes the modulus of the wavelet transform of time series f at time t = b is used to calculate the most important frequency component for the given interval.

The time series of interest in this work, f(t), is constructed as a complex function with its real and imaginary parts corresponding to fluctuations of the bond length and momentum of an OD mode projected along the OD bond, namely (possible unit factors omitted)

f(t) = δr^OD(t) + iδp^OD(t). (4)

Here, the fluctuations are defined as

δq(t) = q(t) − 〈q〉, (5)

where q(t) denotes the instantaneous property (= r^OD(t) or p^OD(t)), and 〈q〉 the corresponding average over the trajectory. Here, r^OD(t) = |r̲^OD(t)| where r̲^OD(t) = r̲_O − r̲_D is the OD bond vector for a given OD bond, obtained from atomic position vectors. Furthermore, the momentum projected along the OD mode is given aswhere m_D, m_O are atomic masses of deuterium and oxygen, v̲_O, v̲_D are atom velocities, and the unit vector along the OD bond. This method is then applied to all the OD modes present in a given system. For further details and recent applications of the method within AIMD, see ref. 43–46.

3 Results

3.1 Vibrational frequency distribution and dynamics

3.1.1 Time-averaged frequency distribution. We first look at time-averaged vibrational frequency distributions P(ω) of OD modes, obtained from time-averaging the time-dependent vibrational frequencies of OD modes resulting from the wavelet transform of f(t). In practice, an averaged vibrational frequency distribution is obtained by dividing the frequency range of 2000–3000 cm⁻¹ into equally spaced bins. If at a given instant the frequency of an OD mode falls within a given bin, the bin height is incrementally increased by unity. This process is repeated for each OD mode and along the entire trajectory.

The time-averaged vibrational frequency distribution of OD modes of the interfacial water on the g-C₃N₄ surface was found to be a distribution ranging from about 2000 cm⁻¹ to about 2900 cm⁻¹ with a mean at 2496 cm⁻¹, as shown in Fig. 2(a). The distribution is similar to a distribution found for bulk water at 300 K using a similar methodology (see Fig. 2(a) of ref. 46, for example), however, with two differences. First, the distribution in our Fig. 2(a) is slightly broader, and, more importantly, a clear shoulder peak now appears in the high-frequency region around 2750 cm⁻¹, which is absent in bulk water. Closer analysis shows that this shoulder corresponds to the free or dangling OD modes, i.e., OD bonds that are not part of a H-bond network and not bound to the surface and stick out into vacuum instead. Similarly, the average vibrational frequency of OD modes on the C₂N surface is 2487 cm⁻¹ (Fig. 2(b)), and the whole time-averaged distribution is similar to the other surface. Also for C₂N, a shoulder peak in the high-frequency region around 2750 cm⁻¹ is found, again mainly due to the free/dangling OD modes of the interfacial water. Due to the non-centrosymmetric environment at the surface, the tetrahedral hydrogen-bond network is distorted for both surfaces and thus the population of free OD modes shows an overall increment. The high-frequency shoulder is more prominently observable for g-C₃N₄, which implies a comparatively higher population of free/dangling OD modes. Interestingly, high-frequency shoulders are also a characteristic feature of vSFG (vibrational sum frequency generation) spectra of other systems with interfacial water.^{24–26,42,47}

Fig. 2 Time-averaged frequency distribution P(ω) (in arbitrary units) of the OD modes of water molecules on the (a) g-C₃N₄ and (b) C₂N surfaces.

3.1.2 Time-resolved frequency distribution and dynamics. The temporal dynamics of the OD modes of water molecules on the surface of g-C₃N₄ and C₂N were studied using time-dependent joint probability distributions. The time-dependent joint probability distribution used can be expressed mathematically as

P(ω₃,t₂,ω₁) = 〈δ(ω(t₂) − ω₃)·δ(ω(0) − ω₁)〉. (7)

Here, δ(a(t) − b) = a(t) − b in analogy to eqn (5), and the average is taken over all OD bonds. The joint probability distribution gives the probability that a given OD mode that was oscillating with the vibrational frequency ω₁ evolves to a frequency ω₃ within a time interval t₂. (Note that for centrosymmetric systems other than those studied here, P(ω₃,t₂,ω₁) is similar to a two-dimensional IR (2D-IR) spectrum at long waiting times t₂.⁴⁶)

In Fig. 3, we have plotted the distributions for the waiting times t₂ = 10, 100, 750 and 3500 fs for OD modes of water molecules on the g-C₃N₄ surface. For small waiting times, i.e., t₂ = 10 fs, the distribution is predominantly a straight line aligned along the diagonal ω₁ = ω₃, indicating that the OD modes keep close to their initial vibrational frequency. However, for the longer waiting times, the distribution eventually gets elongated and evolves into a spherical distribution, which implies that the OD modes undergo vibrational spectral diffusion due to hydrogen-bond rearrangement. Moreover, within 3.5 ps the OD modes have sampled all accessible vibrational states/frequencies (cf.Fig. 2(a)). A key observation in the time-resolved frequency distributions is that the high-frequency shoulder is evidently observable at low as well as long waiting times. Furthermore, the vibrational dephasing of the modes corresponding to the shoulder peak is seen progressing as a vertical slab along 2750 cm⁻¹.

Fig. 3 Joint probability frequency distribution P(ω₃,t₂,ω₁) of OD modes of water molecules on the g-C₃N₄ surface for waiting times t₂ = (a) 10, (b) 100, (c) 750 and (d) 3500 fs.

A similar analysis for the water molecules on the C₂N surface corresponding to the waiting times of t₂ = 10, 100, 750 and 3500 fs is shown in Fig. 4. The vibrational dephasing follows a similar pattern, and overall differences between both surfaces seem small.

Fig. 4 Joint probability frequency distribution P(ω₃,t₂,ω₁) of OD modes of water molecules on the C₂N surface for waiting times t₂ = 0, 100, 750 and 3500 fs.

To more quantitatively address this point, the temporal dynamics of the OD modes of water molecules on the g-C₃N₄ and C₂N surfaces is also analyzed using a frequency time correlation function (FTCF), which is given as

The time-dependent decay of the correlation function of OD modes of water molecules on the g-C₃N₄ and C₂N surfaces is shown in Fig. 5(a).

Fig. 5 (a) Frequency correlation function C_ωω(t) and (b) hydrogen bond correlation function S_HB(t) for OD modes of water molecules on the g-C₃N₄ and C₂N surfaces.

Evidently, the decay is biphasic with a short time regime that extends up to 100 fs and is governed by the ultrafast librational dynamics or concerted collective motion of water molecules,⁴⁸ whereas the long time decay extends up to a few picoseconds and is mainly attributed to the hydrogen bond rearrangement dynamics. The timescale of vibrational dephasing was determined by using a bi-exponential fitting function,

The short time constant τ₀ was found to be 76 fs and the long time decay constant τ₁ was 5.7 ps for water on g-C₃N₄. For C₂N, we find τ₀ = 104 fs and τ₁ = 4.4 ps. We note that for both systems the long time constant, τ₁, is significantly longer than for bulk liquid water.^46,49 According to ref. 46, a computed τ₁ for bulk D₂O is 1.6 ps. Thus, vibrational dynamics of OD modes are slowed down in the proximity of the g-C₃N₄ surface, and slightly less so near C₂N.

Since frequency fluctuations of the OD modes are associated with the rearrangement of the local hydrogen bond network, we have also calculated a “hydrogen bond lifetime” using a continuous HB correlation function given as

Here, h(t) and H(t) are hydrogen bond population labels, defined as follows. For a given pair of hydrogen-bonded, different water molecules at an instant t, h(t) is assigned a value of 1 and otherwise 0. (Time zero refers to the beginning of the segment of the trajectory into which the trajectory was divided.) Furthermore, if the hydrogen bond between a pair of water molecules remains intact from time t = 0 till time t, H(t) is assigned a value of 1 or else 0. The water molecules are considered to be hydrogen-bonded if the intermolecular O⋯D distance is less than 2.45 Å, corresponding to the first minimum in the O⋯D pair correlation function. The temporal decay of the hydrogen bond correlation function indicates the timescale of rearrangement of the H-bond network and is shown for water (D₂O) on g-C₃N₄ and C₂N surfaces in Fig. 5(b). Evidently, also within this measure, the temporal decay of the S_HB(t) correlation function is somewhat faster for C₂N/water in comparison to g-C₃N₄/water, as can be seen by visual inspection of Fig. 5(b). More quantitatively, the long hydrogen bond lifetime τ₁, obtained from bi-exponential fits of S_HB(t), are found to be τ₁ = 1.8 ps for g-C₃N₄ and 1.3 ps for water on the C₂N surface. Calculations of S_HB(t) have been done also for bulk deuterated water and for the water–vacuum interface in ref. 50. The H-bond lifetimes computed there are of similar length, but hard to compare quantitatively to our numbers since they were obtained with another DFT functional and with Car–Parrinello MD rather than Born–Oppenheimer MD. Still, in ref. 50 it was found that the bulk S_HB function decays faster than the interface one. Here, we have the additional effect of the porous surface and we see differences between the two surfaces, g-C₃N₄ and C₂N.

3.2 Vibrational sum frequency generation (vSFG) spectra

Further insight into the water dynamics on the studied 2D materials, more directly accessible by experiment, comes from vibrational spectroscopy. Here we concentrate on vibrational SFG as a particularly surface-sensitive method.

3.2.1 Time-averaged vSFG spectra. We start with time-averaged vSFG spectra. A time-averaged vSFG spectrum provides relevant frequencies and gives information about the orientation profile of molecular dipoles at the surface/interface. The time-averaged vSFG spectrum is determined from the second-order optical susceptibility tensor of the system, with tensor elements^27,28Here, a, b, and c are Cartesian coordinates, α_ab(t) is the (a,b)-component of the polarizability tensor at time t, μ_c(0) is the c-component of the dipole vector at time t = 0, and ω is the frequency. Squared, summed (and weighted) components of χ_abc² give the vSFG spectrum. Here, we concentrate on so-called ssp-polarization spectra, which we approximate by only considering a single χ_abc² element and setting a = b = x and c = z, where x is a direction parallel to the surface and z perpendicular to it. We show real and imaginary parts of the second-order susceptibility, χ_xxz²(ω), whereby the imaginary part can be attributed to the orientation of the interfacial molecules with respect to the surface.^21–23 In practice, we perform calculations here using the ssVVAF method mentioned earlier.^32,47 Then, eqn (11) can be rewritten, for our D₂O/surface vSFG problem (with only O–D stretch vibrations considered), as^32,47Here, r̲^OD_i is the coordinate vector of an OD bond (i) as defined earlier (there are M = 2·N OD bonds for N water molecules), v̲^OD_i = dr̲^OD_i/dt is the corresponding velocity, v^OD_c,i is its c-component, and and are the derivatives of the dipole moments and polarizabilities, respectively. A few things are worth mentioning: the ssVVAF expression (12) is approximate in the sense that only O–D stretching vibrations are considered. Furthermore, we also use a parameterized (linear) version of the dipole moments and polarizabilities without explicitly calculating them – in fact we also do not include the frequency dependence of and and set them as constants. Finally, a sometimes used “quantum correction factor” Q(ω) ³² is also neglected. Still, the approach followed here will give insight into the O–D dynamics and in particular into the dynamics of the interfacial H-bond network. The method has been shown to give for the O–H stretch region of a water/air interface quite similar vSFG spectra compared to a non-parameterized procedure with calculated dipole moments and polarizabilities.²⁸

In Fig. 6, we show the real and imaginary parts of χ_xxz² for D₂O at g-C₃N₄ (a) and C₂N, respectively, again in the frequency range of the O–D stretching vibrations. To gain information on the orientation of the OD bonds, we focus on the imaginary part of χ_xxz². (The real part of the vSFG spectrum is similar to a conventional infrared spectrum and provides information about the frequency required to excite vibrations of the adsorbed molecules.) For D₂O/g-C₃N₄, the imaginary part of the susceptibility is characterized by a broad, (mostly) negative peak for the frequency domain of 2400–2700 cm⁻¹, which refers to the OD modes that are hydrogen-bonded with the other water molecules or with the g-C₃N₄ surface (Fig. 6(a)). We also note a sharp positive peak for the frequency domain of 2700–2800 cm⁻¹, referring to free, dangling, or weakly bound OD units (with D₂O not having four D-bonded neighbours), which are orientated away from the surface. For example, OD bonds pointing into vacuum are of this type. Note that positive intensities refer to dominantly “upward” and negative intensities dominantly to “downward” orientations of OD bonds. We also note that the intensity of this peak remains unchanged relative to the (unsupported) pure water/air interface computed in ref. 30 (Fig. 2 there), while the intensity of bonded, red-shifted OD modes is relatively reduced compared to this reference. This may imply an overall reduction of contributing bonded OD bonds, relative to the free/non-bonded ones, in the case of water on C/N surfaces compared to the water/air interface. This is in agreement with the frequency distribution of water molecules on g-C₃N₄, as shown in Fig. 1(a), where the population of free/dangling OD modes is clearly visible as a high-frequency shoulder.

Fig. 6 Real and imaginary parts of χ_xxz² for D₂O at (a) g-C₃N₄ and (b) C₂N surfaces. The y-axis is in arbitrary units, with the largest positive intensity of imaginary or real parts normalized to 1.

Furthermore, in Fig. 6(b) the second-order susceptibility of the water molecules on the C₂N surface is shown. The imaginary part of the second-order susceptibility for this system is more similar to the pure water–air interface.³⁰ Namely, similar to the latter, the negative-intensity peak covering the spectral domain of 2200–2700 cm⁻¹ for the water/C₂N system is broader and of larger intensity (comparable to the high-frequency peak) compared to the water/g-C₃N₄ surface. (Note that in ref. 28 and 30, non-deuterated water was considered so absolute oxygen–hydrogen frequencies are different).

To summarize, the vSFG spectra of OD modes on the g-C₃N₄ and C₂N surfaces have similar features for the high-frequency region around 2700–2800 cm⁻¹. But for the frequency region of 2200–2700 cm⁻¹, the spectral intensity is depleted in the case of g-C₃N₄, while it is not for water/C₂N, which behaves similar to the water–air interface.²⁸ We attribute this behaviour to a propensity of the C₂N material to stabilize bonded OD bonds.

3.2.2 Time-resolved vSFG spectra. Time-averaged vSFG spectra provide information about the orientation of molecular dipoles at the interface. Nevertheless, dynamics of these molecular dipoles cannot be studied with the conventional vSFG technique. Here, time-resolved vSFG spectra offer more insight. A time-dependent vibrational SFG spectrum can be obtained by modifying eqn (11) asHere, we compute the second-order susceptibility of a given vibrational oscillator provided that the oscillator was vibrating at a frequency ω′ at a time instant t′. By systematically varying the value of the waiting time, T_w in eqn (13), the delay time between a pump pulse and the overlapping IR/vis pulses that generate an SFG signal, we can obtain frequency-resolved, time-dependent susceptibilities and from there, TD-vSFG spectra.³⁰ It is important to note that in the given computational framework, a priori information of the instantaneous frequency of the oscillators/OD modes is needed to obtain the TD-vSFG spectrum. This information is obtained from the wavelet transform as described above. Furthermore, eqn (13) can easily be rewritten in the ss-VVAF form, eqn (12), which we use in the following to compute Im χ⁽²⁾(ω; T_w) (actually, Im χ⁽²⁾(ω_IR; T_w, ω_pump)).

To calculate the TD-vSFG spectra of the hydrogen-bonded OD modes, we have applied a broadband IR pump pulse corresponding to the frequency domain of 2200–2700 cm⁻¹. Technically, this was done by computing the susceptibilities of the water molecules based on their vibrational excitation frequency within this frequency range. The imaginary parts of the second-order susceptibility of the hydrogen-bonded OD modes of the interfacial water molecules on g-C₃N₄ corresponding to the waiting times T_w = 0, 400, 2000, 3000 and 4000 fs are shown in Fig. 7(a). The spectrum for the waiting time T_w = 0 fs, the “impulsive limit”, is dominated by a peak with negative intensity centered around 2600 cm⁻¹. With increasing waiting times, e.g. for T_w = 400 fs, we see a formerly shallow shoulder in the high-frequency region around 2750 cm⁻¹ growing, corresponding to the free/dangling OD modes. In addition, within 4.0 ps, the peak has attained the maximum intensity and becomes saturated. The figure demonstrates the interconversion from hydrogen/surface-bonded (“bound”) OD modes to free OD modes taking place on this timescale.

Fig. 7 Imaginary part of χ⁽²⁾(ω) at different delay times T_w, computed for water on g-C₃N₄ for (a) the H-bonded OD modes and (b) the free/dangling OD modes. For C₂N, the same information is given in panels (c) for bonded OD modes and (d) for free/dangling OD modes. Now, for the y-axes, arbitrary but non-normalized units are used to make the time-evolution of the signal intensities visible.

Similarly, we have calculated spectra of free/dangling OD modes of g-C₃N₄ for waiting times T_w = 0, 200, 1000 and 2000 fs in Fig. 7(b). Here we applied a narrowband IR pump pulse in the range of 2700–2900 cm⁻¹. The spectrum for T_w = 0 fs is a single, positive peak centered around 2750 cm⁻¹. With increasing waiting times of T_w = 200, 1000, and 2000 fs, the frequency domain between 2200 and 2700 cm⁻¹ corresponding to the H-bonded region shows finite-valued, negative intensity, while at the same time the high-frequency, positive peak loses intensity. This indicates interconversion of free OD modes to (H-)bound OD modes. Closer inspection shows, as can also be seen from the shorter timescale shown in Fig. 7(b), that the free → bound interconversion proceeds faster (within about 2 ps) in (b) than the bound → free interconversion (about 4 ps) in (a). This is due to the free OD modes having a strong propensity to form hydrogen bonds and stabilize. When compared to the pure water/air interface,^30,31 it is found that the hydrogen-bond interconversion dynamics for the bonded OD modes are overall slowed down, whereas for the free OD, the timescale remains relatively unchanged in the presence of the g-C₃N₄ surface.

Along the same lines, we have also obtained the spectra of the bonded and dangling/free OD bonds of water on the C₂N surface, as shown in Fig. 7(c) and (d). The spectra of the bonded OD modes are shown for waiting times of T_w = 0, 400, 2000, 5000 and 8000 fs. While the spectral features as seen in the TD-vSFG spectra are similar to previous calculations for the water/air interface,³⁰ a key difference is a very slow temporal evolution of the peak for the free OD modes, which extends up to nearly 8 ps.

The spectra of dangling/free OD modes are shown for waiting times of T_w = 0, 200, 1000, 2000 and 5000 fs. Again, the free OD modes have a strong propensity to form hydrogen bonds due to which the spectral domain corresponding to the bonded OD modes proceeds mostly within 2 ps. The spectra of free OD modes of interfacial water molecules are interestingly independent of the surface, i.e., g-C₃N₄ or C₂N.

To summarize, we note that although the overall vibrational dynamics of OD modes on g-C₃N₄ and C₂N are similar according to the frequency correlation function, Fig. 5(a) (with a tendency of faster H-bond dynamics for C₂N), the interconversion rates of free and bonded OD modes can still be significantly different (and of other order). This refers in particular to the bound → free interconversion of OD bonds, where D₂O/C₂N exhibits a significantly longer (∼8 ps) timescale than D₂O/g-C₃N₄ (∼4 ps). Here we want to emphasize that the time-resolved vSFG spectra of hydrogen-bonded OD modes on the g-C₃N₄/C₂N gives a comparative perspective of stabilization of the OD modes on two surfaces. The observed differences in the interconversion rates for the two systems imply that g-C₃N₄ could be more suitable for catalyzing interfacial reactions involving water, e.g., water dissociation because the OD modes have a relatively higher propensity to break the hydrogen bonds (∼4 ps). In contrast, C₂N may be more suited for processes like desalination or sieving, which require stabilization of water by the surface with an average interconversion rate of nearly 8 ps.

The key differences between both systems is that the g-C₃N₄ has a buckled surface whereas C₂N is flat. Moreover, from the stoichiometric perspective, g-C₃N₄ has a higher percentage of nitrogen compared to C₂N. These factors obviously affect the HB dynamics of adsorbed water. Finally, we note that the present inferences are drawn based on the TD-vSFG spectra of water molecules on a single monolayer of g-C₃N₄ or C₂N. The role of a multi-layered surface as well as the density of water molecules corresponding to a liquid phase in the simulation, may bring more insight into the catalytic properties of the two materials. Work along these lines is in progress in our laboratory.

4 Summary

To summarize, we have studied the vibrational dynamics and (time-dependent) vSFG spectra of deuterated water molecules at g-C₃N₄ and C₂N surfaces, respectively, using DFT-based AIMD simulations. The most important results are as follows.

• The time-averaged frequency distribution for the OD modes of water molecules on both surfaces shows a high-frequency peak corresponding to free OD modes. The “bound” OD bonds, either engaged in the intermolecular H-bond network or interacting with the C/N surface, are usually weakened and therefore red-shifted w.r.t. the high-frequency peak.

• The time-dependent decay of the OD frequency (measured using a frequency correlation function) is quite similar for both surfaces with a timescale of about 4 ps, which is considerably slower than that of the bulk liquid (deuterated) water. Similar trends for the two studied C/N surfaces arise when analyzing the H-bond network via a H-bond measure, S_HB(t), which implies a slightly longer timescale for g-C₃N₄ than for C₂N.

• Time-independent vSFG spectra (measured using χ_xxz²), indicate that the dangling OD bonds pointing upward and away from the water/vacuum interface are dominant features in the spectra and are characterized by positive imaginary parts of χ_xxz² at high frequencies. “Bound” OD modes are usually of lower intensity, have negative imaginary parts of χ² and are red-shifted. As a difference between the two surfaces, the spectral domain around 2000–2600 cm⁻¹ of the vSFG spectrum for D₂O/C₂N is marked by a stronger negative intensity peak compared to g-C₃N₄.

• Detailed insight into OD bond dynamics is provided by TD-vSFG spectra, with interesting trends for the two surfaces. Although the vibrational frequency correlation functions are similar, the time-resolved vSFG of bonded OD modes saturates within about 4 ps for g-C₃N₄, whereas the same process takes nearly 8 ps for the C₂N surface. The “free” → “bound” interconversion is faster (∼2 ps) and more similar for both surfaces. Thus we can infer on the basis of slower interconversion rates for the bonded to free OD state that the C₂N surface has a propensity to stabilize bonded OD moieties.

All of these observations imply different performance of porous C/N-containing materials in either catalysis or sieving. That key factors like pore size, C/N ratio, or the three-dimensional structure of the materials (e.g., flat or buckled) influence the vibrational dynamics of adsorbed water has been demonstrated; however, how these factors act in detail is a matter of further research.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

P. S. acknowledges support by the Deutsche Forschungsgemeinschaft (DFG), through project Sa 547-18, and Ch. P. by Deutsche Forschungsgemeinschaft within Germany’s Excellence Strategy – EXC 2008/1-390540038, UniSysCat. The authors thank the Paderborn Center for Parallel Computing (PC²) for computation time on supercomputer “Noctua 2”.

References

1. O. Björneholm, M. H. Hansen, A. Hodgson, L.-M. Liu, D. T. Limmer, A. Michaelides, P. Pedevilla, J. Rossmeisl, H. Shen, G. Tocci, E. Tyrode, M.-M. Walz, J. Werner and H. Bluhm, Chem. Rev., 2016, 116, 7698–7726.
2. Y. Jung and R. A. Marcus, J. Am. Chem. Soc., 2007, 129, 5492–5502.
3. J. L. Skinner, P. A. Pieniazek and S. M. Gruenbaum, Acc. Chem. Res., 2012, 45, 93–100.
4. P. Jungwirth, J. Phys. Chem. Lett., 2015, 6, 2449–2451.
5. G. Gonella, et al. , Nat. Rev. Chem., 2021, 5, 466–485.
6. J. Zhu, et al. , Angew. Chem., Int. Ed., 2015, 54, 9111–9114.
7. N. G. Hörrmann, Z. Guo, F. Ambrosio, O. Andreussi, A. Pasquarello and N. Marzari, npj Comput. Mater., 2019, 5, 100–105.
8. X. Wang, K. Maeda, A. Thomas, K. Takanabe, G. Xin, J. M. Carlsson, K. Domen and M. Antonietti, Nat. Mater., 2009, 8, 76–80.
9. Y. Wang, X. Wang and M. Antonietti, Angew. Chem., Int. Ed., 2012, 51, 68–89.
10. X. Wang, S. Blechert and M. Antonietti, ACS Catal., 2012, 2, 1596–1606.
11. R. Walczak, B. Kurpil, A. Savateev, T. Heil, J. Schmidt, Q. Qin, M. Antonietti and M. Oschatz, Angew. Chem., Int. Ed., 2018, 57, 10765–10770.
12. J. Heske, R. Walczak, J. D. Epping, S. Youk, S. K. Sahoo, M. Antonietti, T. D. Kühne and M. Oschatz, J. Mater. Chem. A, 2021, 9, 22563–22572.
13. C. Penschke, A. Zehle, N. Jahn, R. E. von Zander, R. Neumann, A. Beqiraj and P. Saalfrank, Phys. Chem. Chem. Phys., 2022, 24, 14709–14726.
14. T. Schäfer, A. Gallo, A. Irmler, F. Hummel and A. Grüneis, J. Chem. Phys., 2021, 155, 244103.
15. H.-Z. Wu, L.-M. Liu and S.-J. Zhao, Phys. Chem. Chem. Phys., 2014, 16, 3299–3304.
16. S. M. Aspera, M. David and H. Kasai, Jpn. J. Appl. Phys., 2010, 49, 115703.
17. X. Zou, M. Li, S. Zhou, Ch Chen, J. Zhong, A. Xue, Y. Zhang and Y. Zhao, J. Membr. Sci., 2019, 585, 81–89.
18. J. Wirth, R. Neumann, M. Antonietti and P. Saalfrank, Phys. Chem. Chem. Phys., 2014, 16, 15917–15926.
19. I.-C. Yeh and G. Hummer, J. Phys. Chem. B, 2004, 108, 15873–15879.
20. N. Kavokine, R. R. Netz and L. Bocquet, Annu. Rev. Fluid Mech., 2021, 53, 377–410.
21. Q. Du, R. Superfine, E. Freysz and Y. R. Shen, Phys. Rev. Lett., 1993, 70, 2313–2316.
22. Y. R. Shen and V. Ostroverkhov, Chem. Rev., 2006, 106, 1140–1154.
23. V. Ostroverkhov, G. A. Waychunas and Y. R. Shen, Phys. Rev. Lett., 2005, 94, 046102.
24. A. Morita and J. T. Hynes, Chem. Phys., 2000, 258, 371–390.
25. A. Morita and J. T. Hynes, J. Phys. Chem. B, 2002, 106, 673–685.
26. T. Ishiyama, T. Imamura and A. Morita, Chem. Rev., 2014, 114, 8447–8470.
27. A. Morita, Theory of Sum Frequency Generation Spectroscopy, Springer, Singapore, 2018.
28. D. Ojha and T. D. Kühne, Molecules, 2020, 25, 3939.
29. E. H. G. Backus, J. D. Cyran, M. Grechko, Y. Nagata and M. Bonn, J. Phys. Chem. A, 2018, 122, 2401–2410.
30. D. Ojha, N. K. Kaliannan and T. D. Kühne, Commun. Chem., 2019, 2, 116.
31. D. Ojha and T. D. Kühne, Sci. Rep., 2021, 11, 2456.
32. T. Ohto, K. Usui, T. Hasegawa, M. Bonn and Y. Nagata, J. Chem. Phys., 2015, 143, 124702.
33. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
34. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
35. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
36. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
37. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
38. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
39. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
40. M. J. Gillan, D. Alfè and A. Michaelides, J. Chem. Phys., 2016, 144, 130901.
41. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
42. G. Melani, Y. Nagata, R. K. Campen and P. Saalfrank, J. Chem. Phys., 2019, 150, 244701.
43. A. Semparithi and S. Keshavamurthy, Phys. Chem. Chem. Phys., 2003, 5, 5051–5062.
44. D. Ojha, A. Henao and T. D. Kühne, J. Chem. Phys., 2018, 148, 102328.
45. D. Ojha, K. Karhan and T. D. Kühne, Sci. Rep., 2018, 8, 16888.
46. D. Ojha and T. D. Kühne, Phys. Chem. Chem. Phys., 2023, 25, 13442.
47. G. Melani, Y. Nagata, J. Wirth and P. Saalfrank, J. Chem. Phys., 2018, 149, 014707.
48. M. Heyden, J. Sun, S. Funkner, G. Mathias, H. Forbert, M. Havenith and D. Marx, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 12068–12073.
49. D. Ojha and A. Chandra, Chem. Phys. Lett., 2020, 751, 137493.
50. J. R. Choudhuri and A. Chandra, J. Chem. Phys., 2014, 141, 194705.

---