Water on hydroxylated silica surfaces: Work of adhesion, interfacial entropy, and droplet wetting

In the last few years, much attention has been devoted to the control of the wettability properties of surfaces modified with functional groups. Molecular dynamics (MD) simulation is one of the powerful tools for microscopic analysis providing visual images and mean geometrical shapes of the contact line, e

The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp with silanol (≡SiOH) groups. It is known that increasing the surface OH density typically results in the enhancement of the wettability. 2 Despite widespread usage of silica for a long time, some properties related to silica when interacting with water are not completely understood. More concretely, it is not straightforward to quantitatively predict the level of wettability enhancement, as silica presents different morphologies (crystalline and amorphous) and polymorphisms (cristobalite, quartz, etc.) and as also the silanol groups can present different conformations (isolated, geminal, and vicinal silanols). 3,4 To determine the surface OH structure and density, infrared spectroscopy, [5][6][7] thermo-gravimetric analysis, [8][9][10][11][12] chemical methods, 5,8 and nuclear magnetic resonance 9,13,14 are used, while the straight link between the OH density and water wettability has not been provided because water molecules adsorbed on the silica surface, i.e., not in the liquid phase but under saturated vapor condition, are measured in these experiments. In addition, there is experimental evidence that the silanol group distribution on the silica surfaces is not homogeneous, but rather in patches, 15 while the effects of such inhomogeniety on the wetting behavior are not clear.
Altogether, these characteristics make it challenging to propose a general model that describes the interaction between silica and water. Given the nature of the system, with the presence of OH groups interacting with water molecules, continuum approaches are usually not enough, and studies at the atomistic level are needed. In this regard, molecular dynamics (MD) simulations have been proven very useful and have been applied to study several phenomena, from biology to materials science, and the study of silica wettability is no exception. The first force field for studying silica properties in the bulk and surface started to be developed already in the 1980s, just after the beginnings of the MD method. [16][17][18] In the past two decades, we have seen a substantial number of studies, considering, for instance, layers of water on top of silica 19 or water molecules confined by layers of silica, 20,21 generally presenting crystalline silica structures created from the beta-cristobalite 19,22,23 or amorphous silica models. 24 Even ab initio calculations have been considered, 25 but the computational cost involved precludes the study of large systems. The properties of the silica-water interface have been discussed in terms of the behavior of the water molecules at the interface, the network of hydrogen bonding, or the effects on the water far from the interface, 19,22,23 but little consideration has been put on obtaining a precise estimation of the water contact angle itself. In addition, many studies were performed at room temperature, in which alpha-cristobalite silica is known to be the most stable form, but the beta-cristobalite silica model was employed, [19][20][21][22][23] due to the early availability of force fields.
More recently, Emami et al. 26 developed a set of very useful force field parameters for simulations of silica, part of the Interface Force Field (IFF) database, 27 addressing many of the limitations observed in early models. They have used several spectroscopic data to determine the bonded parameters (x-ray measurements for bond distance and bond angles and IR/Raman for bond stretching and angle bending constants) and refined the nonbonded parameters to reproduce hydration energies. They were able to reproduce several experimental properties of silica surfaces interacting with water, such as calorimetric heat of immersion, contact angles of surfaces with different silanol densities, and adsorption isotherms. Their parameters are compatible with several commonly used force fields, such as Assisted Model Building with Energy Refinement (AMBER) Chemistry at Harvard Molecular Mechanics (CHARMM), and Condensed-Phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS), and the database contains models for simulating silica surfaces composed of different environments, such as Q2 [geminal silanol groups, = Si(OH) 2 ], Q3 [single silanol groups, (≡SiOH)], and Q4 (no silanol groups, i.e., only siloxane bridges). It also has adjustable degrees of ionization, allowing the simulation of different pH conditions. 28 Two of the authors have applied this force field in a previous study, 29 where the thermodynamic integration (TI) implemented by the dry-surface method was used to calculate the work of adhesion, i.e., the solid-liquid (SL) interfacial free energy relative to the liquid-vapor (LV) and solid-vapor (SV), of water on nonhydroxylated silica surfaces to examine the effects of Coulomb and van der Waals interactions on the SL interfacial energy, where the entropy effects on the SL interfacial energy were also shown. Contemporaneously, Schrader et al. 30 combined experimental and computational efforts to study how heterogeneity can change silica surface hydration properties. They have employed the IFF in their MD simulations, where hydroxylated silica surfaces were generated by using genetic algorithms, and observed that, for water diffusivity, result values can change up to 10% depending on the OH distribution on the silica surface.
From a point of view of fundamental wetting physics unrestricted to silica-water systems, correspondence between the theoretical (or ideal) contact angle calculated by Young's equation and interfacial tensions and the apparent contact angle of the droplet has been discussed. 31 More specifically, the interfacial tensions for the former were obtained either by using thermodynamic methods including the TI, called the thermodynamic route, [32][33][34][35][36][37][38][39] or by mechanical methods obtained from the integration of the stress at the solid-liquid interface, called the mechanical route, [39][40][41][42][43][44] whereas for the latter, the apparent contact angle was typically measured from the average density distribution. 45,46 Recent results showed that the apparent and theoretical contact angles agreed in systems on smooth or relatively homogeneous surfaces, pinning force or line tension should be included in the force balance of non perfectly flat or inhomogeneous surfaces. 38,39,44,[47][48][49][50][51][52][53][54] In order to deepen this discussion about the agreement between the results of the two routes, in the present work, we aimed at extending our previous results, performing MD analysis on the wetting of water on hydroxylated silica surfaces. In particular, we aimed at investigating the effect of the OH groups on the SL interfacial energy by obtaining the free energy through TI. Work of adhesion and entropic effects were analyzed and discussed as well. Then, the contact angle of water interacting with the hydroxylated silica surfaces was calculated and compared with the apparent contact angle obtained from MD simulations of water droplets interacting with silica. To focus only on the effect of the OH density, we employed a neutral and regular silica surface model, where no annealing process was used, i.e., undeformed surfaces, which serve as a first approach, but resulted in a less wettable silica surface than the original model by Emami et al. 26 In the case of TI calculations, here implemented by the phantom wall scheme, a series of simulations of a water film in contact with the silica surfaces is performed, and a repulsive potential wall interacting only with the water film (the phantom wall) is The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp introduced. By displacing the position of this phantom wall, the water film from the surface is lifted quasi-statically. This allows the calculation of the work of adhesion, and if the liquid-vapor interfacial tension is used, it is possible, from the Young-Dupré equation, to evaluate the contact angle. As for the case of the MD simulation of droplets, a larger surface is created and the water film interacting with it will eventually take a droplet shape if the spreading coefficient is negative. 55 Then, by estimating (or, more commonly, by assuming) the position of a contact plane, it is possible to estimate the contact angle. The TI calculations have the advantage of being free of arbitrary choices such as contact plane position or circle fitting procedure, two of the main disadvantages of the droplet simulation; the surface employed in TI calculations can also be smaller, and this contributes to decrease the computational cost, which can be substantial, as a series of MD has to be performed. From TI calculations, however, it is not possible to have information about the dynamics of the droplet on the surface, i.e., Brownian motion and pinning effect cannot be studied. This is a natural advantage of the MD droplet simulations, which allows an analysis of the temporal evolution of the droplet on the surface. Despite their differences, we observed that the results provided by both methods are very close, and although the limitations imposed by our choices caused an overall overestimation of all values, the expected trend is preserved, with the silica surface becoming more wettable as the area density of OH groups increases. We also discussed some physical aspects related to some deviations between the two methods, in terms of behavior of the water molecules close to the surface and the character of the chemical groups involved in the process.

