Solvent-mediated interactions between nanostructures: from water to Lennard-Jones liquid

Solvent-mediated interactions emerge from complex mechanisms that depend on the solute structure, its wetting properties and the nature of the liquid. While numerous studies have focused on the two first influences, here, we compare results from water and Lennard-Jones liquid in order to reveal to what extent solvent-mediated interactions are universal with respect to the nature of the liquid. Besides the influence of the liquid, results were obtained with classical density functional theory and brute-force molecular dynamics simulations which allows us to contrast these two numerical techniques.


I. INTRODUCTION
Interactions between two solids are usually well characterized by their intrinsic physical and chemical properties.
Numerous studies have examined solvent-mediated forces in terms of range, strength and sign using both numerical [10][11][12][13][14][15][16][17][18][19] and experimental techniques [20][21][22][23][24][25][26][27][28][29][30] .In particular, it was found that the sign of solventmediated forces is controlled by the equilibrium contact angle [15][16][17][18]21 . On he one hand, when the solids are solvophilic, solvent molecules are attached to the solid surface.Bringing the two solutes together leads to a perturbation and the removal of this favorable structure causes a strong repulsive hydration pressure [28][29][30] .On the other hand, when the solids are solvophobic, solvent molecules which are located between the two solids are expelled and a more stable vapor cavity emerges thus reducing the overall free energy.This so-called capillary evaporation has been the subject of numerous works [31][32][33][34][35][36] and leads to a solvophobic attraction.Along with the wetting properties of the solid, the role of its geometrical structure was also covered in several studies [37][38][39] .For example, Jabes et al. recently demonstrated that using the same solid composition and size, qualitatively different solvent-mediated forces can be obtained only by changing the solid shape between fullerenes, nanotubes and graphene-like structures 38 .Altogether, this work is integrated to a general mean-field theory of hydrophobicity developed by Lum, Chandler and Weeks 40 and further refined in more recent publications 31,41,42 .When modeling liquids using numerical simulations, two approaches are commonly employed.On the one hand, because water plays a crucial role in most applications and especially in biological systems, atomistic models for water molecules have been developed in order to mimic its thermodynamic properties and some special features including strong hydrogen bonda) Electronic mail: julien.lam@ulb.ac.be ing and ice polymorphism.On the other hand, generic model systems such as hard spheres [43][44][45][46][47][48][49][50][51] and Lennard-Jones potential 17,52,53 are also often used for modeling fluids.Once the model is chosen, solvent-mediated forces can be computed using various types of numerical methods.Monte-Carlo and molecular dynamics simulations are widely employed especially for water modeling while classical density functional theory (DFT) in which molecules are treated as a density field can grant access directly to liquid density profiles and the corresponding free energy 18,54 . DFT i less numerically expensive and avoids using free energy calculation techniques such as thermodynamic integration, transition path sampling and umbrella sampling.However, DFT for water is not as highly developed as for simple fluids [55][56][57] .
While numerous authors have suggested the ability of the Lennard-Jones liquid to reproduce behaviors similar to water regarding solvent-mediated effects [58][59][60] , there is no detailed comparison of solvent-mediated forces obtained with atomistic simulations of extended simple point charge model water and with DFT calculations of Lennard-Jones particles (LJ).In this work, we make a direct comparison of solvent-mediated forces obtained from molecular dynamics simulations of water and DFT calculations of LJ using a very generic system made of two nanometric crystalline slabs immersed in a liquid.Free energy is computed as a function of the interslab distance and we study thoroughly the influence of wettability and of the slab geometrical structure.Our work identifies differences and similarities between atomistic simulations of water and DFT calculations of LJ.Moreover, our results also contribute to the overall understanding of solventmediated forces and discuss more generally to what extent molecular properties of water make it special in comparison to simple fluid models.

