Ab-initio-based Interface Modeling and Statistical Analysis for Estimate 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, which involves developing ab initio prediction scheme of accurately predicting the angle. The scheme requires structural models for the adsorption of liquid molecules on a surface, but their reliability depend on whether the surfaces comprise insulating or metallic materials. Previous ab initio studies have focused on the estimation of the water contact angle on insulators, where the periodic-honeycomb array of water molecules was adopted as the adsorption model for the water on the insulating surface and succeeded in the insulating cases. This study, however, focus on the water contact angle on a metallic surface, and propose a simple ab initio based estimation scheme. We not only adopt the previously proposed structural modeling based on the periodic-honyecomb array, but also consider an ensemble of isolated water oligomers that have different molecular coverage (ML) values. We established a statistic model to predict a contact angle of the water wetting on a Cu(111) surface: The coverage-dependent contact angles obtained from each of the isolated clusters was fit to a quadratic regression, and the contact angle was interpolated by referring to a ML value of water layer in literature. This interpolated value lay within the deviation of experimental angles. In addition, the Boltzmann-average over the isolated clusters was found to agree well with the interpolated one. This indicates that the Boltzmann-average is useful for estimating the contact angle of other metallic surfaces without knowing a ML value a priori.

surfaces comprise insulating or metallic materials. Previous ab initio studies have focused on the estimation of the water contact angle on insulators, where the periodichoneycomb array of water molecules was adopted as the adsorption model for the water on the insulating surface and succeeded in the insulating cases. This study, however, focus on the water contact angle on a metallic surface, and propose a simple ab initio based estimation scheme. We not only adopt the previously proposed structural modeling based on the periodic-honyecomb array, but also consider an ensemble of isolated water oligomers that have different molecular coverage (ML) values. We established a statistic model to predict a contact angle of the water wetting on a Cu(111) surface: The coverage-dependent contact angles obtained from each of the isolated clusters was fit to a quadratic regression, and the contact angle was interpolated by referring to a ML value of water layer in literature. This interpolated value lay within the deviation of experimental angles. In addition, the Boltzmann-average over the isolated clusters was found to agree well with the interpolated one. This indicates that the Boltzmann- * Correspoding authors Email addresses: tak.murono@icloud.com (Takahiro Murono), kenta_hongo@mac.com (Kenta Hongo)

Introduction
Regulating the wettability of surfaces [1,2,3,4] is an important issue that must be addressed to broaden their industrial applications, such as heterogeneous catalysis, corrosion, and electrochemistry. [2,5,6] The wettability of the surfaces of metal catalysts such as Cu has a significant effect on their catalytic activity for electrochemical reactions, [7] 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. [8,9] The angle is determined using the Young's relation [10,11]: where γ SG , γ LG , and γ SL are the surface energies at solid-gas, liquid-gas, and solidliquid 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,2] For example, the contact angle of water on Si [1] 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 [2]. For example, the water contact angle More sophisticated ab initio models based on quantum molecular dynamics (QMD) have been proposed and applied to several insulating systems such as CaCO 3 [13,14] and graphite [15]; for CaCO 3  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, [2,1] an "ice-like bilayer model" [16,17] (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)). Note that the experimental contact angles for the insulating surfaces are definitive due to their small errors and hence the ice-like bilayer modeling can be validated by comparing the resultant values with the experiments.
Similarly, the QMD-based approaches can be validated, even though their computation is heavily demanded.
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 [18,19], implying that the "bilayer model" would not be necessarily appropriate for the ab initio evaluation of the water contact angle on the metal surface. Instead, "isolated models" that do not exhibit the bridging of hydrogen bonds over neighboring unit cells would be rather appropriate.
In addition, experimental values of the water contact angle on the Cu surface deviate from 60 • to 85 • , depending on various experimental conditions that are not completely controllable. [20,21,22] This implies that unlike the insulating surface, the validation of theoretical works on the water-copper wettability is difficult as mentioned before, but only the range of the contact angles can be addressed.
In this study, we aim at establishing the ab-initio-based scheme of estimating the water contact angle on a metallic copper surface, i.e., the metallic Cu(111) surface. We first construct several isolated models (cluster/oligomer) and then peform the ab initio adsorption energy evaluation of the cluster models as well as that of the bilayer model.
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 coverage-dependent contact angles, which is equivalent to the dependence of the adsorption energy on the size/coverage of the oligomers. Referring to a reported coverage, an interpolation between the oligomer and periodic models leads to our ab initio prediction of the water-copper contact angle. We found that (i) the coverage-dependent contact angles monotonously decreases as the coverage increases, (ii) deviation of the coverage dependence is less than that of the experimental angles (25 • ), and hence (iii) the predicted contact angle lies between the experimental deviation. For comparison, we considered the Boltzmann average over the individual contact angles, which was found to be similar to the contact angle predicted from the coverage-dependent angles.
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 the numerical results obtained from our approach, 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.