II. METHOD
We performed MD simulations of water interacting with the silica surface using the LAMMPS package (August-2019 version). 56 In all the performed simulations, the velocity Verlet algorithm was used for time integration, with a time step of 0.5 fs. The visual molecular dynamics (VMD) software 57 was utilized to visualize the trajectories and draw the structures presented in Figs. 1-3 and 7.

A. Preparation of silica surfaces
The modeling of the silica surfaces very closely follows the procedure used in our previous study, 29 but for the sake of consistency, we present here the most relevant details. We used the parameters developed by Emami et al. 26 for the silica force field, which are part of the surface model database in the Interface force field (IFF). 27 These silica models are consistent with those of other pH-sensitive silicates such as tricalcium silicate, clay minerals, and other related aluminate and phosphate minerals, but here, we limited our study to only neutral hydroxylated silica surfaces. Harmonic potentials were used to describe the motions related to bond distances and bond angles.
The IFF provides parameters compatible with several commonly used force fields such as CHARMM and AMBER. In LAMMPS, a family of CHARMM-like pair styles has been implemented, but brings some inconveniences. For instance, "lj/charmm/coul/long" has an incorrect CHARMM implementation of force switching between the inner and the outer cutoff, while "lj/charmmfsw/coul/long" has the correct implementation, but the force is not strictly the derivative of the energy, and also, some packs to accelerate the performance of LAMMPS running in Intel-based central processing units (CPUs) are not available for this particular pair style. Thus, we adopted the usual Lennard-Jones (LJ) and Coulomb interactions as implemented by the more common "lj/cut/coul/long" pair style, which is simple and correct, and also presents accelerated performance in Intel-based CPUs, but the difference between the force fields demands a reparameterization of the unit cell to reduce the residual stress in the structure. For this reparameterization, we used silica bulk simulations to redetermine the unit cell parameters. We started with the αcristobalite structure provided in the supplementary material of Ref. 26, composed of 12 silicon and 24 oxygen atoms. The structure presents the (101) plane perpendicular to the z direction and has the following dimensions: 0.855 × 0.498 × 1.214 nm 3 . We adopted this structure as a basic unit cell and replicated it in the three directions, 4 × 7 × 3 basic cells and 6 × 10 × 4 basic cells, obtaining two bulk structures (labeled small and large, respectively) for the silica. Then, we performed simulations in the NPT ensemble, considering two thermodynamic conditions: room condition at T = 300 K under ambient pressure (1 atm) and the same temperature under zero pressure. These thermodynamic conditions were controlled using Nosé-Hoover style equations, with the thermostat and barostat having three chains each, and damping coefficients of 50 and 500 fs, respectively. In addition, due to the anisotopic nature of the cell, the dimensions in the three directions were controlled independently (aniso option). In all cases, periodic boundary conditions were applied in all directions. The system was thermalized for 0.5 ns, and additional 0.5 ns were used for taking averages. The results presented very similar dimensions, as shown in Table III of Appendix A, and were consistent with our previous study. 29 After the redetermination of the unit cell dimensions using simulations for the bulk case, we move to the construction of the initial silica surfaces. In order to create the initial surface used for the TI calculations, we started creating a bulk structure replicating the α-cristobalite structure (4 × 7 × 2 basic cells). Then, the chemical bonds crossing the z direction were broken. Finally, oxygen atoms were placed to construct the siloxane bridges. This system was used to perform a simulation in the NVT ensemble, at 300 K, where the Langevin thermostat, with a damping coefficient of 1 ps, was employed. Here, periodic boundary conditions were employed in the x and y directions, while the mirror condition was applied in the z direction. The system was thermalized for 1 ns, and an additional run of 1 ns, which outputs configurations every 200 steps, was used to produce the data used to calculate the silica atom mean positions. By using these calculated mean positions, we created an average silica structure to be used as the initial structure in the simulations of silica-water surfaces. As we will detail below, the atoms located at the bottom part of the silica surface were kept fixed in their average position during the simulations. No annealing process was employed, so the structures are very regular ( Fig. 1) and, consequently, artificial to some extent. This does affect the wettability, decreasing it, but allows us to study the hydroxylation effects on the work of adhesion and on the entropy.
In our previous work, 29 we observed an expansion of 11% of the silica crystal structure in the z direction, when compared with the bulk silica. In the present work, we do not observe this, but rather a contraction of 10%. This difference is due to different treatments of the 1-4 interactions between bonded atoms. In our previous study, these interactions were excluded, but they have been included in the present study, following the original model. 26 To hydroxylate the silica surface, we removed the oxygen atom of the siloxane bridges and doped the respective silicon atoms with OH groups, producing pairs of silanol groups. In the model of Emami et al., 26,27 the total charge of each OH group is half of that of removed oxygen atoms, as shown in Table I, so this procedure preserves the neutrality of the system and avoids introducing counterions into the simulation box. This also means that the silica surfaces have a zero ionization level, which would correspond to a situation of acid pH.
The pairs of OH groups present, at the closest, an O-O distance of ∼0.4 nm, so there is no hydrogen bond formation among them, and they can be considered as isolated silanol groups rather than vicinal ones. We consider four levels of hydroxylation. In the case of fully hydroxylated silica surface, all O atoms of siloxane bridges at the surface were replaced by OH groups (i.e., ρ OH A = 4.700 nm −2 ). For the intermediate cases, we started with the fully hydroxylated surface and randomly deleted pairs of OH groups and replaced them with an oxygen atom at the suitable position, restoring the siloxane bridge. This procedure produced OH distributions without any symmetric or homogeneous appearance, which can be regarded as more natural, as there is some experimental evidence that the silanol groups do not appear homogeneously distributed over the silica surface, but rather in patches. 15 Three intermediate cases were considered as partially hydroxylated surfaces, where 75%, 50%, and 25% of the siloxane bridges were replaced by silanol groups, resulting in ρ OH A of 3.525, 2.350, and 1.175 nm −2 , respectively. The snapshots of the partially hydroxylated surfaces are shown in Fig. 1.
These surfaces were used for the TI calculations, whereas for the water droplet simulations, we replicated the surface eight times in the x direction, producing a wider system with an x-width of about 30.8 nm, as shown in Fig. 2. We opted to extend the surface through replication instead of randomly generating a new OH distribution in order to keep the consistency between both methods, aware of the potential change in the results caused by differences in the OH distribution. 30 The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp

