Interface Modeling for the Ab Initio Evaluation of the Water Contact Angle on a Metallic Cu(111) Surface

Controlling the water contact angle on a surface is important for regulating its wettability in industrial applications. Therefore, it is crucial to develop ab initio evaluation methods that can accurately predict this angle. The ab initio predictions require an adsorption structure model for the adsorption of liquid molecules on a surface, but the construction of this model depends on whether the test surface comprises an insulating or metallic material because the surface reconstruction takes quite a di ﬀ erent form in each case. Previous studies have focused on the estimation of the water contact angle on insulators; however, this study eluci-dates the water contact angle on a metallic surface, Cu(111). Because the feasibility of ab initio evaluations depends on the approximation of liquid–gas interface energy, which can be roughly estimated through the interface energy of crystal ice, it is natural to use the periodic-honeycomb array of water molecules as the adsorption model for the water on the surface. However, despite the successful application of the periodic model for ab initio prediction of the water contact angle on insulating surfaces, applying this model to metallic surfaces has not provided satisfactory predictions that reproduce experimental values. Therefore, in this study, we propose the use of models with isolated water oligomers for the ab initio prediction of the water contact angle on a metallic surface, which achieved an accurate prediction. The ambiguity of the models based on the size and coverage of the oligomers was small ( ∼ ± 10 ◦ ), which was averaged out to give a plausible value based on the Boltzmann weight with the adsorbing energies. The proposed procedure can be used to estimate the wettability of the surfaces of other


I. INTRODUCTION
Regulating the wettability of surfaces [1][2][3] is an important issue that must be addressed to broaden their industrial applications, such as heterogeneous catalysis, corrosion, and electrochemistry. [1,4,5] The wettability of the surfaces of metal catalysts such as Cu has a significant effect on their catalytic activity for electrochemical reactions, [6] wherein the receding contact angle of water on the metal surface corresponds to the oxygen reduction reactions of the catalyst. This indicates that the efficiencies of the electrochemical reactions of metal catalysts can be enhanced by controlling the water contact angle with the metal surface.
The contact angle (θ) is the primary measure of wettability; a larger contact angle (θ > 90 • ) indicates hydrophobicity, i.e., low wettability, whereas a smaller contact angle (0 • < θ < 90 • ) indicates hydrophilicity, i.e., high wettability. Experimentally, the water contact angle on metal surfaces can be measured directly by observing the angle captured by the cameras. [7,8] The angle is determined using the Young's relation [9,10]: where γ SG , γ LG , and γ SL are the surface energies at solid-gas, liquid-gas, and solid-liquid interfaces, respectively. Although these energies can be computed by various approaches, the proposed study focuses on ab initio evaluations (the comparison with other approaches is provided later).
The ab initio prediction of water contact angles on surfaces was recently introduced, and it has been successfully applied to estimate this angle on the surfaces of insulators. [1,11] For example, the contact angle of water on Si [11] was predicted via ab initio evaluation (88 • ), which was close to the experimental value (91 • ). In addition, ab initio evaluations have also been applied to estimate the water contact angle on transition metal oxides [1]. For example, the water contact ously reported studies, the γ LG of water was substituted with that of crystal ice because water and crystal ice have similar surface energies [12].
Previously reported studies based on the ab initio evaluations of contact angles of water at different surfaces have considered insulators as the surfaces. In previous studies, [1,11] an "ice-like bilayer model" [13,14] (periodic-honeycomb structure of ice) has been adopted as the adsorption model for conducting the ab initio evaluations of the contact angles of water on various insulator surfaces, which is a natural choice for the γ LG = γ ice model (see Fig.1 (a)). However, when estimating the water contact angle on the surface of metals, the structural modeling to describe the adsorption of water molecules on the metallic surface should be different from that used to describe the adsorption of water molecules on the insulator surface. Previous studies have reported the presence of isolated water hexamers on the surface of metals such as Cu [15,16], indicating that the "bilayer model" is not appropriate for the ab initio evaluation of the water contact angle on the metal surface. Therefore, "isolated models" that do not exhibit the bridging of hydrogen bonds over neighboring unit cells would be more appropriate.
In this study, we construct an isolated model for conducting the ab initio evaluations of the water contact angle on a metallic surface. First, we investigated the advantages offered by the modeling method. We observed that the type of N-mer/oligomer (ranging from monomer to hexamer) and number of N-mers within a unit cell should be specified. Particularly, we examined the dependence of the ab initio prediction on the specific modeling choice. We confirmed that (i) the contact angle of water molecules on the metallic Cu(111) surface, predicted using the ab initio evaluations with isolated modeling, was closer to the experimental value than that obtained with the periodic honeycomb modeling. Further, (ii) the freedom of selection [(N-mer)⊗(coverage)] and the averaging with the Boltzmann weight based on the adsorption energy were useful.
The paper is organized as follows. The "Computational Methods" is divided into two subsections. The "Interface modeling" subsection provides a detailed description of our "isolated cluster" models. The "Computational details" subsection specifies our ab initio methods of evaluating the energies based on the models. The "Results and Discussion" section shows that the numerical results are comparable with available experiment results, and deals with pros and cons of our approach by comparing the proposed approach with other approaches. Finally, the findings are summarized in the "Conclusion" section.