Interface modeling
Based on the previous studies, [1,2] we estimated 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,2] 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 [2]: 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. [23] For the ease of comparison with previous works, [1,2] we used the same exchange-correlation function, GGA-PBE, [24] used in previous studies. Because the valence shell of Cu consists of d-electrons, the Hubbard U correction [25,26,27] may also be important.
However, this correction matters when the orbital becomes more localized where the self-interaction cancellation is damaged. [28] For our metallic Cu surface, the orbital is delocalized, and hence, we did not consider such corrections, as in most previous studies. [29] In addition, we did not consider the exchange-correlation functional with van der Waals (vdW) corrections due to Tkatchenko and Scheffler [30], because the vdW-DFT is know to fail to reproduce the delocalized nature of electrons. [13,31] Accurate reproduction of the non-covalent molecular interactions involves other approaches such as quantum Monte Carlo (QMC), [32,33] but QMC is not applicable to our target system including water clusters on the delocalized metallic surface due to its much heavier computation than DFT. Norm-conserving pseudopotentials [34] were used to describe the ionic cores. The detailed computational conditions for each calculation is summarized in Table 1.
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 periodichoneycomb model (bilayer) [16,17] [panel (a) in Fig. 1] and those by the isolated molecular models (buckled) [19,35] [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 [36] (planar model) 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 [35,36] [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. [19,35] The dependence on the coverage, i.e., on the number of N-mers placed within a unit cell, was examined up to N = 3.

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 Cu-struct-opt 8×8×8 800 Cu-surf-relax 11×11×1 800 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 struc-  Fig. 1], we used the "bilayer model" employed in preceding studies. [16,17] The structure was generated from the "H-parallel" structure, [36] 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. [19,35] To evaluate the contact angles in the four-layer slab models (as shown later in Fig. 2), water molecules were placed on the models using the on-top alignment geometry [35,36] with the H-parallel orientation [36]; herein, the vertical height was optimized by DFT.   It is found from Fig. 2 that the coverage-dependent contact angles monotonously decreases as the coverage increases. In order to determine the computed contact angle from the dependence, we need to determine an appropriate value of coverage for water layer on Cu(111). The knowledge of the coverage enables us to estimate the contact angle from the coverage-dependence of the angle. In the present study, we adopt the coverage of ∼ 0.4 ML that is estimated from a previous study [40] based on a Langmuir adsorption isotherm. [41] Referring to the fitted curve, an interpolation into the 0.4 ML leads to the contact angle of 74 ± 6 • . As can be seen from Table 2 In addition, we evaluate the Boltzmann-averaged contact angle over the isolated cluster models, leading to its contact angle of 78 ± 5 • . The Boltzmann-averaged angle agrees well with both the isolated-cluster-based and interpolated ones, though the former is slightly larger than the latter. The Boltzmann-average scheme is a straight-  [36]) were larger than those obtained in this study. The significant increase in values for the dimer can be attributed to the difference in definition. In their definition, the hydrogen bonding interactions between water molecules were included in the E watar ads , and hence, it increased with an increase in the size of N-mer.

Water coverage on metallic surfaces
Here we make some remarks on the water coverage that is a key quantity of estimating the contact value from the coverage-dependence, though it is not necessary for the Boltzmann-average-based estimation. As a matter of fact, the Cu(111) surface has a lower coverage than other metallic surfaces, e.g., Pd(111), Pt(111), and Ru(0001). [42,43]  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. [17] The lower coverage for Cu(111) was derived in a previous 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. Because the Langmuir adsorption isotherm [41] is basically applicable to solid-gas (steam) interfaces, applying it to a solid-liquid interface requires careful consideration.
Applying it in this case would be justified when (i) the intermolecular interactions between water molecules are negligible, and (ii) the absorption process of the water molecules on the solid surface is equilibrated and then none of the absorbed molecules are dissociated into H + and OH − , so the molecules in the liquid water can then be regarded as molecules in steam. The isolated model satisfies (i), and the resultant lower coverage seems self-consistent. For (ii), the dissociation energy of the water molecules on Cu(111) has been reported to be 1.23 [eV/H 2 O], [44] being much larger than the energy scale of the adsorption energy, ∼352 [meV/H 2 O]. [36] 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.