B. Interaction potentials
As in our previous study, 29 we employed the SPC/Fw model for water, 58 as several properties calculated with this model present good agreement with experimental data. 59 Harmonic potentials were used to keep both, bond distance and bond angles, allowing changes in the dipole moment. This treatment can be regarded as more realistic, especially when dealing with polar surfaces. For all silica-water systems considered in this study, the liquid film was composed of 3600 water molecules.
We used the Lennard-Jones (LJ) and Coulomb interactions as the potential functions. The combination rules of Lorentz-Berthelot were employed. For the van der Waals interactions, a strict cutoff of 1.2 nm was used. For the Coulomb potential, the particle-particle particle-mesh (PPPM) method is applied to treat long-range electrostatic interactions. Relative accuracy is set to be 10 −6 , and a strict cutoff of 1.3 nm was used. As the system is non-periodic in the z direction, the slab option implemented in LAMMPS was used in this direction in order to prevent the appearance of spurious interactions between the real system and its replicas. 60 The system temperature was maintained at 300 K. The silica temperature was controlled by using a Langevin thermostat, 61 while the liquid film was controlled by the thermostat proposed by Bussi et al., 62 as it has been shown that it correctly reproduces the water diffusive behavior independent of the damping coefficient. 63 For both thermostats, a damping coefficient of 1 ps was used. The Langevin thermostat was used only in an intermediate region of silica surface, comprising 0.8 nm of height, starting about 0.5 nm away from the top silica layer, as shown in Fig. 2(a).
One additional standard liquid film system, with dimensions of 3 × 3 × 18 nm 3 and consisting of 1800 water molecules, was created without any silica surface to measure the liquid-vapor interfacial tension, which can be calculated from the components of the pressure tensor, where the pairwise interaction form of the virial pressure is used. 64 Other simulation conditions are the same as those of the film systems. This system was equilibrated for 1.0 ns, and the time average of the following 30 ns was used to extract the pressure tensor. The resulting liquid-vapor interfacial tension γ LV was (59.29 ± 0.37) × 10 −3 N/m. This value was used in the Young-Dupré equation to obtain the contact angle from the TI calculation as mentioned below.

