On the effects of induced polarizability at the water–graphene interface via classical charge-on-spring models

Abstract

Molecular models of the water–graphene interaction are essential to describe graphene in condensed phases. Different challenges are associated with the generation of these models, in particular π–π and dispersion interactions; thus quantum and classical models have been developed and due to the numerical efficiency of the latter, they have been extensively employed. In this work, we have systematically studied, via molecular dynamics, two polarizable graphene models, denominated CCCP and CCCPD, employing the charge-on-spring model of the GROMOS forcefield, both being compatible with the polarizable water models COS/G2 and COS/D2, respectively. These models were compared with non-polarizable graphene and SPC water models. We focused the study on the water–graphene interface in two distinct systems and under the influence of an electric field: one composed of graphene immersed in water and the other composed of graphene with a water droplet above it. In the former, the orientation of water close to the graphene layer is affected by polarizable graphene in comparison to non-polarizable graphene. This effect is emphasised when an electric field is applied. In the latter, carbon polarizability reduced water contact angles, but graphene retained its hydrophobicity and the computed angles are within the experimental data. Given the significant extra computational cost, the use of polarizable models instead of the traditional fixed-charged approach for the graphene–water interaction may be justified when polarizability effects are relevant, for example, when applying relatively strong fields or in very anisotropic systems, such as the vacuum–bulk interface, as these models are more responsive to such conditions.

---

1 Introduction

Graphene is a one-atom thick material formed by a honeycomb lattice of aromatic sp2 carbons. A wide range of materials derived from it¹ have multiple potential applications, such as desalinization², drug-delivery³ and antiseptics,^4–6 among others. Most of their applications are within aqueous environments; therefore, suitable models of the interaction between water and graphene at the molecular level are essential for molecular simulation studies and the interpretation of experimental data.^7,8

The binding of molecules on surfaces can be either by chemical reactivity or by direct physical contact without the formation of any covalent or ionic bonds, i.e. chemisorption or physisorption, respectively. Due to the non-reactive nature of the water–graphene interface, the binding of water to these materials pertains to the latter process.⁹ Briefly, physisorption can emerge by dipolar interactions of permanent, induced or fluctuating dipoles, the latter are the well known dispersion or van der Waals forces.¹⁰ In this regard, the apolar nature and the delocalized π-electron cloud of graphene imply that, ideally, a graphene molecular model should include both dispersion interactions and induced polarizability.

Formally, a rigorous molecular description of the water–graphene interface requires quantum mechanical methods (QM). In this regard, the aforementioned dominance of dispersion interactions requires the use of QM models that incorporate electron–electron correlation terms; indeed, density functional studies of graphene, more often than not, include empirical factors that account for London dispersion forces.^10–14 All together, the applicability of these models is, in general, limited to systems of a couple of thousands of atoms and, in many cases, only single point calculations in the gas phase are numerically affordable.

The simpler and numerically cheaper classical force field (FF) models, in particular those employed in molecular dynamics (MD) simulations, offer the capacity to cover larger systems and take into account the system entropy.^15,16 Traditionally, classical models include a dispersion term through a mean-field approximation via the use of the Lennard-Jones (LJ) type of potentials but induced polarizability is commonly neglected, these are the so-called fixed-charge models,^17–20 even though many polarizable potentials have been developed since.^21–28 In the context of the GROMOS FF,²⁹ induced polarizability has been implemented as an iterative and self-consistent charge-on-spring (COS) model,³⁰ also known as the Drude oscillator³¹ or shell model.³² In this way, two GROMOS-compatible water models have been reported: COS/G2³³ and COS/D2²⁵. The COS/G2 force field was set to reproduce the gas phase dipole moment of water, thus “G” in COS/G2 stands for “gas”, and the 2 is the second iteration of the model. For the COS/G2 model, molecular polarizability and oxygen–oxygen LJ parameters were optimized to reproduce the heat of vaporization and density of liquid water at room temperature and pressure. This model has a slightly too large dielectric permittivity; therefore, COS/D2 addresses this issue by damping the polarizability induced by large electric fields. For the COS/D2 model, molecular polarizability and oxygen–oxygen LJ parameters were calibrated against the heat of vaporization, density, static dielectric permittivity and the position of the first peak in the oxygen–oxygen radial distribution function at room temperature and pressure. The “D” stands for damped and the 2 is the second iteration of the model.

The majority of studies that simulate graphene, models it as a honeycomb lattice of sp² carbons with a fixed charge of zero.^34–37 This inability to model electronic polarisation may lead to an inaccurate description of the graphene–water interaction.³⁸ For this reason, the effects of including explicit electronic polarization in the context of classical FF need to be assessed. Recent MD studies have compared the interaction between different water models, including a polarizable water model (SWM4_DP³⁹), with non-polarisable and polarisable graphene.^35,40 Multiple analyses were done, whose results suggested that the water structural properties were not dependent on the graphene polarizability. Interestingly, the use of charged graphene influenced water dynamics,^2,35,40 suggesting that the use of an external electric field can affect the interactions between water and graphene.

A fundamental challenge that arises in the parametrization of the water–graphene interaction is the existence of accurate experimental data to be calibrated with; a common approach is to adjust the water–carbon LJ potential to reproduce water contact angle (WCA) measurements.^41,42 The WCA is defined as the angle formed between a water droplet and the surface; WCA is derived from Young's equation and is the result of the forces between the water surface tension and the interfacial tension between water and graphene.⁴³ WCA values define the chemical nature of the surface, with higher and lower values describing the hydrophobicity or hydrophilicity of the surface, respectively. Regarding graphene (and derived materials) the latter strategy has two main issues; (i) a precise WCA value of water over pristine graphene is under debate due to different experimental problems. These include how graphene is synthesized,⁴⁴ the structure in which the graphene layer is supported on,⁹ the WCA measurement methodology⁴⁵ and the graphene purity.⁴⁶ All these factors influence the physicochemical nature of the surface, thus it is not strange that the reported water–graphene or water–graphite contact angles range from 33 to 143;⁴⁴ (ii) molecular simulations, in general, are on a microscopic or nanoscopic scale in comparison with experimentally determined macroscopic WCA (millimetres). For nanodroplets the contact angle depends on the line tension, which is a constant term added to the Young's equation to correct the nanoscopic WCA⁴⁷; still, the existence, the sign, and the value of line tension are subject to debate.^48,49 In this regard a method based on the estimation of the work of adhesion via rigorous free-energy calculations circumvents this finite-size effect⁵⁰; recently Sresht and collaborators⁵¹ have shown that this method and the one that directly simulates nanoscopic droplets, earlier proposed by others,^41,42 converge to similar values when the first two water layers of the droplet are removed. Employing the work of adhesion method, polarizable graphitic surfaces parameterized against QM calculations highlighted that induced polarizability influences WCA values.⁵²