Comparison with other methodologies
Herein, we compare the computational evaluations of the contact angle obtained herein with those obtained previously because it will be helpful for deeper understand- The ab initio simulations with interface modeling appropriate for insulating surfaces ("ice-like bilayer" model) have proven to accurately reproduce the contact angles. [1,2] Recently, a more sophisticated approach based on ab initio molecular dynamics (AIMD) simulations has been proposed by Lu et al. and applied to insulators. [13,14,15] The AIMD approach takes into account relaxation of water molecules and then agree well with experiments. This approach is, however, infeasible for metallic surfaces: The AIMD simulation is done on the Γ-point, so its simulation cell is larger than that of usual band-structure simulations in order to treat condensed matters.
Because of this requirement for a large supercell, the AIMD approach is inadequate for our delocalized Cu surface that requires a much larger supercell than insulating surfaces. Instead, as explained before, the present study proposed "isolated cluster" models, which relies on the experimental fact about a lower coverage on the Cu surfaces. [48] In addition, experimental contact angles on the Cu surface deviate widely, depending on their sample preparations and so on. This indicates that we cannot make a simple comparison between experiment and theoretical contact angles, and hence validate reliability of our proposed interface modelling. It is to be noted, however, that deviation of the coverage-dependent contact angles is less than that of experimental angles.
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 [49,50], dry-surface method [51], and interface wetting potential method [52,53] are combined with MD simulations; the test-area perturbation approach [54,55] and the excess surface free energy approach [56,57,58] 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 [59,60,45,61] 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 [47], and (ii) the direct approach provides more general estimations applicable to a wide range of wettability phenomena [45].
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 [62], 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. [63] The LJ potentials require a set of parameters for each molecular interaction; major parameter sets are available for water [63,64], but individual parameter fittings are necessary for different surfaces by comparing with ab initio simulations or experiments. [65,66] MD simulations with various types of the force fields have been applied to evaluate water contact angles on metallic surfaces such as Pt [67] and face-centered cubic Cu [68], but the accuracy of their predictions has been reported to vary depending on the force fields adopted. [60] 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, [69] although its prediction should be calibrated by the experimental validation. [60] Further, ab initio MD (AIMD) has been applied to surface/interface properties of water and evaluated the surface tension [70] and the contact angle [71].
The AIMD overcomes (semi)empirical force fields, but has a remarkable limitation on system sizes; hundreds of water molecules in a previous study [71] were insufficient to identify the "droplet shape," as more than 500k water molecules were required for that purpose [61]. 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 case-by-case parameter fittings. From the viewpoint of complexity of simulation procedures, our DFT-based indirect approach can be regarded as being simpler than the MD-based direct approach.

Limitations in comparison with reality
The present model achieves the maximum possible toward the reality within the available computational feasibility, but 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. [72] 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,2] 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 mentioned in section "Computational Methods", dispersion or vdW interactions is another point for further investigations. The exchange-correlation functionals capable of capturing the dispersion interactions have gradually be underway, but they are typically a number of times heavier than the conventional ones in their computational cost as a result of the extra processing (e.g., non-locality) required. [73] 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. In addition, the current vdW-DF approach is not appropriate for delocalized systems, [13,31] implying other approaches beyond DFT are necessary for the present case.

Conclusion
In this study, we established two schemes for estimating the water contact angle on the Cu (111) surface; one is a regression-based interpolation of the coverage-dependent contact angles at the reference value of the coverage, the other is the Boltzmannaveraged contact angle over all the isolated clusters with different values of the coverage. The Boltzmann-averaged estimate was found to agree with the regression-based one within their error bars. We found that experimental contact angles in literature remain ambiguous, and the uncertainty of experimental contact angles is larger than our two models' error bars. Comparing our two models, the regression-based estimate involves the reference molecular coverage a priori, whereas the Boltzmann-averaged estimate is more straightforward than the regression-based one. This indicates the Boltzmann-averaged estimate would be applicable to other metallic surfaces without knowing a molecular coverage value a priori.