C. Thermodynamic integration through the phantom wall method
In general, the TI method is applied to obtain the difference in the free energy between two equilibrium states of a certain system, by introducing a parameter variable that is analytically differentiable in the Hamiltonian H and that connects both states through a quasi-static path. More specifically, if this parameter is λ, then the Hamiltonian is given by H(Γ, λ) as a function of the phase-space variable Γ consisting of all positions and momenta of the constituent molecules and λ. If we set that λ = 0 and 1 are, respectively, the reference and target systems in the NVT ensemble, the difference ΔF in the Helmholtz free energy F between them is given by the following statistical mechanical relationship: One should note, however, that the brackets, within the statistical mechanical definition, correspond to the ensemble average, while in an MD simulation, this is replaced by the time average, and by performing the calculation for many independent equilibrium systems between λ = 0 and 1, we can numerically integrate the expression. Actually, in most cases, including this study, it is more convenient to use the reversible path, going from the target system to the reference one. In this case, formally, the integration is performed in the most right-hand side of Eq. (1). In this study, a phantom wall, 33 parallel to the solid-liquid interface, was introduced, and its vertical position was expressed via λ. If the phantom wall only interacts with the liquid film molecules, and the interaction is short-range, then the λ = 1 state can be obtained by having the phantom wall sufficiently below the liquid film. On the other hand, the phantom wall can also be quasi-statically moved upward in an isothermal process to a position where the interaction between the solid and the liquid film becomes very close to zero. Due to the long-range Coulomb interaction, it is not possible to obtain a zero value, so we assumed the state to be λ = 0 when very small values for the force exerted by the wall were obtained (between 0.2 and 0.4 ×10 7 N/m 2 ). In this study, the vertical positions of the phantom wall in the reference system and the target system are defined as z 0 and z 1 , respectively, and the position of the phantom wall zw(λ) is expressed by the formula The interaction between the virtual wall and water can be set arbitrarily if the above conditions are satisfied, but numerical integration is facilitated by applying only repulsive forces.

ARTICLE scitation.org/journal/jcp
Thus, in order to perform the numerical integration, we performed a series of MD simulations of the system represented at the left side of Fig. 2. The phantom wall started 0.1 nm above the region being thermostated, i.e., at least 0.3 nm distant from the water liquid film for the non-hydroxylated silica surface, and more than 0.4 nm for the hydroxylated silica surfaces (Fig. 2). An MD simulation in the NVT ensemble is performed with 0.1 ns of thermalization and 0.5 ns of sampling. The wall is then displaced by 0.01 nm in the z direction, and the procedure is repeated. This is done for a total displacement of 0.7 nm, i.e., 71 simulations were performed for each of the five surfaces, when it became clear that the interaction between water molecules and the silica surface was sufficiently small. This small displacement (and, consequently, large number of simulations) is necessary to obtain smooth curves of the force exerted by the wall as a function of z, which will be integrated to calculate the work of adhesion. We remove the linear momentum from water molecules every 1000 steps in order to prevent the water liquid film from drifting away from the silica surface. The interaction between the liquid film and the phantom wall is implemented through the usual 12-6 Lennard-Jones potential. The phantom wall has parameters with values similar to those used to the water oxygen atom, i.e., σ = 0.316 nm and ε = 0.16 kcal/mol. In order to make it repulsiveonly, the cutoff is set to 0.3546 nm, so there are no attractive forces, and the phantom wall is set to interact only with the water oxygen atoms.
The obtained curves are then integrated by using the trapezoidal rule, resulting in the value for ΔF. Next, we use the fact that the work of adhesion (W adh ) is −ΔF per unit of area, i.e., In our case, we assumed the size of the area A as being the same as that of the xy-plane, i.e., 11.91 nm 2 (see Fig. 2). Finally, once we have W adh , we can apply the Young-Dupré equation, and obtain the contact angle θ YD .