In this work we have explored the effects of induced polarizability, in light of the influence of external electric fields, for the water–graphene interaction via MD simulations in the context of the GROMOS FF. Regarding graphene, we tested two polarizable COS models, denominated CCCP and CCCPD, being compatible with the COS/G2 and COS/D2 water models, respectively. CCCP or CCCPD account for the polarizability of each carbon, with CCCPD damping over-polarization. These models where compared with the traditional zero charges models for graphene immersed in water and combinations of these. In detail, structural and dynamic descriptors such as atomic densities, dipolar orientations and diffusivities were compared; furthermore, nanoscopic WCA for water droplets of different radii was computed, allowing for the estimation of macroscopic WCA which was then compared to the WCA reported in the literature. Our results indicate that polarizable systems do not render relevant differences with respect to fixed charged models, under zero field and bulk conditions, with most of the properties of the graphene–water interaction mainly dependent on the water model; the use of an applied electric field did generate rather small changes in the water structure and dynamics. WCA did suffer substantial changes with the use of COS/G2 and COS/D2, producing much higher WCA values when compared to the traditional SPC model. On the contrary, polarizable graphene (CCCP or CCCPD) reduced WCA values (for all water models) with respect to non-polarizable graphene (CCC). Macroscopic WCA values for both polarisable systems are in line with the experimental WCA values. In light of these results, the use of more sophisticated and numerically costly COS models in the graphene–water interaction may be justified for cases in which increased electrostatic responsiveness is important, like WCA calculations.

2 Methods

2.1 Modelling and simulation

2.1.1 Polarizable COS models. The self-consistent damped charge-on-spring (COS) formulation²² models induce dipoles as separations of point charges. An induced dipole [small mu, Greek, vector]^ind_i on each atom is determined bywhere α_i is the polarizability of the virtual site on atom i and E⃑_i is the electric field (from the permanent atomic charges and [small mu, Greek, vector]^ind_i) taken at the location of the COS charge; p_i is a parameter that defines the damping level and E_0,i is a certain field strength which serves as a truncation parameter. When E_0,i = E_i the COS models reduce to the undamped formulation. The self-polarization contribution to the potential energy is thenInitially, we used the minimized coordinates of non-polarizable systems (CCC-SPC) as the starting point. We created two polarizable graphene, CCCP and CCCPD, models compatible with the polarizable water models COS/G2³³ and COS/D2,²⁵ respectively. CCCP has a virtual COS charge of −8 e and an α of 1.1 × 10⁻³ nm³ at each carbon atom.⁴⁰ CCCPD has the same virtual COS charge and polarizability. For COS/D2 and CCCPD models E₀ equals to 140 (kJ mol⁻¹ nm⁻³)^1/2 and p was equal to 8. Both graphene models employ the SP2 carbon LJ parameters of the 53A6 GROMOS FF,²⁰ for more details on the specific parameters for all models see Tables S1 and S2 (ESI†). To explore the effects of graphene polarizability, α values ranging from 0 to 1.5 × 10⁻³ nm³, with a granularity of 0.15 × 10⁻³ nm³, were tested. A combination of polarizable and fixed charge models was also tested.

2.1.2 Simulated systems. A periodic (in the x and y plane) graphene layer (CCC) was generated with a bond length of 0.14 nm using the Carbon Nanostructure Builder plugin of VMDv1.93.,⁵³ with an edge length of 3.1 nm × 3.3 nm × 4.0 nm. The system was solvated, containing a total of 1213 water molecules (SPC water model⁵⁴, COS/G2³³ or COS/D2²⁵), resulting in a density of 999 kg m⁻³ (see Fig. 1).

Fig. 1 Representation of a graphene layer (van der Waals representation in grey) immersed in water (surface representation in light blue).

For WCA computations, we constructed a periodic graphene layer in a box of 14.2 nm × 14.2 nm × 18 nm. Subsequently, a water droplet resembling a disk with a radius of 3 nm was placed over the graphene layer. The water droplet was composed of 2594 SPC water molecules, for more details please refer to Fig. 2. For line tension estimations, 5 droplets with radii ranging approximately from 3 nm to 4 nm were also constructed.⁴¹

Fig. 2 System of a water droplet over graphene after 2 ns of simulation. Water molecules and graphene are represented using a van der Waals representation.

2.2 Molecular dynamics simulations