A. Studied system
Our calculations make use of two types of molecules: slab molecules and liquid molecules.The slabs are composed of rigid arrangements of solid molecules while the liquid is treated dynamically.We held constant temperature and density of the liquid while varying the solid properties.The slabs are made of three square layers of 40 + 41 + 40 = 121 atoms which are kept fixed in a face centered cubic structure with the (100) face exposed and with the lattice spacing, a.The interaction between slab and liquid particles is modeled via a Lennard-Jones potential parametrized by its length scale, σ wall and well depth, wall .When varying a, σ wall is also modified using: σ wall = a a0 σ liq where a 0 is the zerotemperature FCC equilibrium lattice spacing equal to 1.5424 σ liq 61 .With such a model, wetting properties as defined by the contact angle are driven by the ratio between the liquid/liquid and liquid/solid attractions.In practice, we varied wall while keeping the liquid properties constant and measured the corresponding contact angle, θ.We note that additional complexities that also influence the surface solvophobicity including functionalization and polarity effects can not be captured with our present model 38 .
Finally, results are shown in physical units.For water, energy is displayed in kcal/mol and distances in angstroms.When computing a, we used σ liq = 2.75 Å as it is the approximate size of a water molecule.For LJ, the potential parameters are denoted LJ and σ LJ .We worked at a temperature of k B T = 0.8 LJ and we chose: k B T = 0.593 kcal/mol at 300 K in order to rescale LJ to real units.Concerning the distances, we imposed σ LJ = 2.75 Å as well.

B. Molecular dynamics simulation of water
The extended simple point charge model (SPC/E) is used for water 62 .Bonds in water molecules are constrained using the SHAKE algorithm and long-range Coulombic interactions are computed with the Particle-Particle-Particle-Mesh solver with a precision tolerance equal to 10 −4 and a real space cutoff equal to 9.8 Å.At the initialization step, water molecules are arranged on a simple cubic lattice structure with a lattice spacing equal to 3.1 Å. Solids are modeled with rigid molecules made of three face-centered cubic layers.Solid molecules and oxygen atoms of water interact via a truncated and shifted Lennard-Jones (LJ) potential with a cutoff equal to 9.8 Å.The LAMMPS package 63 is used for all the simulations.From there, two types of calculations are performed: (i) Droplet equilibration to measure the wetting properties of the solid and (ii) Solvent-mediated forces between two slabs.

Droplet equilibration and contact angle
A hemisphere of water with radius equal to 50 Å is initially deposited onto the solid surface.On top of the spherical cap, a cage made of fixed atoms is also placed to help the droplet equilibration and prevent it from leaving the solid surface at the initialization stage.These atoms only interact with oxygen atoms via a LJ potential ( cage = 1 kcal/mol and σ cage = σ O−O = 3.166 Å).The entire simulation box measures 300 Å×300 Å×200 Å which is large enough to avoid the influence of periodic images.For the equilibration protocol, the timestep is set equal to 0.5 fs.NVE simulations are performed during 5 ps then NVT simulations are performed during another 5 ps at 300 K. From there, cage atoms are removed to allow for droplet shape relaxation and the timestep is changed to 2 fs.After an equilibration run during 500 ps, snapshots are taken every 2.5 ps during 500 ps.Density profiles of oxygen atoms are averaged through time [See Fig. 1].From the density profiles, a liquid/gas interface is obtained as depicted in Fig. 1.A linear fit with all the points located below 8 Å is then used to compute the contact angle [See Fig 1].Uncertainties are evaluated by the standard deviation measured every 100 ps for 5 independent runs then an additional factor of two is incorporated to account for error in the method for contact angle extraction.In addition, the duration of the simulation is sufficient to reach equilibration as observed in Fig. 4.a.

Calculation of the solvent-mediated interactions
Two nanoslabs which are made of 40 + 41 + 40 = 121 atoms are positioned parallel to each other.The entire simulation box measures 52 Å×52 Å×60 Å and contains 5344 water molecules [See Fig. 5].The solids are first disposed on top of each other and NVT simulations are performed at 300 K during 60 ps with a timestep equal to 2 fs.After this equilibration procedure, the solids are instantaneously moved apart by 0.25 Å with the distance measured as the difference in height between the center of mass of both slabs.For each separation denoted z, the system is equilibrated for 50 ps and production run is done during another 50 ps.The free energy as a function of z is then given by numerical integration of the forces: where f 1 and f 2 are the forces between water molecules and solid atoms with 1 and 2 designating respectively the upper and the lower solids.u z is a unit vector along the z direction going upward.The difference in forces used in the integration scheme Eqn.II B 2 is shown in