D. Droplet simulations
We have also performed simulations of water droplets on silica surfaces to obtain the apparent contact angle θapp. As explained in Sec. II B, the silica surfaces employed in these simulations were the same surfaces used for the TI calculations, but replicated four times on both +x and −x directions (Fig. 2, bottom). The same initial configuration of water molecules used for the TI calculations was employed. Due to the system dimension, this initial configuration becomes a quasi-2D droplet, with the cross section of a cylinder, eliminating the effects from the line tension on the wettability. To obtain accurate values for θapp, typically long simulations are necessary. We performed 1 ns of thermalization followed by 10 ns of sampling. The results for the contact angle value were checked to be converged by analyzing the time evolution of the average value (see Appendix B). To facilitate the comparison of the results, the thermodynamic conditions are, as much as possible, the same as those used for the TI calculations. The simulations were performed using the NVT ensemble, at 300 K. Similar to the TI calculation, the Langevin thermostat was used only in the same 0.8 nm long intermediate region of the silica surfaces. However, for the water molecules comprising the droplet, no thermostat was used, in order to avoid any influence on the droplet shape. The trajectory generated by 10 ns simulation was used for obtaining the mean density distribution of the water. Using this distribution, it is possible to fit a circle to an isodensity contour, and by knowing (or, more commonly, by assuming) the position of the contact line, to obtain θapp.

A. Work of adhesion and entropy loss
We start discussing the results for the thermodynamic integration. Figure 3 presents the curves that describe how the force exerted by the wall changes with the wall position, obtained for the five cases. The curves are very smooth, as a result of choosing a sufficiently small value (0.01 nm) for the phantom wall displacement. As one would expect, as ρ OH A increases, the interaction of water molecules with the surface becomes larger, so the phantom wall must exert a larger force to lift the liquid film. This is very clear in Fig. 3: The curve for the non-hydroxylated silica surface (0 nm −2 , purple full circles) presents a peak around 30 × 10 7 N/m 2 ; the peak for the subsequent cases is displaced to larger values of z and presents a monotonic increase, with values of ∼35, 40, and 45 × Due to entropy-energy compensation in thermodynamic integration, 29,65 and in accordance to Eqs. (1) and (3), W adh can be expressed as where the subscript indicates the interaction between the phantom wall and the liquid (denoted by PL hereafter), and all the other components cancel out. For the PL interactions, u PL and TΔs PL are the corresponding mean energy and entropy per unit area. Note that the signs on the right-hand side depend on the integration direction. As we are more interested in the overall interfacial properties, than only just the PL interactions, we can freely add solid-liquid components back into Eq. (5), where u PL | λ=1 ≈ u PL | λ=0 = 0 was used, as the PL interactions at the reference and target systems are either non-existent or very weak (we estimated it as being on the order of 10 −5 mN/m for the non-hydroxylated surface, for instance). The mean solid-liquid energy per unit area −u SL is defined using the average of the total solid-liquid interaction potential energy Φ SL (Γ, λ) as The Note that the values of Φ SL (Γ, 1) and Φ SL (Γ, 0) can be calculated using the LAMMPS "group/group" command for the target and reference systems, and note also that −u SL is positive with this definition. On the other hand, the interfacial entropy difference Δs INT per unit area, consisting of both solid-liquid and PL components, is defined as Considering that TΔs INT is negative because the liquid at the SL interface with a restricted molecular motion at λ = 1 has lower entropy than the liquid at the interface with the phantom wall, interfacial entropy loss −TΔs INT given by will be used to describe the entropy component as opposed to the interfacial entropy gain to intuitively understand the entropy contribution to the relation between the work of adhesion and interfacial energy. The expression in Eq. (9) is identical to that of the dry-surface method used in our previous work, although the derivation is slightly different. 29 Figure 4 shows the work of adhesion W adh , average solid-liquid energy per unit area u SL , and entropy loss −TΔs INT for water on silica with various ρ OH A . As one can see, the three thermodynamic quantities increase linearly as ρ OH A increases, though all of them present different slopes. The linear relation between W adh and ρ OH A can be seen as a consequence of a small entropic effect, because the OH groups are flexible but short. A recent work has not observed such behavior when the surface is covered by long chains. 66 We have also noted that the entropy loss increases faster than the work of the adhesion, so the relation between the two quantities deserves a further discussion. Figure 5 displays the relations between W adh and u SL , as well as between −TΔs INT and u SL . Black lines and points represent the data obtained in the present work and, once again, show a linear increase in the W adh , as also of the entropy loss, with an increase in ρ OH A . In addition, the data obtained in our previous work 29 are also presented for comparison. Keep in mind that both, previous and present results, were obtained using thermodynamic integration, but in the present study, we used the phantom wall method, 33 The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp while the previous results were obtained using the dry-surface method. 35 In the dry surface, instead of introducing a phantom wall into the system and adopting its position as the parameter for the Hamiltonian, changes are made in the interaction potential between the silica surface and water molecules through the introduction of a set of parameters (κ, λ), which modify the Lennard-Jones (LJ) potential and Coulomb potential equations, as follows: In our previous study by strengthening the interactions, we observed an increase in the entropy loss, which was greater when increasing electrostatic interactions. Now, comparing this result with the present data, where the potential is fixed at the usual value, but ρ OH A has been changed, we can observe that the contribution of the OH groups to the entropy loss is large, making −TΔs INT increase even faster. A partially hydroxylated surface of ρ OH A = 1.175 nm −2 is already able to bring −TΔs INT to the same magnitude value of W adh , and at ρ OH A = 2.350 nm −2 , the −TΔs INT is already larger than W adh . This will ultimately influence the behavior of the water molecules near the surface, enhancing the wettability.