MD simulations were carried out with the GROMOSv1.3.0 software package (http://www.gromos.net/).⁵⁵ The simulations were run at constant volume and temperature. The geometry of molecules was constrained by using the LINCS algorithm⁵⁶ with an order degree of 4. The temperature was weakly coupled⁵⁷ to a bath of 298.15 K with a relaxation time of 0.1 ps. The van der Waals and electrostatic interactions were calculated using a triple-range cutoff pairlist with radii of 0.8 and 1.4 nm. The short-range interactions were calculated every time step by updating the molecular pair list for distances smaller than the first cut-off radius of 0.8 nm. For the intermediate range of distances between 0.8 and 1.4 nm, the pair list was updated every 10 fs. The long-range electrostatic interactions, beyond the outer cut-off of 1.4 nm, were represented by a reaction field⁵⁸ with a relative dielectric permittivity (ε_RF) of 78.5 and 61 for polarizable and SPC water, respectively. The equations of motion were integrated using the leapfrog algorithm with a time step of 2 fs. The velocities of the atoms at the beginning of the simulation were assigned from a Maxwell distribution at 298 K. During the runs, configurations of the system were saved every 0.5 ps. After 500 ps, production runs lasted on average 1 ns each.

To explore the responsiveness of the COS models, external electrical fields of 0.001, 0.005, 0.01 and 0.05 e nm⁻² were applied. These fields produce a potential difference across the simulation box of 5.8, 28.8, 58 and 288 mV, respectively. These potentials are well within biologically occurring potentials.⁵⁹

For the water droplet over a graphene layer system, we ran an equilibration simulation of 500 ps and a production simulation of 2 ns using the previously described simulation parameters, except for ε_RF which was set to 1 in order to simulate vacuum conditions. To avoid water evaporation and overheating, translations and rotations of the water molecules were thermostated independently.⁶⁰

2.3 Analysis

A plethora of analyses were done using the GROMOS++⁶¹ package, covering different structural and dynamical properties of water.

2.3.1 Number density function. We calculated the oxygen and hydrogen density along the normal of the graphene plane. The oxygen density profile permitted us to determine the first hydration layer used in the following analyses. The first hydration layer corresponds to water molecules with a distance below 0.509 nm from the center of the carbon atoms. This value corresponds to the basin between the two peaks of the density profile created by the first two water layers close to graphene (see Fig. 3).

Fig. 3 Water oxygen (ρ_O) and hydrogen (ρ_H) density and dipolar orientation (P1) profiles along the graphene normal, as function of polarizability (α). In A and C, ρ_O for the CCCP-COS/G2 and CCCPD-COS/D2 systems under zero field conditions, respectively; B and D present ρ_O under the effects of an external electric field. In E and G, ρ_H for the CCCP-COS/G2 and CCCPD-COS/D2 systems under zero field conditions, respectively; F and H present ρ_H under the effects of an external electric field. In I and K, P1 for the CCCP-COS/G2 and CCCPD-COS/D2 systems under zero field conditions, respectively; J and L present P1 under the effects of an external electric field. Left and right sides, represent the top and bottom layers, respectively.

2.3.2 Dipolar orientation. We characterized the orientation distribution of the dipole vector long the z axis by calculating the first Legendre polynomial

P₁ = 〈cos(θ_z)〉, (3)

where θ_z is the angle between the graphene normal, along the z axis, and the dipole moment.

2.3.3 Mean residence times. To quantify the dynamics of the first hydration layer of graphene, we calculated the survival probability of water molecules,⁶²where P_i is equal to 1 if a particle i is located inside the hydration layer at time t_k. This is calculated for N loaded molecules, in which the angular brackets denotes the average probability that a particle remains inside the hydration layer at time t₀ + t, with multiple time origins t₀. This function decays exponentially with time and can be fitted using:where τ_s is the mean residence time.

2.3.4 Diffusion coefficient. We calculated the mean-square displacement (Δ(t)_xy) of water molecules moving in the x − y plane within the first hydration layer of graphene, according towhere r_i corresponds to the x and y components of the position vector of each oxygen atom, of each water molecule i, and N_m is the total number of molecules. The average is over multiple origins t₀ and molecules in the computational box which were within the first hydration layer. Finally, the self-diffusion coefficient D_xy was obtained from the long-time limit of Δ(t)_xy according to the Einstein relation,

2.3.5 Rotational relaxation. To get insights on the rotational relaxation of water molecules, the re-orientational correlation function, C^μ₁, was calculated:

C^μ₁(t) = 〈μ_i(t₀)·μ_i(t₀ + t)〉_τ (8)

which corresponds to the autocorrelation function of the ith water dipole vector μ. The average is over multiple origins t₀ and for all molecules that were within the first hydration layer. In general C^μ₁(t) exhibits an exponential decay over time which is fitted using:where τ_μ is the rotational relaxation time for the dipole vector.

2.3.6 Water contact angle. The water contact angle (WCA) of a water droplet over a graphene layer was estimated in three different steps, as done by others.^41,42 The first step is to calculate the surface boundary of the droplet over the surface. For this reason, it is necessary to compute the radial density of different horizontal layers, whose reference point is the centre of geometry of the droplet in the x–y plane. Each layer was moved upwards by 0.05 nm in the z-direction, starting from the bottom and finishing at the top of the droplet; the average density of water molecules was radially binned every 0.3 nm afterwards and to obtain the point where the density drops to zero the radial density values of each layer were fitted using a curve which represents the decrease of the radial density:where ρ(r) is the radial density, a and b are proportional to the bulk liquid density and the thickness of the density decrease, respectively, and r₀ is the surface boundary. The second step is to obtain the circumference of the droplet using the boundary surface of the droplet. Fittings were done using a loss function to reduce the influence of outliers on the solution of the non-linear least squares problem and to avoid the noise at the top of the droplet, see Fig. S4 (ESI†) for representative plots of this procedure. Finally, the WCA was obtained by computing the angle of the tangential line at the intersection point of the circumference with the surface. This was done for each simulation frame.⁴¹

Lately, Sresht et al.⁵¹ suggested that better agreement with the work of adhesion method⁵⁰ is achieved by removing the first two water layers; in our case, it is expected that these two layers are the most affected by graphene polarizability, thus to avoid the noise produced by density oscillations at the boundaries, instead of not including these two layers in the density calculations, we employed the aforementioned loss function in the non-linear fitting procedures (see Fig. S4, ESI†).

2.3.7 Macroscopic contact angle θ_∞. The macroscopic contact angle θ_∞ is related to the microscopic contact angle θ through the modified Young's equation⁴⁹where γ represent the surface tension between different phases (solid, S; liquid, L; and vapor phase, V), τ is the line tension, θ the contact angle and r_B is the droplet base radius. For larger droplets the third term tends to 0. Accordingly, the macroscopic contact angle θ_∞ can be defined as:

cos θ_∞ = (γ_SV − γ_SL)/γ_LV. (12)

If we subtract by γ_SL and divide by γ_LV on both sides of eqn (11), and then replace the left side using eqn (12) and rearrange the terms, we can obtain a relationship between the microscopic water angle θ with the macroscopic water angle θ_∞:

This permits us to estimate θ_∞ by plotting contact angles as a function of r_B⁻¹ for water droplets of different sizes. cos θ_∞ is the y-intercept of a linear fit of the previous plot with the slope being .

3 Results and discussion

In this section the structural and dynamical effects of polarizable models are systematically evaluated. Firstly, we determined the responsiveness, upon the application of an external electric field, of different graphene polarizability levels (α); following this, combinations of polarizable and fixed-charge models (for both water and graphene) are explored; finally WCA are computed for the same combinations of polarizable and non-polarizable models.

3.1 Polarizability and its response to external electric fields

For a graphene layer immersed in water (see Fig. 1), we tested the influence of induced polarizability on the localization and structuring of the COS/G2 and COS/D2 water models in the presence of an external electric field. Thus, we created CCCP-COS/G2 and CCCPD-COS/D2 systems with different polarizability levels of the graphene, that is α values ranging from 0 to 1.5 × 10⁻³ nm³. An applied electric field strength of 0.05 e nm⁻² was exerted along the positive z axis, that is along the normal of the graphene top layer. This field magnitude, which renders a potential of around 250 mV, is close to the limit where no statistically significant water dipolar alignment was observed for a fixed charged model.⁶³ Consequently, any orientational preference is mainly due to water and graphene polarizabilities.

In Fig. 3 the oxygen density (ρ₀), hydrogen density (ρ_H) and dipolar alignment (P1, see eqn (3)) along the graphene normal for the aforementioned systems are presented. Under zero field conditions (see first and third columns of Fig. 3) there are no significant effects due to the choice of α, only a slight increment in ρ_H at close proximity, a consequence of water dipolar alignment. In this way, P1 values exhibit a moderate alignment with the graphene normal, which is augmented as α increases; no substantial differences between the CCCP-COS/G2 and CCCPD-COS/D2 systems exist, only a slight increase of P1 values for the CCCP-COS/G2 system at distances around 0.2 nm for the upper range of α values.

On the other hand, the external electric field (see second and fourth columns of Fig. 3) breaks the symmetry between the top and bottom sides, with changes in ρ_H due to an enhanced interaction and reduced P1 at the top and bottom layers, respectively. Qualitatively, dipolar alignment (at the top-layer) within the graphene–water interface builds up, with an enhanced structuring for higher α values; these values tend to converge for α above 1.0 × 10⁻³ nm³; for the bottom layer, dipolar alignment with the graphene normal is substantially reduced. Polarizable atoms have an associated virtual charge to them, that will move according to the electric field of the system. In this case, the graphene negative virtual charge will move against the field direction; accordingly, an electric field applied towards the positive z axis generates a surface in which the bottom side will concentrate a negative charge density and the upper side a positive one, thus explaining the above results. The CCCP-COS/G2 and CCCPD-COS/D2 profiles are very similar, only at the first P1 peak the CCCPD-COS/D2 system presents a minimal reduction due to polarisability damping.

Overall, under zero field conditions there are no substantial differences between fixed charged models (α = 0 × 10⁻³ nm³) and the polarizable ones; upon the application of a relatively weak field and as expected, there is an enhanced electrical responsiveness, which reaches a plateau; for these field strengths, both polarizable models present almost equal results. In light of the aforementioned, the undampened models, with α = 1.1 × 10⁻³ nm³, for graphene were chosen for further analyses.

In the next section, we compare polarizable, fixed charged and combinations of these models to quantify the influence of graphene and water polarizability, respectively.

3.2 Polarizable vs. fixed charged models

3.2.1 Dipolar structuring. As shown in the previous section, P1 values at the graphene–water interface exhibited a decrease for α = 0 × 10⁻³ nm³ at close distances from the graphene layer when compared to the upper limits of induced polarizability, α > 1 × 10⁻³ nm³ (see green traces vs purple traces in Fig. 3, panels I and K), which may suggest that this structuring is due to induced polarizability; nonetheless, the COS/G2 and COS/D2 models have different VdW interactions when compared to SPC water (see Table S2, ESI†).²⁵ In this manner, in Fig. 4A we compared P1 values for a combination of polarizable and fixed charged models; for completeness, the ρ_O and ρ_H profiles are presented in Fig. S1 (ESI†). The P1 profiles of the CCCP-SPC and CCCP-COS/G2 pairs (yellow and red traces, respectively) exhibit the highest values at the first hydration shell; the latter indicates that the graphene polarizability significantly influences water dipolar structuring; quite interestingly, only using the SPC model (yellow and blue traces) shows higher P1 values with respect to the COSG2 model (green and red traces).

Fig. 4 Water dipolar alignment along the graphene normal for simulated systems; (A) P1 values for polarizable, fixed charged models and combinations of these. (B) P1 values for the hybrid CCC-COS/G2 system at different electric field strengths. (C) P1 values for the fully polarizable CCCP-COS/G2 system at different electric field strengths. Left and right sides, represent the top and bottom layers, respectively.

Additionally, we explored whether these descriptors were sensitive to different thermostats by running simulations (under zero field conditions) employing the Nosé–Hoover algorithm as shown in Fig. S5 and S6 (ESI†). Overall there are no relevant effects when switching thermostats; density profiles are almost identical and P1 values only exhibit a small increase (0.005) for the CCCP-COS/G2 pair, which does not alter the conclusion that graphene polarizability increases dipolar alignment within the first hydration shell.

3.2.2 Interfacial water kinetics. To explore the effects of induced polarizability on the water dynamics at the graphene–water interface, we computed mean residence times (τ_s), diffusivity along the graphene plane (D_xy) and dipolar relaxation times (τ_μ) of interfacial waters – placed within the second and first density peaks of ρ_O shown in Fig. 3 – for a combination of polarizable and fixed charged models. Results are plotted in Fig. S2 (ESI†) and summarized in Table 1.

Table 1 Mean residence times (τ_s), plane diffusivity (D_xy) and dipolar relaxation times (τ_μ) of interfacial water for combinations of polarizable and fixed charge models

	SPC	COS/G2
τ s (ps)
CCC	15.5 ± 1.0	19.6 ± 1.6
CCCP	15.3 ± 0.7	19.8 ± 1.3
D xy (nm^2 ns^−1)
CCC	4.6 ± 0.6	3.2 ± 0.5
CCCP	4.7 ± 0.6	3.0 ± 0.4
Bulk	4.5 ± 0.3	2.3 ± 0.2
τ μ (ps)
CCC	3.5 ± 0.2	7.3 ± 0.6
CCCP	3.4 ± 0.2	7.3 ± 0.5
Bulk	3.5 ± 0.2	7.4 ± 0.7

---

The τ_s, D_xy and τ_μ values of interfacial waters are mostly dependent on the water model. The only noticeable effect of graphene when compared to bulk conditions involves D_xy; in detail, results show that the diffusion coefficients of water molecules close to graphene are slightly higher, in particular for the COS/G2 model, when compared to bulk conditions²¹ due to the smoothness and almost frictionless graphene surface, an effect that has been thoroughly described for water fluxes within narrow carbon nanotubes.^64,65 This increase is more significant for the polarizable model as the potential energy surface of the bulk is rougher when compared to that of traditional fixed charged models, a fact revealed by the longer τ_μ values for the COS/G2 model.

We explored the effects, on the same dynamic descriptors (τ_s, D_xy and τ_μ), of a constant electric field applied along the graphene normal. Results are plotted in Fig. S2 (ESI†) and summarized in Table 2. In this case, as the dipolar alignment along the field breaks the symmetry between the top and bottom layers, τ_s, D_xy and τ_μ were independently estimated for both layers. Noticeable effects arise only for the CCCP-SCP pair; for example, when comparing the top layer against the bottom one, a consistent increase in τ_s and decrease in D_xy for the highest field strength is observed; furthermore and for both layers, D_xy tends to increase with the field intensity, as well; for the CCC-SCP pair, this tendency is also present but it is not statistically significant. On the other hand, the COS/G2 models fully nullify this effect. Overall, all these descriptors, under zero field conditions and in the presence of weak external electric fields, are mainly dependent on the water models. Polarizable graphene leads to a slight increase in τ_μ for higher fields when employing a polarizable water model as well, but always within the values observed in the bulk. The use of combined fixed charge and polarizable models (at higher fields) significantly increases the asymmetry in τ_s and to some degree D_xy; these effects are cancelled by the use of polarizable water.

Table 2 Mean residence times (τ_s), plane diffusivity (D_xy) and dipolar relaxation times (τ_μ) of interfacial water for combinations of polarizable and fixed charge models at different electric field strengths applied along the graphene normal

	E field [e nm^−2]	SPC		COS/G2
		Top	Bottom	Top	Bottom
τ s (ps)
CCC	0.001	15.5 ± 0.9	15.5 ± 1.0	19.7 ± 0.7	20.2 ± 1.2
	0.005	15.6 ± 0.8	15.1 ± 0.6	19.5 ± 1.0	20.1 ± 1.0
	0.01	15.6 ± 1.1	15.4 ± 0.5	20.3 ± 1.8	19.8 ± 1.3
	0.05	16.3 ± 0.8	14.6 ± 0.9	20.0 ± 1.8	19.3 ± 1.0
CCCP	0.001	15.1 ± 0.6	16.0 ± 1.0	20.4 ± 1.2	20.6 ± 1.3
	0.005	15.7 ± 0.9	15.1 ± 0.6	20.3 ± 1.1	19.7 ± 1.4
	0.01	16.2 ± 0.8	14.9 ± 0.5	20.0 ± 1.4	20.0 ± 1.3
	0.05	19.0 ± 1.0	13.8 ± 0.7	20.4 ± 1.3	19.5 ± 1.1
D xy (nm^2 ns^−1)
CCC	0.001	4.8 ± 0.8	4.7 ± 0.6	2.9 ± 0.3	3.0 ± 0.3
	0.005	4.7 ± 0.8	4.7 ± 0.6	3.2 ± 0.4	3.2 ± 0.3
	0.01	4.6 ± 0.6	5.1 ± 0.8	3.1 ± 0.5	3.0 ± 0.4
	0.05	4.5 ± 0.6	5.0 ± 0.7	3.2 ± 0.6	3.2 ± 0.5
CCCP	0.001	4.6 ± 0.5	4.8 ± 0.6	3.2 ± 0.3	3.0 ± 0.3
	0.005	4.5 ± 0.4	4.6 ± 0.7	3.0 ± 0.4	3.2 ± 0.5
	0.01	5.0 ± 0.5	4.8 ± 0.5	2.9 ± 0.4	3.1 ± 0.4
	0.05	5.7 ± 1.0	6.3 ± 0.8	3.3 ± 0.6	3.5 ± 0.4
τ μ (ps)
CCC	0.001	3.6 ± 0.2	3.6 ± 0.2	7.2 ± 0.5	7.7 ± 0.8
	0.005	3.5 ± 0.2	3.5 ± 0.2	7.6 ± 0.5	7.4 ± 0.6
	0.01	3.5 ± 0.2	3.5 ± 0.3	7.5 ± 0.6	7.3 ± 0.6
	0.05	3.5 ± 0.2	3.5 ± 0.2	7.2 ± 0.8	7.3 ± 0.4
CCCP	0.001	3.5 ± 0.2	3.4 ± 0.3	7.3 ± 0.6	7.1 ± 0.5
	0.005	3.5 ± 0.2	3.5 ± 0.2	7.5 ± 0.7	7.2 ± 0.5
	0.01	3.3 ± 0.2	3.5 ± 0.2	7.6 ± 0.6	7.3 ± 0.5
	0.05	3.6 ± 0.2	3.8 ± 0.2	7.8 ± 0.5	7.9 ± 0.7

---

3.3 Water contact angle

In the previous analyses, structural and dynamic descriptors were mainly functions of the employed water model, with negligible effects when induced polarizability was included for the graphene model. We extended our analyses in a situation in which the responsiveness of polarizable models may produce relevant effects and also a comparison with experimental data could be carried out. In this way, we computed the water contact angle (WCA) from MD simulations of a water droplet placed on top of a graphene layer, as shown in Fig. 2, for fixed charge, polarizable models and for combinations thereof. We also tested the dampened CCCPD model, as strong polarization at the water–graphene–vacuum interface may produce overpolarization. The results are presented in Fig. S3 (ESI†) (WCA time series) and summarized in Table 3.

Table 3 Microscopic water contact angles

	SPC	COS/G2	COS/D2
CCC	77.8° ± 0.29	98.1° ± 0.68	91.2° ± 1.28
CCCP	68.3° ± 0.12	92.8° ± 0.92	—
CCCPD	67.5° ± 0.42	—	85.5° ± 0.54

---

Overall, the WCA is more affected by the water model than the graphene model. In detail, polarizable water models have a higher contact angle (∼90°) than non-polarizable water model (∼70°), reflecting that induced polarizability (and the reduced LJ parameters of the COS/G2 and COS/D2 models) reduces the water–graphene interaction; COS/D2 has a slightly lower WCA value than COS/G2 which is explained by the higher ε value for the LJ potential of the dampened model (see Table S2, ESI†). When exploring the influence of induced graphene polarizability, WCA values are reduced; briefly, there is a decrease of ∼10 for SPC water molecules going from a CCC to CCCP or CCCPD; regarding the COS/G2 and COS/D2 water models, this decrease is approximately ∼6. These results indicate that polarizable graphene interacts more favourably with water than non-polarizable graphene.

Others have reported that WCA values are mainly dependent on the surface tension (γ) of the water model.⁵⁰ In this regard, γ values have been reported for SPC, COS/G2 and COS/D2.²⁵ The SPC model has a lower surface tension (γ = 48.5 mN m⁻¹) when compared to COS/G2 (γ = 59.0 mN m⁻¹) and COS/D2 (γ = 63.6 mN m⁻¹) models and also presents a lower WCA. On the other hand, both polarisable models essentially present a very similar γ, this is expected as their main difference is the polarisability dampening to avoid over-polarization, which under bulk and at zero field conditions is negligible (see Fig. 3J and K); we can conclude that polarisable water models, which exhibit higher γ values, enhance water–water interactions and therefore increase the WCA with graphene.

We additionally explored whether relatively strong electric fields along the graphene normal of +0.1 e nm⁻² and −0.1 e nm⁻² (which produce potentials of around 1.5 and −1.5 V, respectively) for the CCCP-COS/G2 pair influence the WCA values. The fields rendered WCA values of 93.1° ± 1.21 and 93.3° ± 0.89 for the positive and negative fields, respectively; thus, no significant effects on WCA values were observed for these field strengths, as reported by others.^66,67

To compare our results with experimental observations, the macroscopic contact angle (θ_∞) was estimated by computing (microscopic)-WCA for water droplets of different sizes (for the CCCP-COS/G2 and CCCPD-COS/D2 pairs); for more details please refer to the methods section and eqn (11)–(13). Accordingly and in Fig. 5, the cosine of WCA as a function of the inverse droplet radius is shown, with the slope and intercept of these curves employed to estimate the line tension and θ_∞, respectively which are shown in Table 4. The line tension for the systems CCCP-COS/G2 and CCCPD-COS/D2 are −5.3 × 10⁻¹¹ J m⁻¹ and −4.3 × 10⁻¹¹ J m⁻¹, respectively, and θ_∞ estimations are 109.5° and 97.3°, respectively. These line tension values are close to the experimental consensus, reported to be in the order of 10⁻¹¹ J m^{−1 48, 49} with the negative sign in line with reported values for hydrophobic surfaces;^48,78 the θ_∞ are on the upper range of experimentally reported WCA,⁴⁴ thus reflecting the hydrophobic nature of the CCCP and CCCPD graphene models, when combined with polarizable water models. Given the aforementioned spread of experimentally derived WCA, it becomes interesting to qualitatively compare our results with recent similar computational studies. Therefore in Table 5 a set of reported WCA values of similar systems (in sizes and physicochemical nature) are presented. All in all, the addition of induced polarizability renders the WCA much more in line with the majority of these reports, when compared to the non-polarizable SPC-CCC pair, which produces a somewhat hydrophilic graphene surface. A direct quantitative comparison can only be carried out for θ_∞; in this regard, the θ_∞ estimations of 109.5° and 97.3° for the CCCPD and CCCP models, are within the range of the reported MD-derived θ_∞ values.

Fig. 5 Cosine of the WCA as a function of the multiplicative inverse of the droplet radius r_B for the systems CCCP-COS/G2 (red dots) and CCCPD-COS/D2 (green dots). Simple linear regressions were done for these systems (CCCP-COS/G2, red line; CCCPD-COS/D2, green line); the slope and intercept of which were used to calculate the tension line and the macroscopic water contact angle, respectively.

Table 4 Line tension (τ) and θ_∞

	CCCP-COS/G2	CCCPD-COS/D2
τ	−5.3 × 10^−11 J m^−1	−4.3 × 10^−11 J m^−1
θ ∞	109.5°	97.3°

---

Table 5 Reported WCA (θ) and macroscopic WCA (θ_∞) for similar systems

Computational system and methodology	θ, θ∞/θexp	Ref.
Water over graphene; Born–Oppenheimer QMD	87°, —	68
Flexible SPC over graphene; MD	—, 86°	69
CVD single-layer graphene on Cu substrate; exp.	86.2°	69
SPC/E over graphene monolayer; MD	104°, —	70
CVD graphene on Si-based substrate; exp.	<95°	70
SPC/E over graphene; MD	91°, —	71
SPC/E over graphene; MD	100.7°, —	72
SPC/E over double-layer graphene; MD	94.3°, —	73
SPC/E over graphene-coated copper; MD	92.2°, —	74
SPC/E over double-layer graphene; MD	97.1°, —	75
TIP4P/2005 over graphene; MD	31°–147°, 24°–150°	48
SPC/E over single-layer graphene; MD	94.9°, —	76
Captive bubble on graphene; exp.	42°± 3	77

---

The COS/G2 and COS/D2 models are derived from the primary SPC water model; hence, this systematic comparison between these three models gives the closest view of what the effect of induced polarization is. Nonetheless, there are iterations of the SPC model, with the SPC/E version⁷⁹ being widely employed in the literature; thus it becomes interesting to compare our WCA estimations with recently reported values employing SPCE/E. Table 5 shows that microscopic WCA estimations range from 91° to 104° for SPC/E over non-polarizable graphene or similar materials; these are in line with our WCA estimations for the COS models (see Table 3). The SPC/E was re-parametrized to better reproduce the water dielectric constant without any difference in the LJ parameters, consequently it has a bigger dipole moment than SPC (2.27 and 2.35 Debye, respectively); the latter renders a γ of 63.6 mN m⁻¹ for SPC/E which is quite similar to the ones calculated for both COS models and in good agreement with the experimental data.^25,80 This is in line with the larger reported WCA values for SPC/E when compared to SPC over a fixed charge graphene model. Quite interestingly, later reports suggest that graphene is partially hydrophylic^45,77,81,82 with a WCA of around 40°. In this regard, Table 3 reveals that the primary SPC model renders a WCA closer to this later experimental evidence (mainly due to its lower γ), particularly when including polarizability in the graphene layer because of enhanced dipolar interactions (see Fig. 4A).

To investigate how popular fixed charge water models are affected by graphene polarizability, we carried extra MD simulations of bulk SPC/E and TIP3P⁸³ with a CCC or a CCCP graphene layer and computed ρ_O, ρ_H and P1 values along the graphene normal, which are shown in Fig. S7 (ESI†). There are no significant differences in oxygen and water densities when including induced polarizability and using SPC/E or TIP3P; as expected, the CCCP model further increases dipolar alignment of water with graphene when compared with the CCC-SCP/E or CCC-TIP3P pairs (see Fig. S7C and F, ESI†) and this alignment is enhanced with respect to SPC (see Fig. 4A) which can be explained by the larger dipole moments of SPC/E and TIP3P than SPC. This in combination with the polarizable graphene model deepens the dipole–dipole interactions between water and graphene; therefore, it is expected that the SPC/E-CCCP pair will render WCA values below the ones reported in the literature (see Table 5), akin the tendency presented in Table 3 for all the studied water models herein. This whole analysis suggests that an efficient and rather simple approach to calibrate the water–graphene interface in light of the recently reported hydrophilicity of graphene and in the context of well established fixed charge bio-molecular force fields is to include induced polarizability in the graphene model while retaining the fixed charge formulation for water.

4 Conclusions

In this work we have systematically explored the effect of induced polarizability, via COS models, of the graphene–water interface. Our results indicate that the water–graphene interaction of a solvated graphene layer using polarizable systems does not render relevant differences when compared to non-polarizable models, under zero field conditions. Most of the properties of the graphene–water interaction mainly depend on the water model, except for those related with dipolar orientation at the interface with graphene. Only under the influence of an applied electric field we have observed changes in the water structure and dynamics.

On the other hand, when estimating WCA values, there are significant effects, which are not only due to the employed water model. Briefly, the COS/G2 and COS/D2 models render higher WCA values when compared to the SPC model. Moreover, the use of polarizable graphene (CCCP or CCCPD) reduced WCA values (for all water models) with respect to non-polarizable graphene (CCC). Both polarizable models allowed us to estimate macroscopic WCA values, which are within experimentally determined values and in line with other similar MD studies. Indeed the importance of polarizability on WCA estimations has been previously pointed out by Misra and Blankschtein,⁵² whereby employing the work of the adhesion method to compute WCA over polarizable multilayered graphene paramertrized against QM calculations has shown that induced polarizability greatly affects interfacial entropy, thus evidencing the importance of dipolar induction at the graphene–water interface.

Regarding the underlying mechanism of the graphene–water interaction, our results clearly show that graphene polarizability increases water dipolar alignment along the graphene normal (see Fig. 4); this implies that the water–graphene interaction is enhanced when including polarizability into graphene, not regarding the water model. The latter predicts that WCA values when using the CCCP and CCCPD models should be lower when compared to the CCC models, which is exactly what was obtained (see Table 3). Therefore, the main effects of graphene's induced dipoles are an enhancement of the water graphene dipole–dipole interactions at the first water layer, not regarding the nature of the water model, which macroscopically renders lower WCA values.

Overall, the more simplistic and widely employed fixed charge models for the graphene–water interaction, when carefully calibrated, produce similar results under bulk and zero field conditions with respect to the more sophisticated and numerically costly COS models only differing in the water dipolar alignment at the first water layer when including polarizabiity in graphene. Thus, only for conditions in which an increased electrostatic responsiveness is important, such as WCA calculations, the use of COS models is justified and may provide some further insight into the molecular mechanisms at the graphene–water interface.

Note added after first publication

This article replaces the version published on 16 March 2022, which contained errors in Table 1.

Conflicts of interest

There are no conflict of interest to declare.

Acknowledgements

FONDECYT grant 1180987. CONICYT-PFCHA/DOCTORADO BECAS NACIONAL/2020-21201020 (to M. B. U.). The Centro Interdisciplinario de Neurociencia de Valparaiso (CINV) and The Millennium Nucleus in NanoBioPhysics (NNBP) are funded by ANID-Millennium Science Initiative Program-, projects IACE210014 and NCN2021_021. Access to the supercomputing infrastructure of the National Laboratory for High-Performance Computing was provided through Grant ECM-03 (Powered@NLHPC).

References

1. A. K. Geim and K. S. Novoselov, Nat. Mater., 2007, 6, 183–191.
2. A. Striolo, A. Michaelides and L. Joly, Annu. Rev. Chem. Biomol. Eng., 2016, 7, 533–556.
3. H. Yue, W. Wei, Z. Gu, D. Ni, N. Luo, Z. Yang, L. Zhao, J. A. Garate, R. Zhou, Z. Su and G. Ma, Nanoscale, 2015, 7, 19949–19957.
4. O. Akhavan and E. Ghaderi, ACS Nano, 2010, 4, 5731–5736.
5. W. Hu, C. Peng, W. Luo, M. Lv, X. Li, D. Li, Q. Huang and C. Fan, ACS Nano, 2010, 4, 4317–4323.
6. S. Liu, T. H. Zeng, M. Hofmann, E. Burcombe, J. Wei, R. Jiang, J. Kong and Y. Chen, ACS Nano, 2011, 5, 6971–6980.
7. V. C. Sanchez, A. Jachak, R. H. Hurt and A. B. Kane, Chem. Res. Toxicol., 2012, 25(1), 15–34.
8. Y. Zhang, C. Wu, S. Guo and J. Zhang, Nanotechnol. Rev., 2013, 2, 27–45.
9. S. Singla, E. Anim-Danso, A. E. Islam, Y. Ngo, S. S. Kim, R. R. Naik and A. Dhinojwala, ACS Nano, 2017, 11, 4899–4906.
10. E. Voloshina, D. Usvyat, M. Schütz, Y. Dedkov and B. Paulus, Phys. Chem. Chem. Phys., 2011, 13, 12041.
11. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
12. R. R. Freitas, R. Rivelino, F. D. B. Mota and C. M. De Castilho, J. Phys. Chem. A, 2011, 115, 12348–12356.
13. J. Klime and A. Michaelides, J. Chem. Phys., 2012, 137, 120901.
14. P. Singla, M. Riyaz, S. Singhal and N. Goel, Phys. Chem. Chem. Phys., 2016, 18, 5597–5604.
15. W. F. van Gunsteren and H. J. C. Berendsen, Angew. Chem., Int. Ed. Engl., 1990, 29, 992–1023.
16. A. Amirfazli and A. W. Neumann, Adv. Colloid Interface Sci., 2004, 110, 121–141.
17. W. L. Jorgensen and J. Tirado-Rives, J. Am. Chem. Soc., 1988, 110, 1657–1666.
18. W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell and P. A. Kollman, J. Am. Chem. Soc., 1995, 117, 5179–5197.
19. A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin and M. Karplus, J. Phys. Chem. B, 1998, 102, 3586–3616.
20. C. Oostenbrink, A. Villa, A. E. Mark and W. F. Van Gunsteren, J. Comput. Chem., 2004, 25, 1656–1676.
21. H. Yu, T. Hansson and W. F. Van Gunsteren, J. Chem. Phys., 2003, 118, 221–234.
22. A. P. E. Kunz and W. F. Van Gunsteren, J. Phys. Chem. A, 2009, 113, 11570–11579.
23. I. V. Leontyev and A. A. Stuchebrukhov, J. Chem. Theory Comput., 2012, 8, 3207–3216.
24. P. E. Lopes, J. Huang, J. Shim, Y. Luo, H. Li, B. Roux and A. D. Mackerell, J. Chem. Theory Comput., 2013, 9, 5430–5449.
25. S. J. Bachmann and W. F. Van Gunsteren, J. Chem. Phys., 2014, 141, 22D515.
26. K. Vanommeslaeghe and A. D. Mackerell, Biochim. Biophys. Acta, Gen. Subj., 2015, 1850, 861–871.
27. Y. Shi, P. Ren, M. Schnieders and J. P. Piquemalc, Rev. Comput. Chem., 2015, 2, 51–86.
28. J. A. Lemkul, J. Huang, B. Roux and A. D. MacKerell, Chem. Rev., 2016, 116, 4983–5013.
29. A. P. E. Kunz, J. R. Allison, D. P. Geerke, B. A. C. Horta, P. H. Hünenberger, S. Riniker, N. Schmid and W. F. Van Gunsteren, J. Comput. Chem., 2012, 33, 340–353.
30. T. P. Straatsma and J. a. McCammon, Mol. Simul., 1990, 5, 181–192.
31. P. Drude, The Theory of Optics, Longsmans, 1901.
32. M. Born and K. Huang, Dynamical Theory of Crystal Lattices, 1954.
33. H. Yu and W. F. Van Gunsteren, J. Chem. Phys., 2004, 121, 9549–9564.
34. M. C. Gordillo and J. Martí, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 1–5.
35. T. a. Ho and A. Striolo, Mol. Simul., 2014, 7022, 1–11.
36. C. A. Jimenez-Cruz, S.-g. Kang and R. Zhou, Wiley Interdiscip. Rev.: Syst. Biol. Med., 2014, 6, 329–343.
37. J. Liu, Z. Yang, H. Li, Z. Gu, J. A. Garate and R. Zhou, J. Chem. Phys., 2014, 141, 22D520.
38. V. H. Rusu, S. Bachmann and W. F. van Gunsteren, Mol. Phys., 2016, 114, 845–854.
39. G. Lamoureux, A. D. MacKerell and B. Roux, J. Chem. Phys., 2003, 119, 5185–5197.
40. T. A. Ho and A. Striolo, J. Chem. Phys., 2013, 138, 054117.
41. E. R. Cruz-chu, A. Aksimentiev and K. Schulten, J. Phys. Chem. B, 2006, 110, 21497–21508.
42. T. Ingebrigtsen and S. Toxvaerd, J. Phys. Chem. C, 2007, 111, 8518–8523.
43. T. Werder, J. H. Walther, R. L. Jaffe, T. Halicioglu, P. Koumoutsakos and V. Sunny, J. Phys. Chem. B, 2003, 107, 1345–1352.
44. A. Ashraf, Y. Wu, M. C. Wang, N. R. Aluru, S. A. Dastgheib and S. W. Nam, Langmuir, 2014, 30, 12827–12836.
45. A. Kozbial, C. Trouba, H. Liu and L. Li, Langmuir, 2017, 33, 959–967.
46. Z. Li, Y. Wang, A. Kozbial, G. Shenoy, F. Zhou, R. McGinley, P. Ireland, B. Morganstein, A. Kunkel, S. P. Surwade, L. Li and H. Liu, Nat. Mater., 2013, 12, 925–931.
47. D. R. Tobergte and S. Curtis, Springer Ser. Surf. Sci., 2013, 53, 1689–1699.
48. J. Włoch, A. P. Terzyk and P. Kowalczyk, Chem. Phys. Lett., 2017, 674, 98–102.
49. B. M. Law, S. P. McBride, J. Y. Wang, H. S. Wi, G. Paneru, S. Betelu, B. Ushijima, Y. Takata, B. Flanders, F. Bresme, H. Matsubara, T. Takiue and M. Aratono, Prog. Surf. Sci., 2017, 92, 1–39.
50. F. Leroy, S. Liu and J. Zhang, J. Phys. Chem. C, 2015, 119, 28470–28481.
51. V. Sresht, A. Govind Rajan, E. Bordes, M. S. Strano, A. A. Pádua and D. Blankschtein, J. Phys. Chem. C, 2017, 121, 9022–9031.
52. R. P. Misra and D. Blankschtein, J. Phys. Chem. C, 2017, 121, 28166–28179.
53. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graph., 1996, 14, 33–38.
54. H. Berendsen, J. Postma, W. van Gunsteren and J. Hermans, Intermolecular Forces, Dordrecht, Germany, Pullman. B, Reidel, 1981, pp. 331–342.
55. N. Schmid, C. D. Christ, M. Christen, A. P. Eichenberger and W. F. Van Gunsteren, Comput. Phys. Commun., 2012, 183, 890–903.
56. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
57. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
58. I. G. Tironi, R. Sperb, P. E. Smith and W. F. van Gunsteren, J. Chem. Phys., 1995, 102, 5451–5459.
59. B. Roux, Biophys. J., 2008, 95(9), 4205–4216.
60. J. A. Garate, T. Perez-Acle and C. Oostenbrink, Phys. Chem. Chem. Phys., 2014, 16, 5119–5128.
61. A. P. Eichenberger, J. R. Allison, J. Dolenc, D. P. Geerke, B. A. Horta, K. Meier, C. Oostenbrink, N. Schmid, D. Steiner, D. Wang and W. F. Van Gunsteren, J. Chem. Theory Comput., 2011, 7, 3379–3390.
62. B. Mukherjee, P. K. Maiti, C. Dasgupta and A. K. Sood, J. Chem. Phys., 2007, 126, 124704.
63. J. Garate, A. Bernardin, Y. Escalona, C. Yanez, N. English and T. Perez-Acle, J. Phys. Chem. B, 2019, 123, 2599–2608.
64. J. Kofinger, G. Hummer and C. Dellago, Phys. Chem. Chem. Phys., 2011, 13, 15403–15417.
65. J.-A. Garate, N. English and J. MacElroy, Molecular Simulation, Taylor & Francis, 2009, vol. 35, pp. 3–12.
66. A. Bateni, S. Laughton, H. Tavana, S. S. Susnar, A. Amirfazli and A. W. Neumann, J. Colloid Interface Sci., 2005, 283, 215–222.
67. V. Vancauwenberghe, P. D. Marco and D. Brutin, Colloids Surf., A, 2013, 432, 50–56.
68. H. Li and X. C. Zeng, ACS Nano, 2012, 6, 2401–2409.
69. J. Rafiee, X. Mi, H. Gullapalli, A. V. Thomas, F. Yavari, Y. Shi, P. M. Ajayan and N. A. Koratkar, Nat. Mater., 2012, 11, 217–222.
70. C. J. Shih, Q. H. Wang, S. Lin, K. C. Park, Z. Jin, M. S. Strano and D. Blankschtein, Phys. Rev. Lett., 2012, 109, 176101.
71. C.-J. Shih, M. S. Strano and D. Blankschtein, Nat. Mater., 2013, 12, 866–869.
72. F. Taherian, V. Marcon, N. F. A. Van Der Vegt and F. Leroy, Langmuir, 2013, 29, 1457–1465.
73. Q. Liu and B. Xu, AIP Adv., 2015, 5, 67150.
74. C. T. Nguyen and B. H. Kim, Int. J. Precis. Eng. Manufact., 2016, 17, 503–510.
75. J. Włoch, A. P. Terzyk, P. A. Gauden, R. Wesołowski and P. Kowalczyk, J. Phys.: Condens. Matter, 2016, 28, 495002.
76. J. E. Andrews, S. Sinha, P. W. Chung and S. Das, Phys. Chem. Chem. Phys., 2016, 18, 23482–23493.
77. A. V. Prydatko, L. A. Belyaeva, L. Jiang, L. M. C. Lima and G. F. Schneider, Nat. Commun., 2018, 9, 4185.
78. I. Szleifer and B. Widom, Mol. Phys., 2006, 75, 925–943.
79. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys. Chem., 1987, 91, 6269–6271.
80. C. Vega and E. de Miguel, J. Chem. Phys., 2007, 126, 154707.
81. D. Parobek and H. Liu, 2D Mater., 2015, 2, 032001.
82. G. Hong, Y. Han, T. M. Schutzius, Y. Wang, Y. Pan, M. Hu, J. Jie, C. S. Sharma, U. Müller and D. Poulikakos, Nano Lett., 2016, 16, 4447–4453.
83. W. Jorgensen, J. Chandrasekhar, J. Madura, R. Impey and M. Klein, J. Chem. Phys., 1983, 79, 926–935.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1cp05573a

---