C. Density functional theory calculation of Lennard-Jones
For this second method, liquid particles interact via a LJ potential with LJ and σ LJ as energy and length parameters.The cutoff distance is equal to 3σ LJ .The density and the temperature are respectively ρ LJ = 0.7 σ −3 LJ and k B T = 0.8 LJ which is located between the triple point and the critical temperature.This corresponds to the liquid density for a chemical potential supersaturation equal to ∆µ = 0.27k B T 18 .While the value of ∆µ influences quantitatively the solvent-mediated forces 17 , the supersaturation is chosen in this work to match the ratio of pressure between water coexistence pressure and atmospheric pressure (P coex = 1/20P atm ).Interactions between liquid molecules and solid atoms are also modeled with LJ potential truncated at 3σ LJ .Within the DFT framework, free energy is expressed as a functional of the liquid density.For LJ interaction, the potential is separated in two parts, the repulsive part modeled with the White Bear functional 64 and the attractive part treated in mean field.The density is computed on a discretized three dimensional grid with 8 lattice points per unit of σ LJ and the free energy is obtained through minimization with respect to the density field.In order to match MD calculations that are made in the NVT ensemble, DFT calculations are also run with a fixed number of particles rather than a fixed chemical potential.The DFT method is described in greater details in our previous contributions 18,54,65 .Accuracy of the DFT treatment is discussed in this review 54 .From there, droplet equilibration results were taken from our previous work 18 .For solvent-mediated interactions, we used the same system as with molecular dynamics simulation of water except that there is no equilibration protocol and the free energy is obtained directly through DFT.
In recent works regarding solvent mediated forces, calculations are performed in µVT 17 or NPT 13,38 ensembles in order to supply particles during the drying transition.To evaluate if our system is large enough to cope with this issue, we also ran calculations in the µVT with ∆µ = 0.27k B T [See Fig. 8.b].Results are not significantly different from those obtained in NVT thus justifying our approach.

III. RESULTS AND DISCUSSION
Figures 5 and 6 show typical results obtained respectively for water and Lennard-Jones.
In both cases, when walls are solvophilic, the gap between the two slabs is filled with liquid even at small distances [See Figs.(5b, 6b)].For solvophobic walls, this happens only for large enough separations.In addition, structuring can be observed in the vicinity of the slabs especially when looking at the density profiles obtained by DFT of Lennard-Jones particles [See Fig. 6].The structure is more pronounced for solvophilic walls.

A. Influence of the wall lattice spacing
In this first study, we worked at a fixed value of wall while changing the wall lattice spacing so that both the structure and the hydrophobicity are modified.In Figs .(7a, 7b), excess free energy is plotted for a moderate value of wall (0.1 kcal/mol and 0.2 LJ ) that in both cases lead to 35 • when a = a 0 .An almost linear decrease is observed which is consistent with previous works on solvophobic attraction 10,18,38 .In particular, near contact, one can show that the slope depends solely on bulk properties using a capillary model 17 .At intermediate distances, the slope is also influenced by the wall solvophobicity since the presence of a meniscus leads to a non trivial shape of the gaseous phase 17,18 .When a is reduced, the walls are denser and thus becomes less solvophobic.Therefore, both the range and the height of the solvent-mediated interaction are reduced.The results obtained with LJ DFT and with water MD are qualitatively similar and we demonstrate that in this case, LJ can be used to reproduce water-mediated interactions.For solvophilic walls, the results of the comparison are not so close [See  Figs .(7c, 7d)].In both systems, oscillations in the free energy are observed due to layering of the liquid near the walls.However, the oscillation amplitudes vary significantly between water and the LJ fluid.This is likely due to the asymmetry of water molecules: as they pack together to form denser layers near the wall, their interlayer distance does not depend solely on their average size but rather on their size in some particular directions 14 .