A. Interface modeling
Based on the previous studies, [1,11] we estimated Top view Side view FIG. 1. Structural models describing the adsorption of water on metal surface for the ab initio evaluation of the contact angles. Periodic honeycomb structure known as the "bilayer model." [13,14] Panel (a) shows the bilayer model, which is used for the estimation of the contact angles of water on insulator surfaces. [1,11] Panel (b) and (b') show the "isolated molecular (cluster) models," which were found to be more appropriate for the estimation of the water contact angle on metal surfaces.
where E water ads is the adsorption energy of water molecules on the surface (A is the unit area within which E water ads is defined). For the liquid-gas interfaces, which was evaluated using the surface energy of ice, γ ice . [1,11] Then, the formula for the contact angle can be provided as follows: which can be reduced to the ab initio energy evaluations. E water ads is given as [1]: where E CuSlab , E water (onSurf), and E tot , are the energy of the Cu slab, energy of water molecules on the surface of the adsorbed structure, and total energy of the system with the molecules on the surface, respectively. The surface energy γ ice can be evaluated as: which is a measure of the stabilization of the ice surface, where E ice (n; Slab) and E ice (n; Bulk) are the energies of the slab and the bulk of the ice composed of n water molecules, respectively.
The water contact angle on the surface is estimated using five quantities: E ice (n; Bulk), E ice (n; Slab), E CuSlab , E water (onSurf), and E tot , which can be evaluated using ab initio density functional theory (DFT) calculations. Considering Cu(111) as the surface, we calculated the energies using a DFT package, CASTEP. [17] For the ease of comparison with previous works, [1,11] we used the same exchange-correlation function, GGA-PBE, [18] used in previous studies. In addition, norm-conserving pseudopotentials [19] were used to describe the ionic cores. The detailed computational conditions for each calculation is summarized in Table I.
The detailed information used for modeling the geometries of the adsorbed molecules is provided in the subsequent subsection: Computational Details. The main comparison was conducted between the predictions achieved by the periodic-honeycomb model (bilayer) [13,14] [panel (a) in Fig. 1] and those by the isolated molecular models (buckled) [16,20] [panel (b)], which were evaluated using a Cu slab comprising nine atomic layers. Furthermore, to investigate the existence of the considerable bias on the choice of N-mer and coverage, we compared the predictions with N = 1, 2, 3, 4, and 6; however, this comparison was conducted with reduced cost and complexity of the computation, namely with an H-parallel model [21] (planar model) [panel (b') in Fig. 1] and with reduced number of layers (four). N = 5 was excluded because we limited the possible geometry to satisfy the "on-top alignment," which is supported by the findings of previous studies [20,21] [i.e., N = 5 cannot accommodate the molecules within a unit cell, such that they have an "on-top alignment" geometry]. For the oligomers, N > 2, they can either assume the chain or circular form. Therefore, we considered both possibilities for N = 3 and 4 and considered only the circular form for N = 6, as supported by previous studies. [16,20] The dependence on the coverage, i.e., on the number of N-mers placed within a unit cell, was examined up to N = 3.  [22] and the plane wave energy cutoff, respectively. These values were determined by the convergence of the total energies. 'Cu-struct-opt' and 'Cu-surf-relax' indicate the optimized Cu bulk and Cu-slab structures, respectively.

Computational details
For the energies required to evaluate the contact angles, we prepared the geometries of the Cu slab (for E CuSlab ), ice bulk slab (for E ice (n; Bulk) and E ice (n; Slab)), and the water molecules adsorbed on the Cu slab (for E water (onSurf) and E tot ).
A Cu bulk with an initial lattice constant of 3.6147 Å was prepared to construct the Cu-slab structure. The constant was optimized under the bulk structure as 3.728 Å; subsequently, the bulk was cleaved along the (111) plane and used as the initial structure of the slab. Then, a nine(four)-  Fig. 1], we used the "bilayer model" employed in preceding studies. [13,14] The structure was generated from the "H-parallel" struc- ture, [21] where all the H 2 O molecular planes were parallel to the slab surface. Then, the positions of the oxygen atoms were adjusted such that their heights from the slab surface were alternated in every molecule to form the buckled structure. [16,20] To evaluate the contact angles in the fourlayer slab models (as shown later in Fig. 2), water molecules were placed on the models using the on-top alignment geometry [20,21] with the H-parallel orientation [21]; herein, the vertical height was optimized by DFT.

III. RESULTS AND DISCUSSION
As shown in Table II, the periodic-honeycomb model (bilayer) [1,11] is the most inferior model from the viewpoint of the similarity of results to the experimental values obtained for the water contact angle on the metallic Cu(111) surface. [25,26] In contrast, the isolated molecular models (buckled and planar) [14,16] were more appropriate for predicting the water contact angle on metal surfaces. The difference in their prediction accuracy can be attributed to the bridging of the hydrogen bonds over neighboring unit cells. In addition, because the lattice mismatch between the Cu surface and honeycomb network was −2.1% [13], the difference in the accuracy of the predictions cannot be attributed to the lattice mismatch between the Cu surface and honeycomb network.
As mentioned in the Introduction section, the freedom of choice, with respect to [(N-mer)⊗(coverage)] was achieved for the isolated molecular model. We examined the dependence of the predictions on various models, but with low-computational-cost (four-layer) and complex (planar) structures; the results are shown in Fig. 2. Consequently, the contact angle changes from 79.59 • (buckled/nine layers) → 77.13 • (buckled/four layers) and → 77.29 • (planar/ four -layers), which is considerably small to justify the simplification. Fig. 2 shows the dependence of the predictions on various four-layer models. The dependence range within ±10 • was observed to be around their average value. Therefore, the importance of average value is signified, which considers all the possible choices using the Boltzmann weight based on the adsorption energy, E water ads , as a measure of the stability. The average value (78.47 • ), is similar to the experimental values [25,26], owing to the higher contact angles predicted for smaller N-mers (monomer and dimer).   [29] and Grimme [30] schemes, whereas for the dipole corrections, we applied NSC (nonself-consistent) [31] and SC (self-consistent) [32], as they are available in CASTEP. As shown in Table III, the dispersion corrections lead to a larger magnitude of E water ads as a consequence of the shorter bonding length to the surface. The dispersion correction is then found to reduce the contact angle by approximately 14 • , whereas dipole corrections cause a smaller reduction of approximately 2 • . When we applied both corrections, the contact angle is corrected by approximately δθ ∼ −15 • , thus deviating even further from the experimental value. However, the final conclusion of the present study, namely, that the isolated model is appropriate, would not be affected by these corrections. The reduction δθ is mainly due to the contraction of the bonding length, i.e., the interaction vertical to the surface, whereas the difference between the isolated and bilayer models lies in the intermolecular interactions in the horizontal directions on the surface. The vertical interactions are common to both isolated and bilayer models and hence do not change the conclusion that the isolated model provides values closer to the experimental ones. Because the valence shell of Cu consists of d-electrons, the Hubbard U correction [33][34][35] may also be important. However, this correction matters when the orbital becomes more localized where the self-interaction cancellation is damaged. [36] For our metallic Cu surface, the orbital is delocalized, and hence, we did not consider such corrections, as in most preceding studies. [37] B. Structural model and the evaluation method Our conclusion that the isolated model is appropriate for the metallic Cu(111) surface contradicts the cases of Pd(111), Pt(111), and Ru(0001) surfaces. [38,39] For these surfaces, molecular coverage with horizontal intermolecular interactions is predicted to be the stable surface structure.
This contrast can be attributed to the difference in coverage depending on the adsorption energy.
As explained below, the coverage for Cu(111) becomes lower than those for Pd, Pt, and Ru. Although the lower coverage for Cu is well described by the isolated model, higher coverage is better described by a network of molecules with interactions. The lower coverage is due to the weaker adsorption of molecules on Cu(111), namely, ∼ 352 [meV/H 2 O], as estimated by TPD experiments, [21] which is approximately half of the theoretical estimations for Pd, Pt, and Ru. [38] The molecules adsorbed on Cu then more easily desorb again from the surface, leading to the lower coverage than those for Pd, Pt, and Ru. This lower coverage was actually derived in a preceding study [40] by using a Langmuir adsorption isotherm. [41] The dependence on the adsorption energy appears through the adsorption equilibrium constant as a function of the Boltzmann factor.
Corresponding to the lower coverage, isolated monomers, dimers, trimers, and tetramers have actually been observed in STM (scanning tunneling microscope) experiments on several metallic surfaces. [14] Because the Langmuir adsorption isotherm [41] is basically applicable to solid-gas (steam) . [21] A previous DFT study [40] also concluded that the adsorption of water molecules on Cu(111) without dissociation is much more stable than that with dissociation.
In Eq. (2) of the present study, the left-hand side denotes the balance of forces horizontal to the surface acting as the boundary of the liquid-covered region. By achieving coverage, the system becomes stabilized by E ads , as indicated by the right-hand side of the equation. However, not only this stabilization term but also another term, γ intMol , contributes to pulling the boundary inside toward liquid-covered region owing to intermolecular interactions. This tension term can be omitted when the intermolecular distance increases under the lower coverage, but when the coverage increases, the approximation in Eq. (2) makes the the contact angle less accurate. To justify using this equation in prior studies, [1,11,43] the structural models had to be "sufficiently sparse" to omit γ intMol . Because our isolated model is more sparse than the honeycomb models in those studies, the approximation in Eq. (2) is even more reliable in the present work. Furthermore, using GGA-PBE without considering intermolecular interactions is more consistent with the omission of γ intMol . If we include the dispersion correction, γ intMol would become more pronounced, which would increase the contact angle solely by pulling the boundary inside the wet region. The results shown in Table III, however, are the opposite, i.e., considering the dispersion correction actually reduces the angle. This reduction is caused by another mechanism, as explained above, i.e., the shortening of the vertical bond to the surface. At least within the exchange-correlation potential we used, the horizontal dispersion interactions between molecules thus seem smaller than those for the vertical molecular-surface binding. This is also consistent with omitting γ intMol from Eq. (2). The ab initio simulations with interface modeling appropriate for insulating surfaces ("ice-like bilayer" model) have proven to accurately reproduce the contact angles. [1,11] Furthermore, the present study demonstrates the applicability of the ab initio approach to metallic surfaces. As explained in the present study, our interface modeling ("isolated cluster" model) relies on the experimental facts about how water clusters are absorbed into surfaces. Because water clusters on various metallic surfaces have been investigated comprehensively [47], the interface modeling appropriate for individual metallic surfaces can be achieved without any difficulties. Although the accuracy of ab initio DFT energies generally depend on the exchange-correlation function, but as far as the water absorption on metallic surfaces is concerned, DFT is expected to work efficiently.
As for the molecular levels, several combinations of interface modeling and the molecular simulations have been developed. Within the framework of the molecular simulation, the interface energy is thermodynamically described via fundamental quantities that can be evaluated by MD/MC sampling. The phantom-wall method [48,49], dry-surface method [50], and interface wetting potential method [51,52] are combined with MD simulations; the test-area perturbation approach [53,54] and the excess surface free energy approach [55][56][57] employ MC simulations.
As explained later, MD/MC requires molecular force fields appropriate for individual metallic surfaces; therefore, its applicability strongly depends on available force fields.
(II) The direct approach usually employs MD-based simulations to construct molecular droplets and attempts to gauge the contact angles from the MD snapshots of the droplet geometries. Several schemes [44,[58][59][60] have been proposed to identify the "droplet shape" from the molecular condensation; then, the contact angle is evaluated via elemental differential geometry. Note that this direct approach (II) involves more uncertainty than the indirect one (I), but the direct one is applicable even for determining wettability on rough surfaces and dynamic wettability. This implies that (i) the indirect approach provides more accurate (static) contact angle estimations than the direct one, assuming the interface energies to be accurate [46], and (ii) the direct approach provides more general estimations applicable to a wide range of wettability phenomena [44].
In both the direct and indirect approaches, the accurate prediction of the contact angle strongly depend on that of molecular interactions at each interface. Molecular force fields used in MD/MC simulations are crucial for reproducing the molecular interactions [61], whereas ab initio DFT approaches usually predict the energies accurately. In MD simulations, Lennard-Jones (LJ) potentials associated with electrostatic ones are mostly selected as the force field. [62] The LJ potentials require a set of parameters for each molecular interaction; major parameter sets are available for water [62,63], but individual parameter fittings are necessary for different surfaces by comparing with ab initio simulations or experiments. [64,65] MD simulations with various types of the force fields have been applied to evaluate water contact angles on metallic surfaces such as Pt [66] and face-centered cubic Cu [67], but the accuracy of their predictions has been reported to vary depending on the force fields adopted. [59] Recently, several attempts have been made to address the force field issue in classical MD.
A straightforward way is to directly improve the force field accuracy by fitting DFT data, [68] although its prediction should be calibrated by the experimental validation. [59] Further, ab initio MD (AIMD) has been applied to surface/interface properties of water and evaluated the surface tension [69] and the contact angle [70]. The AIMD overcomes (semi)empirical force fields, but has a remarkable limitation on system sizes; hundreds of water molecules in a previous study [70] were insufficient to identify the "droplet shape," as more than 500k water molecules were required for that purpose [60]. Within the framework of the direct approach, ab initio evaluations of the contact angle are still a challenge for AIMD, whereas classical MD involves a number of caseby-case parameter fittings. From the viewpoint of complexity of simulation procedures, our DFTbased indirect approach can be regarded as being simpler than the MD-based direct approach.
Finally, we consider the temperature dependence of the contact angles. For the water wetting, its contact angles at ambient/low temperatures (less than a boiling point) remain almost constant, whereas those at high temperatures suddenly decrease. [71] For example, the contact angles of water on an aluminum surface are independent of the temperature at 20 < T < 120 • C being 90 • , while the temperature dependence is given as 157.4 − 0.55T at T ≥ 120 • C. Looking at our scheme, we took the Boltzmann averaging over several types of water clusters to evaluate the contact angle, which might be expected to describe the temperature dependence. However, our scheme provides almost constant values of the contact angles within 0.2 • for a wide range of temperatures 0 < T < 150 • C. These values are consistent with those at low temperatures, whereas they are inconsistent with those at high temperatures. This is because no evaporation is assumed in our scheme, while in reality, water evaporates at high temperatures, thereby reducing the surface tension as the temperature increases. [63,72] Because our scheme assumes that there is no intake/outtake of water molecules, the applicability of our scheme is restricted to low/ambient temperature regions. We thus conclude that our scheme would successfully predict the contact angles at low temperatures.

D. Limitations in comparison with reality
We note that the comparison with experiments shown in Table II should be made only for the limited purpose of demonstrating the superiority of our structural model over the conventional honeycomb model. [13,14] Though the present model achieves the maximum possible toward the reality within the available computational feasibility, there are several unsatisfactory factors to describe the reality as we can point out below. To obtain a quantitatively reliable estimate for the interaction energy between water and the surface, one would wonder that the larger coverages than one bilayer as well as structural disorder in the water layer should be considered. To capture such factors, the larger simulation cells are required. The feasibility of practical ab initio calculations is, however, subject to limitations on simulation sizes. The computational cost typically scales as the cube of the system size. [73] The available memory capacity is far smaller than that required to store the dynamic updates of such a larger system.
Not only in the present model, but also in all the current ab initio treatments, γ LG is approximately modelled by the surface energy of Ih-ice. [1,11] Though there are several possible criticisms against this approximations, this modelling is actually the keystone for the treatment, otherwise it becomes infeasible to describe the contact angle in the ab initio framework within the practical computational cost. At least within the current feasibility, we can either (a) compromise between the limitation of the modeling and the universality of ab initio framework to capture the trends with the same formalism, or (b) abandon the universality, just focusing on a specific system to pursue the descriptiveness of reality by using qualitative model-theoretical treatments.
As we discussed in Table III, the dependence of the predictions on the exchange-correlation functionals is another point for further investigations. The functionals capable of capturing the dispersion interactions have gradually be underway, but they are typically 10 times heavier than the conventional ones in their computational cost as a result of the extra processing (e.g., nonlocality) required. [74] Hence, from the viewpoint of the computational cost, it is a trade-off between the direction toward capturing the disorder effects and that toward describing the dispersion interactions in detail.
It should be noted, again, that the possible criticism on the incompleteness of the modeling enumerated above is not specific to the present modeling, but is commonly applied to all the previous ab initio studies. Although the common criticism still remains unsolved, the superiority over the conventional geometrical modeling has been performed for the metallic surfaces within the available computational resources at the best effort, as emphasized again.

IV. CONCLUSIONS
In this study, an ab initio evaluation of the water contact angle on Cu(111) surface was performed by determining the surface energy of water (∼ γ LG ) using that of ice crystal, which has been widely used in previous studies. [1,11] The results revealed that the isolated water oligomers used in this study as the adsorption model afforded a more accurate prediction of the water contact angle on the metallic Cu(111) surface than that obtained by the periodic-honeycomb lattice of ice crystals, which has been widely used for estimating the water contact angle on the surfaces of insulating materials. [1,11] The results were supported by several experimental observations. [14,16,21] For the freedom of the model choice (namely, the size of the oligomers and the coverage), we revealed that averaging out the possible choices using the Boltzmann weight based on the adsorption energies leads to a fairly reasonable prediction compared with the experimental values.