B. Water droplet simulation and apparent contact angle
As discussed in Sec. II B, we have obtained the interfacial tension γ LV = 59.29 × 10 −3 N/m for the SPC/Fw water model. Using this value and W adh in Fig. 4, we can obtain the contact angle θ YD through Eq. (4). On the other hand, from the simulations of the water droplet on different hydroxylated surfaces, we obtained the average density distribution around the center of mass (a discussion about the center of mass calculation is provided in Sec. III C) of the droplet, as shown in the top panels of Fig. 6. As one can see, as ρ OH A increases, the height of the droplet decreases, while the width at the bottom part of the droplet increases, which is a clear indication that the interaction between the surface and water becomes stronger as we add more OH groups, consequently decreasing the contact angle. In order to determine the apparent contact angle θapp, we fit a circle to a isodensity contour and measure the angle with respect to the contact plane at the interface. Based on the results presented on the top part of Fig. 6, we adopted a density value of 0.5 g/cm 3 for the isodensity contour.

FIG. 6.
Comparison between the contact angle values obtained from the thermodynamic integration calculations and the Young-Dupré equation (cos θ YD ), and the values obtained from the apparent shape of the droplet in MD simulations (cos θapp). Error bars for cos θ YD were obtained by performing the error propagation of γ LV , while error bars for cos θapp were calculated using the error from the circle fitting procedure. Top: the average density distribution of the water droplet for the five considered cases.