B. Influence of the wall energy
Solvent-mediated interactions are plotted at a fixed value of a = a 0 but for different values of wall in Fig. 8.When comparing results from MD water and DFT LJ, several similarities can be identified.First, when the walls are solvophobic, free energy monotonically increases as the nanoslabs are pulled apart with an almost linear behavior.Then, when the walls are solvophilic, damped oscillations are observed because of the emergence of structured layers near the wall.Also, the lowest energy state is always at contact meaning that the nanoparticle would preferentially stay near the wall as long as it overcomes the intermediate free energy barrier.Finally, the range of the depletion force does not go beyond 20 Å which corresponds to approximately 7 liquid layers.These similarities were already raised in the literature [58][59][60] and our work allows for a more direct comparison as we studied the same system (ie.two nanoslabs made of the same structure) while only changing the liquid nature.The solvophobicity is increased as coloring goes from blue to red with wall going from 0.50 kcal/mol (0.65 LJ ) to 0.05 kcal/mol (0.05 LJ ) for water (for LJ).This corresponds to contact angle ranging from 0 • to 180 • .Each color is separated by 0.05 kcal/mol (by 0.1 LJ ) for water (for LJ).Dash lines correspond to results obtained with LJ/DFT in the µVT ensemble using ∆µ = 0.27kBT .
In order to quantitatively compare results from water MD and LJ DFT, we define two positions: (i) z = 5 Å gives ∆F contact and (ii) the position, denoted z mid , at which the most solvophilic interaction reaches its maximum is used as an intermediate value called ∆F mid [See Fig. 8].In Fig. 9, results are reported for different contact angles that are determined after the equilibration of sessile drops.For the highest degrees of solvophilicity (and solvophobicity), droplets are not stable and the contact angle is trivially 0 • (and 180 • ).Therefore, when reporting ∆F mid and ∆F contact as a function of the corresponding contact angle, not all the data from Fig. 8 are considered.As the contact angle is increased, ∆F mid decreases almost linearly.Water MD and LJ DFT curves have similar slopes and the constant difference between the two curves is roughly of 24 kcal/mol.However, the sign of ∆F mid is different which indicates that qualitatively different behaviors are expected.In the water case, this intermediate distance is less energetically favorable than having the nanoslabs far from each other.For ∆F contact , LJ DFT and water MD also lead to qualitatively different results.Indeed, while for LJ, ∆F contact , like ∆F mid , decreases linearly, for water, ∆F contact is non-monotonic and peaks around 80 • .As already raised in the previous section, LJ is not well-adapted to model water at contact because water has orientational order especially for solvophilic walls which can not be seen with LJ.Furthermore, while there is an intermediate range of solvophibicity (θ ∈ [70 • : 110 • ]) where a good agreement for ∆F contact is found, the signs of ∆F mid are different as raised above.

C. Hysteresis in solvent-mediated forces
Throughout this study, results on the solvent-mediated forces were obtained as the two solutes are pulled apart from each other.Yet, another possibility concerns the case when the slabs are disposed far from each other and then brought together.This leads to the question of reversibility of the interaction.In Fig. 10, the solvent mediated forces are plotted in this second approach.In the case of water MD, qualitative agreement is found when comparing results from Fig. 8 with solvophobic walls.However, strong repulsive interactions are observed with solvophilic walls.This results from water molecules that can be trapped between the two plates if they are brought together too rapidly.Furthermore, the capillary evaporation which is not observed when the solutes are pulled apart, is not only driven by the solute interdistance and additional order parameters such as the solvent density between the solutes can be used [32][33][34] .Essentially, the time for gas to nucleate between the two walls is so large that we can not obtain the equilibrium state with brute-force molecular dynamics simulations 32 .This apparent hysteresis is not found in LJ DFT since the technique enables to circumvent any of these kinetic issues and leads directly to the most stable state in which the gap between the walls is emptied of liquid.

IV. CONCLUSIONS
In summary, solvent-mediated forces were measured in a very generic case of two nanoslabs embedded in a liquid.Two different models for the liquid along with two different methods for measuring free energy were employed.Such a direct comparison of these two approaches allowed us to identify both similarities and differences.On the one hand, oscillations for solvophilic walls and a linear decrease for solvophobic walls were observed with the two liquids.In addition, the range of the depletion force and the presence of a minimum at contact are two additional features that seem to support the idea of a universal behavior of the solvent-mediated forces.On the other hand, no region of solvophobicity seems to show a quantitative agreement between water and LJ.In particular, amplitudes of the oscillations and the resulting sign of the free energy for intermediate distances are different.Also, the value of the free energy at contact does not have the same behavior as the contact angle is changed.Ultimately, using LJ or water in order to model solvent mediated forces should depend on the desired level of accuracy and our results provide a benchmark that quantifies the error made if one wishes to use LJ instead of water.

Figure 1 .
Figure 1.Water density profiles for an equilibrated water droplet on top of a wall with wall equal to 0.1 kcal/mol (a), 0.2 kcal/mol (b) and 0.3 kcal/mol (c).Red dots represent the liquid-gas interface which is used to compute the contact angle and red lines are the associated linear fit.The image size is equal to 30 Å×60 Å.

Figure 2 .
Figure 2. Contact angle as a function of the wall energy depth obtained with water/MD (a) and LJ/DFT (b).Dotted lines are obtained through linear fitting.

Fig 3 .
Duration of the simulation time is considered sufficient to reach equilibration as assessed by Fig. 4.b.Error bars in this figure are computed as the standard deviation obtained with 5 independent runs.

Figure 3 .
Figure 3. Solvophobicity influence on the forces obtained at a = a0 with SPC/E Water.The solvophobicity is increased as coloring goes from blue to red with wall going from 0.50 kcal/mol to 0.05 kcal/mol.This corresponds to contact angle ranging from 0 • to 180 • .Each color is separated by 0.05 kcal/mol.

Figure 4 .
Figure 4. Influence of equilibration time over (a) the contact angle and (b) the excess free energy for different degrees of hydrophobicities.For contact angle calculations, red and blue data are obtained with wall respectively equal to 0.42 kcal/mol and 0.1 kcal/mol.The simulation time employed in this work is represented with black dots (1000 ps for contact angle) and with plain lines (100 ps for free energy calculation).

Figure 5 .
Figure 5. Liquid density profiles for water molecules confined between two nanoslabs obtained for a = a0 at different distances and degrees of solvophobicity.Contact angles of 35 • and 120 • are obtained with wall respectively equal to 0.42 kcal/mol and 0.10 kcal/mol.Each image is a slice of width 1 Å passing through the middle of the solute and measuring 48 Å×48 Å.

Figure 6 .
Figure 6.Liquid density profiles for Lennard-Jones particles confined between two nanoslabs obtained for a = a0 at different distances and degrees of solvophobicity.Contact angles of 35 • and 120 • are obtained with wall respectively equal to 0.4 LJ and 0.2 LJ .Each image is a slice of width 0.275 Å passing through the middle of the solute and measuring 68.75 Å×68.75 Å.

Figure 7 .
Figure 7. Lattice spacing influence on the excess free energy obtained with SPC/E Water (left) and with Lennard-Jones (right).Calculations are run at a value of wall for which the contact angle is θ = 120 • (top) and at θ = 35 • (bottom) for a = a0.Values of wall are given in captions of Fig.5 and Fig.6.

Figure 8 .
Figure 8.Wall energy influence on the excess free energy obtained at a = a0 with SPC/E Water (a) and with Lennard-Jones (b).The solvophobicity is increased as coloring goes from blue to red with wall going from 0.50 kcal/mol (0.65 LJ ) to 0.05 kcal/mol (0.05 LJ ) for water (for LJ).This corresponds to contact angle ranging from 0 • to 180 • .Each color is separated by 0.05 kcal/mol (by 0.1 LJ ) for water (for LJ).Dash lines correspond to results obtained with LJ/DFT in the µVT ensemble using ∆µ = 0.27kBT .

Figure 9 .
Figure 9.Comparison between results of LJ DFT and water MD using the contact angle dependence of ∆F mid (a) and ∆Fcontact (b).In grey, results for ∆F mid obtained with LJ DFT is vertically shifted in order to show that the difference with water MD is only a constant.

Figure 10 .
Figure 10.Solvent-mediated forces at a = a0 with SPC/E Water (a) and with Lennard-Jones (b) obtained as the solutes are brought closer to each other.Color designations are described in Fig. 8.
The work of JL was funded by the European Union's Horizon 2020 research and innovation program within the AMECRYS project under grant agreement no.712965.JFL thanks the European Space Agency (ESA) and the Belgian Federal Science Policy Office (BELSPO) for their support in the framework of the PRODEX Programme, contract number ESA17 AO-2004-070.Computational resources have been provided by the Consortium des Equipements de Calcul Intensif (CECI) and by the Federation Lyonnaise de Modelisation Sciences Numeriques (FLMSN).