ARTICLE scitation.org/journal/jcp
For the contact plane, it is worth noting that, while for a more simplistic model, a definition for the plane is straightforward, for surfaces presenting roughness and (specially) inhomogeneity, as in our case, this is not necessarily trivial, and some compromise has to be made. Here, we adopt the contact plane as the following: for the non-hydroxylated surface, due to its hydrophobicity, we can observe the formation of a thin gap between the droplet and the surface. Considering the average position of O atoms at the topmost silica layer and the peak of the first adsorption layer, the distance between the liquid film and the surface is ∼0.26 nm. The contact plane should lie somewhere inside this gap. We decided positioning it in the middle, i.e., around 0.13 nm above the average position of O atoms at the topmost silica layer. Coincidentally, it would be roughly the position where the OH groups would be found if it was a hydroxylated surface. It is also around this position that the closest H atoms of the water molecules can be found. For the hydroxylated surfaces, the definition presents a more sensitive choice: the major part of the water molecules at the base of the droplet stay on top of the plane defined by the OH groups, so, as it would be natural, we assume the contact plane as being there, despite, in some points, being possible to observe some water molecules below this plane, due to the strong interaction and the holes in the OH distribution (especially when ρ OH A is low). We will see below how this limitation in the definition of the contact plane affects the results.
The values of θ YD are listed in Table II. A fast inspection on Table II reveals that the values of the contact angle decrease as ρ OH A increases, with the values of the extreme cases differing by a factor of 2, i.e., the contact angle that is 102.4 ○ for the non-hydroxylated surface decreases to 50.5 ○ for the fully hydroxylated surface, showing, as expected, that the presence of the OH groups enhances substantially the wettability of the silica surface. Even for the smallest level of hydroxylation considered here (ρ OH A = 1.175 nm −2 ), a significant increase in the wettability can be observed, with θ YD decreasing to 87.4 ○ . This is important from a manufacturing perspective, as it can be difficult to assure that a silica surface is fully hydroxylated, but partially hydroxylated silica surfaces also present significant better wettability. Figure 6 (red circles) presents the values for the cosine of the apparent contact angle (cos θapp). As one could expect, the observed behavior is very close to the one obtained with the thermodynamic integration calculation: as ρ OH A increases, the cos θapp values increase almost linearly. Consequently, the values for θapp also decrease almost linearly (Table II), with the value for the nonhydroxylated surface and the fully hydroxylated surface differing again by a factor of 2, i.e., 103.8 ○ for the non-hydroxylated surface and 51.6 ○ for the fully hydroxylated one. The proximity of the values presented on Table II shows the consistency between the two theoretical methods employed, but a second inspection of Fig. 6, now comparing the values for cos θapp (red circles) and cos θ YD (black squares), shows that while for most cases the values for both methods are in fair agreement, for two of the partially hydroxylated surfaces (ρ OH A = 1.175 and 2.350 nm −2 ), they are significantly different.
To understand the reason for this difference in the values, we have to look back at the approaches and assumptions involved in both methods. One hint we can get from Fig. 6 is that the values for cos θ YD seem to present an approximately linear behavior, while the values for cos θapp seem to deviate from such behavior, suggesting that some limitation when obtaining the latter could cause the difference observed in the values.
Here, we come back to the aforementioned difficulty of establishing the contact plane for a surface that presents inhomogeneity. Due to the presence of the OH groups at random sites of the silica surface, we have very inhomogenous surfaces when ρ OH A is 1.175 and 2.350 nm −2 , with regions with high density of OH groups, while in other regions, the OH groups will be completely absent (Fig. 1). In such cases, the holes in the OH distribution allow water molecules to have more space to interact closer to the OH groups, while being repealed by the slightly hydrophobic character of the siloxane bridges.
The right side of Fig. 7 presents the top views of several snapshots of the droplets on one of the partially hydroxylated surfaces (ρ OH A = 1.175 nm −2 ), displaying, for clarity, only the top layer of the silica surface and the water molecules (here, in blue color) at the adsorption layer of the droplet. Here, we can see that the water molecules tend to be positioned near the OH groups while avoiding regions with a great number of siloxane bridges.
This natural chemical behavior of the system makes determining the contact plane more complex: we assumed the plane as being on the top of the OH groups, but there is some considerable occurrence of water molecules below the plane. When ρ OH A is 1.175 and 2.350 nm −2 , this is significant, making the values obtained for cos θapp and cos θ YD considerably different. When ρ OH A is 3.525 nm −2 , while the inhomogeneity is, in some sense, the same as that when ρ OH A is 1.175 nm −2 , the size of the holes in the OH distribution are smaller, so the water molecules basically stay on top of the OH groups, making the choice for the contact plane suitable and the error neglectable, and then, good agreement was observed between the values obtained with both methods. This is an advantage of methods based on thermodynamic integration, based on the phantom wall 33 or dry surface, 35 as they can provide the value for the free energy and ultimately the contact angle, without any arbitrary choice of a contact plane, as has been discussed in a recent review. 31 It is also worth mentioning that there are some attempts to solve the above limitation, which suggests transforming the simple problem of finding a suitable plane for the contact line into a more complex problem of describing the contact plane as a 3D object in order to account for the roughness and inhomogeneity, 67  Finally, we discuss the absolute values of the contact angle themselves. Inspection of Table II shows that the obtained contact angle values for all the five cases considerably overestimate the experimental values. In particular, for a completely hydroxylated silica surface, i.e., a hydrophilic surface, it would be expected to have a contact angle of zero, while for the nonhydroxylated surface, i.e., the hydrophobic one, experimental measures place the value between 40 ○ and 80 ○ , depending on the technique. 68,69 Emami et al. 26 have already obtained theoretically values in very good agreement with those measurements using the IFF.
Our main objective was not the reproduction of the experimental values but rather to discuss the equivalence of the thermodynamic and mechanical routes and the effects of the hydroxylation on the work of adhesion and entropy. Because of this, despite the fact that the IFF allows the construction of very specialized surfaces, we opted for employing surfaces where the main difference resides on ρ OH A only, without considerable deformations. Thus, no annealing process was used, and the initial structures are very regular, even artificial to some extent (see Fig. 1). Ionization of the surfaces was not considered either: they are all neutral, which is a good representation of an acid pH case (between 2 and 4), but measurements are usually made with pH larger than 5. 68 These more idealized surfaces serve as a first approach for silica, but tend to have a decreased wettability. However, they were suitable for our purpose, as we were able to obtain clear and very well-defined relations between ρ OH A , work of adhesion, and entropy loss. Additionally, to make an assertive comparison with experimental measurements, we would have to take into account the possibility of other effects, for instance, the size dependence of the 2D droplet and the influence of the droplet shape itself, as quasi-2D droplets tend to present larger contact angle values when compared to 3D ones, because of the absence of the line tension. 38

C. Mobility of the water droplet on the silica surface
The droplet simulations also provide information about the mobility of the water droplet over the surface. As the temperature control is only performed at the middle part of the silica, the Brownian motion of the droplet on top of the surface is expected to not be affected by the thermostating process. As the droplet shape is expected to show some fluctuation during the simulation, a simple and clear approach to track the water droplet displacement on top of the surface is to follow the center of mass position of the water droplet during all the simulations. Particularly, in our case, as we are simulating a quasi-2D system, it is sufficient to know the coordinate of the center of mass along the axis of the motion, in the present case, x-axis (i.e., xcm). Computing the center of mass is non-trivial for a periodic system, and in our case, we have used the algorithm proposed by Bai and Breen. 70 The left side of Fig. 7 presents the time evolution of xcm for all the five studied cases. As it would be intuitive, the inhomogeneous presence of OH groups decreases the mobility of the droplet due to the hydrogen bond formation. In all hydroxylated surfaces, xcm stays around the same value, varying up to 1 nm on both +x and −x directions. The snapshots on the right side of Fig. 7 present the xcm position at different time steps (purple sphere) for the case of ρ OH A = 1.175 nm −2 . In addition to the initial (t = 0) and final steps (t = 10 ns), where the xcm is found at −1.818 and −1.717 nm, we also present snapshots when xcm is considerably displaced to +x (t = 1 ns) and -x (t = 7 ns) directions, being found at −1.052 and −2.413 nm, respectively. This exemplifies how the mobility of the water droplet is very limited when we have OH groups on the silica surface, as the graphic on the left side of Fig. 7 shows. For the non-hydroxylated surface, the situation is clearly different: with no hydrophilic group to prevent the droplet displacement, xcm can move as further as 4 nm The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp from its initial position within 10 ns of simulation (left side of Fig. 7, purple line). The pinning level would also be affected by the inhomogeneity of the surface. We have estimated the pinning force of flat Lennard-Jones systems in the past; 54 however, this estimation can be more complex for surfaces presenting roughness and local charges and is beyond the scope of the present work. Here, it suffices to say that the pinning of the contact line would be expected to be larger on inhomogeneous surfaces, explaining the decrease in the mobility of the water droplet presented in Fig. 7 in the case of partially hydroxylated surfaces. In addition, the pinning force could be used in a modified version of the Young-Dupré equation to obtain the contact angle values including the pinning effect on the contact line. 54 The differences in the pinning force due to different levels of inhomogeneity of the surface could also be a factor for the differences observed in the contact angle values of partially hydroxylated surfaces obtained from the Young-Dupré equation and the apparent shape of the droplet (Fig. 6). This will be investigated in the future.

IV. CONCLUSION
We have performed a study about the interaction of water with hydroxylated silica surfaces. Different OH area density values were considered, and two theoretical methods were used: thermodynamic integration implemented by the phantom wall and MD simulations of water droplets on silica surfaces. The TI method consists of a series of MD simulations of a water film interacting with the silica surface, where a repulsive potential wall is introduced to quasi-statically strip the film from the surface. Its results are free of arbitrary choices present in droplet MD simulations, such as the contact plane position and circle fitting procedures, but do not give information about Brownian motion or pinning, which are obtainable from MD simulations. Contact angle values obtained from both methods are consistent and very similar, and point out, as expected, that the hydroxylated silica surfaces have an enhanced wettability, despite considerably overestimating experimental values and results from previous simulations, due to the employed approaches used in the construction of the surface to isolate the hydroxylation effect.
Some differences of values observed for the partially hydroxylated silica surfaces seem to arise from the inhomogeneity, which makes difficult a straight definition of the contact plane, necessary to measure the angle value from the apparent shape of the droplet. Further works can reduce this difference by trying to convert the simple determination of the plane into a 3D problem, which is still a state-of-the-art problem.
We have also analyzed the interfacial entropy loss and found a linear relation with the OH area density. This can be seen as a consequence of a small entropic effect because the OH groups are flexible but short. Nevertheless, the presence of OH groups has a great impact on the interfacial entropy loss values. An OH area density of 1.175 nm −2 is already enough to have work of adhesion and entropy loss at the same magnitude.
Finally, we evaluated the mobility of the droplet on the silica surfaces and observed its decrease in inhomogeneous surfaces. This decrease can be connected with a strong pinning of the contact line and will be a subject of future investigations. ACKNOWLEDGMENTS C.B. and Y.Y. were supported by the JST CREST (Grant No. JPMJCR18I1), Japan. Y.Y., D.S., and H.K. were supported by the JSPS KAKENHI under Grant Nos. JP18K03978, JP20K14659, and JP20J20251, Japan, respectively. Numerical simulations were partly performed on the Supercomputer system "AFI-NITY" at the Advanced Fluid Information Research Center, Institute of Fluid Science, Tohoku University. We appreciate fruitful discussions with Takeshi Omori and Kotaro Oda.

APPENDIX A: EVALUATION OF CELL SIZE PARAMETERS
In Table III, we present the results for the evaluation of the cell parameters, as explained in Sec. II A. Two thermodynamic conditions were employed, as also as two different system sizes, and essentially, all the results point out the same values. They are also very close to the values that were obtained in our previous work. 29

APPENDIX B: CONVERGENCE OF THE CONTACT ANGLE VALUES
In Fig. 8, we present the time evolution of the contact angle values. Except for the case of ρ OH A = 4.700 nm −2 , all the contact angle values seem to be converged from the middle of the run. For ρ OH A = 4.700 nm −2 , additional 5 ns were performed, but the value was basically the same, so we used the value obtained within 10 ns for consistency.