Large effect of lateral box size in molecular dynamics simulations of liquid-solid friction

Molecular dynamics simulations are a powerful tool to characterize liquid-solid friction. A slab conﬁguration with periodic boundary conditions in the lateral dimensions is commonly used, where the measured friction coefﬁcient could be affected by the ﬁnite lateral size of the simulation box. Here we show that for a very wetting liquid close to its melting temperature, strong ﬁnite size effects can persist up to large box sizes along the ﬂow direction, typically ∼ 30 particle diameters. We relate the observed decrease of friction in small boxes to changes in the structure of the ﬁrst adsorbed layer, which becomes less commensurable with the wall structure. Although these effects disappear for lower wetting cases or at higher temperatures, we suggest that the possible effect of the ﬁnite lateral box size on the friction coefﬁcient should not be automatically set aside when exploring unknown systems.


I. INTRODUCTION
Nanofluidic systems [1,2] show great promise for applications such as blue energy [3][4][5] and water desalination [6][7][8]. When liquids are confined at nanometric scales, surfaces play an increasingly important role over bulk liquid properties, and it is therefore critical to better understand and control the behavior of liquid-solid interfaces. In particular, both experiments and molecular simulations have shown that liquids can slip on some surfaces at the nanoscale [9], which can enhance the performance of nanofluidic systems [10][11][12][13][14]. Wall slip is controlled by the liquid-solid friction, as initially discussed by Navier [15,16], who wrote that the viscous shear stress in the liquid at the wall η∂ z v| z=z w (where η is the shear viscosity, z is the direction normal to the interface, z w is the wall position, and v is the tangential velocity) is equal to an interfacial fluid friction stress, τ = λv s (where v s is the slip velocity, the velocity jump at the interface, and λ is the fluid-solid friction coefficient). This relation is different from the Coulomb friction coefficient defined as the ratio of the friction force to the normal force, usually applied in the context of solid friction. Navier's boundary condition can be rewritten as the so-called partial slip boundary condition, v s = b ∂ z v, defining the slip length b, which can be expressed as a function of the liquid viscosity and the interfacial friction Because standard slip lengths lie in the nanometric range, experimental characterization of slip is very challenging [17][18][19][20]. Molecular dynamics (MD) simulations of liquids confined between parallel walls (slab configuration) or in cylindrical pores represent an interesting alternative for slip or friction characterization , offering good control of the simulated setups. Moreover, MD can help us to understand the mechanisms underlying interfacial friction by providing access to atomistic information [42][43][44][45][46][47][48][49], and quantum effects on interfacial friction were also examined recently through ab initio simulations with massive computing power [47]. However, MD simulations also pose a number of specific issues. For instance, in direct simulations of transport properties by nonequilibrium MD, the thermostating algorithm required to maintain the system at constant temperature can affect the results and must be carefully chosen [50][51][52].
Simulations are also limited by the finite size of the simulation cells. Although the use of periodic boundary conditions can, in some aspects, efficiently mimic a macroscopic system, finite size effects can still impact the measured quantities. For instance, the self-diffusion coefficient strongly depends on the size of the simulation box because of hydrodynamic interactions with periodic images, an effect well captured by continuum hydrodynamics [53][54][55][56]. In slab simulations, two types of finite size effects should be considered. On the one hand, the effect of the confinement height has been explored in previous work [27][28][29]36,37,[43][44][45]. On the other hand, the effect of the lateral size of the simulation box has been little explored. Huang and Szlufarska [35] compared the predictions of two methods to extract the friction coefficient from equilibrium simulations, using different lateral sizes. However, we are not aware of tests of the effect of the lateral box size on the friction coefficient, regardless of the measurement method. Previous MD work on liquid-solid friction cited above typically used lateral box sizes ranging (a) Density (black line) and velocity (red dots) profiles for three lateral box sizes, l x = 3.13, 6.27 and 18.8 nm; the red dotted line is a linear fit of the velocity profile in the bulk region. As illustrated in the three panels, the slip velocity v s is measured as the liquid velocity at the position of the first density peak. (b) Snapshots of the systems, from left to right: perspective view of the small system and side views of three systems with increasing lateral size l x . from ∼5σ to ∼50σ , with σ being the molecular size. In this work we would like to assess the validity of such box sizes by examining finite size effects on liquid-solid friction. To that aim we will consider a model Lennard-Jones (LJ) system and use nonequilibrium shear flow simulations. We will show that significant finite size effects can persist up to large box sizes for very wetting liquids close to the melting temperature, and we will show the atomic origin of these effects.

II. SYSTEMS AND METHODS
We simulated a liquid slab confined between parallel walls [see Fig. 1(b)]. The walls were fcc crystals, each made of eight atomic layers and exposing a (001) face to the liquid. Every pair of nearest neighbors in the walls was bound through a harmonic potential h (r) = k/2 (r − r eq ) 2 , with r being the interparticle distance, r eq = 0.277 nm, and k = 46.8 N/m. Interactions between liquid particles and between liquid and solid particles were modeled by a LJ potential LJ 6 ], with i and j being L for liquid particles and S for solid ones, where the LJ interaction was truncated at a cutoff distance of r c = 3.5σ and quadratic functions were added so that the potential and interaction force smoothly vanished at r c [57]. We used the following parameters: σ LL = 0.340 nm, σ LS = 0.345 nm, ε LL = 121 K × k B , and ε LS = α ε LL , where the "wetting coefficient" α controls the wetting properties and can be related to the contact angle (see Appendix A). We used α = 0.774, corresponding to a complete wetting case, unless otherwise specified. Finally, the atomic masses were m L = 39.95 u and m S = 195.1 u.
The total system height (walls included) and the total liquid height along the z direction were ∼10.8 and 8 nm, respectively, to ensure that an ∼2-nm-thick region of bulk liquid (with a constant density) existed in the center between the layered liquid structures formed near the two walls. We used periodic boundary conditions along the x and y lateral directions. The box size along the y direction l y was 3.13 nm, and we varied the size along the x direction l x between 3.13 and 25.04 nm. To that aim, we simply replicated the smallest system along the x direction, so that for all systems the number of liquid particles per wall surface area was strictly identical.
The liquid was sheared by imposing opposite velocities of ±10 m/s along the x direction to wall particles belonging to the outmost layer of both walls. The temperature was set to 85 K = 0.703ε LL /k B , i.e., very close to the melting point. To that aim we applied a Langevin thermostat to wall particles belonging to the second outmost layer, only to the y degrees of freedom perpendicular to confinement and flow in order to minimize the impact of the thermostating procedure on the dynamics of the liquid and of the liquid-solid interface. In particular, with the present cutoff distance, liquid particles do not directly interact with the thermostated wall particles. We checked that the liquid temperature remained within 1 K of the set value. To control the system pressure at 4 MPa, we used the following procedure: first, we used the top wall as a piston, without shear; second, we applied shear, still using the top wall as a piston, and finally, we fixed the vertical wall position and continued shearing the system. The equations of motion were integrated using the velocity-Verlet algorithm, with a time step of 5 fs. The production run lasted between 50 and 250 ns depending on the measured quantity. Figure 1(a) shows the density and velocity profiles obtained for different lateral sizes l x . The density profiles are independent of the box size and indicate that for the temperature and the liquid-solid interaction energy we chose, the liquid is strongly layered at the walls: the maximum density in the first adsorption layer (hereafter called the contact layer) is ∼7 times the bulk value, and at each surface around eight peaks can be observed before the density reaches the bulk value. Note that because we used the top wall as a piston during the preparation stage, the position at which the top wall was fixed during the production stage varied slightly from one simulation to the other. However, the resulting difference was negligibly small, e.g., about 4×10 −3 nm between systems with l x = 3.13 nm and l x = 18.8 nm.

III. RESULTS AND DISCUSSION
The shear velocity profiles are linear throughout almost all the liquid's thickness, except in the contact layer, where the shear rate slightly decreases, which can be understood as a small increase of the local viscosity due to the large local density [58,59]. In this work, we determined the slip velocity v s as the difference between the wall velocity and the liquid velocity at the position of the first adsorption peak. In Fig. 1 it is already clear that v s changes with the box lateral size. In order to extract the friction coefficient from the simulations, we also measured the friction stress τ as the friction force F x between the liquid and the wall (top or bottom) divided by the interface area: τ = F x /(l x l y ), where F x was calculated as the total force in the x direction between wall particles of all atomic layers and fluid particles. Note that with the present cutoff distance mentioned in Sec. II, there was no direct interaction between liquid particles and thermostated wall particles or with outermost wall particles moving at a constant velocity. The friction coefficient was then given by λ = τ /v s . Figure 2 displays the relation between the lateral size l x and friction coefficient λ. It shows that the system lateral size has a dramatic impact on the friction coefficient for the present system of interest at T = 85 K and α = 0.774: first, the amplitude of the effect is large (λ decreases by a factor of ∼6 in the smallest box), and second, the effect remains significant up to quite large simulation boxes (∼10 nm, i.e., ∼30 particle diameters). As detailed later, this dependence disappeared for systems with a smaller wetting coefficient or at higher temperature.
In order to understand the mechanism underlying the lateral size dependence of friction, we investigated the evolution of the distribution of surface density of the contact layer for different box sizes. The surface density ρ ads n was computed by counting the number N ads of particles between the wall and the first minimum of the liquid density profile (separating the first and second adsorption layers) and dividing it by the interface area: ρ ads n = N ads /(l x l y ). Note that at each time N ads is an integer, so that the instantaneous surface density can take only discrete values. The probability distribution of the surface density is shown in Fig. 3 for different box sizes. While the average surface density remains constant, as expected from the similar vertical density profiles (Fig. 1), the probability distribution of surface density is qualitatively different between the largest box and boxes with sizes below ∼10 nm. In the largest box, the probability follows a Gaussian-like smooth bell curve centered around the average value; in contrast, in the small boxes the probability is strongly peaked at a value higher than the average surface density.
To connect the structural differences reported in Fig. 3 and the differences in friction coefficient presented in Fig. 2, we focused on the smallest system (l x = 3.13 nm), and we explored the correlation between the instantaneous friction coefficient and the instantaneous surface density. To that aim we extracted from the simulation the average values of friction stress and slip velocity for the different values of the instantaneous surface density, τ (ρ ads n ) and v s (ρ ads n ) , and we computed the friction coefficient at a given surface density from their ratio: λ(ρ ads n ) = τ (ρ ads n ) / v s (ρ ads n ) . The result is presented in Fig. 4. The liquid-solid friction is much lower at the largest surface densities, which are more probable for small systems, explaining why the friction coefficient is overall lower in small systems.  As a final step to explain the correlation between friction and surface density, we computed the two-dimensional Fourier transform of the structure of the contact layer at the two most probable surface densities (7.33 and 7.64 nm −2 ) for the smallest system: with k being a two-dimensional wave vector and x j (t ) being the two-dimensional position of atom j in the contact layer, among a total of N (t ) atoms. Figure 5(b) presents the time average of |c(k, t )|. However, averaging in time smooths out some information on the two-dimensional structure of the liquid in the contact layer because the structure has a rotational degree of freedom. Therefore, we also applied a rotation procedure to synchronize the instantaneous |c(k, t )| before averaging them so that the instantaneous Fourier transform at each step could have the same principal axes [ Fig. 5(c)]. Figure 5(c), highlighting the averaged intrinsic structure of the contact layer with rotational diffusion removed, shows that at both surface densities particles arrange along a hexagonal lattice, although the peaks in the Fourier transform are less pronounced for the low density, indicating that the structure has more defects. This is confirmed by visual examination of snapshots of the contact layer shown in Fig. 5(a). Generally, we found that higher densities corresponded to more ordered structures, with thinner peaks in the Fourier transform. In the hexagonal packing observed at large surface density, the firstneighbor distance between fluid molecules estimated from the wave vector of the peaks in the Fourier transform [ Fig. 5(c)] was ∼0.391 nm; this distance is commensurate with the box lateral size and suggests that the hexagonal structure is promoted by the finite size of the simulation cell.
The decrease in friction when the surface density goes from 7.33 to 7.64 nm −2 can then be explained by the dependence of friction on commensurability between the solid surface structure and the contact layer structure [45]. Since the intrinsic contact layer structure is hexagonal and the wall surface structure has a square symmetry, they are incommensurable, and disorder can only increase the commensurability. As a consequence, the high density contact layer with fewer defects generated a smaller friction than the low density contact layer. The key role of commensurability can also be seen in Fig. 5(b), which represents the directly averaged structure of the contact layer. It appears that the defects in the hexagonal structure of the contact layer at small surface densities let the contact layer adapt to some extent to the underlying wall structure, so that new peaks appear in |c(k, t )| with the square symmetry of the wall surface; indeed, the new peaks in the Fourier transform of the low surface density structure [ Fig. 5(b)] correspond to the lattice parameter of the underlying solid fcc crystal, i.e., √ 2r eq . In contrast, at large density the hexagonal symmetry of the contact layer structure is preserved. As a last hint pointing to the lower commensurability at high surface densities, we observed that the contact layer was slightly farther from the wall at high densities (see Appendix B), which also contributes to a decrease in friction since the amplitude of friction decreases very fast as a function of the distance to the wall [46,60,61].
Finally, we explored the robustness of the observed finite size effects by changing the temperature and the wetting properties. Indeed, we saw that slightly increasing the temperature (from 85 to 90 K) or decreasing the wetting coefficient (from 0.774 to 0.387) could kill the size dependence of friction (see Fig. 2). In Appendix C, we also show that reducing the wetting coefficient has a consistent impact on the probability distribution of surface densities, which goes from a peaked distribution to a smooth bell curve. Therefore, while finite size effects on friction can be very large and persist even for long boxes, it is reassuring to see that in most cases, where the liquid is far from its melting temperature or when the surface is not very wetting, such finite size effects should be negligible in practice.

IV. CONCLUSION
To measure liquid-solid friction in molecular dynamics simulations, a slab configuration is commonly used, where a liquid is confined between parallel walls. Periodic boundary conditions are used in the lateral directions to mimic an infinite slab, but the finite lateral size of the simulation box in the periodic directions could still impact, in principle, the measured friction coefficient.
Here we explored such finite size effects by simulating a generic Lennard-Jones liquid confined between parallel walls.
We focused on the case of a very wetting liquid close to its melting temperature. We observed a large effect of the box size along the flow direction on the measured friction coefficient, persisting up to a box size of ∼30 particle diameters.
We have then connected the change in friction to the surface density of the first adsorption layer. While the box size has no significant effect on the average surface density, we have shown that in small boxes the distribution of surface density in the first liquid adsorption layer became peaked at values larger than the average. We have also highlighted that, for these larger values of surface density, the first adsorption layer became more structured and less commensurable with the wall structure, thus explaining the decrease in friction.
Finally, we have shown that for less wetting cases or liquids further from their melting temperature, the strong finite size effects observed here disappeared and that boxes down to ∼9 particle diameters could be safely used. Nevertheless, we suggest that when exploring unknown liquid-solid interfaces, care should be taken with possible effects of the lateral box size on the measured friction. of a hemicylindrical droplet on a flat solid surface, where the potential model including the wetting coefficient α, the solid wall with fcc structure, its lattice parameter, and the exposing face of (001) were the same as in the present shear-flow system in Fig. 1. Eight layers of solid wall were located on the bottom (−z direction) of the rectangular simulation cell of l x × l y × l z = 39.2×3.92×16.0 nm 3 with periodic boundary conditions set in the x and y directions and a mirror boundary condition at the top boundary in the z direction. The temperature of wall particles in the second layer from the bottom was controlled by using the Langevin thermostat in all directions at the same control temperature of 85 K = 0.703ε LL /k B as the present shear-flow system. A hemicylindrically shaped droplet was formed on the solid surface with the droplet axis parallel to the y axis as an equilibrium state under liquid-vapor coexistence. The cylindrical liquid droplet was first equilibrated away from the solid surface, and after its automatic adsorption onto a solid surface followed by a further equilibration run of 20 ns, an initial equilibrium hemicylindrical droplet was obtained. The average of 40 ns thereafter was used for the analysis, where we obtained two-dimensional density distributions in the frame of reference relative to the center of mass of the droplet, considering that the droplet showed a random Brownian motion on the smooth and flat solid surface in this study.
The density distributions for three α values are shown in Fig. 6. As observed before [57,[62][63][64], multiple adsorption layers with a thickness about 1-2 nm were formed at the liquid-solid interface, and a hemicylindrical liquid-vapor in- terface with a uniform curvature was observed above. This indicated that the liquid-vapor interfacial tension was uniform except around the contact line. The contact angle θ of the droplet was defined by the angle between the least-squares fitting circle to a density contour line of ρ = 400 kg/m 3 and the fluid-solid interface, where for the former, we used the density contour only above the height where the density in the droplet was considered constant and the liquid-vapor interface did not overlap the adsorption layers. The details of the choice of the fluid-solid interface position based on a mechanical interpretation is described in Ref. [62]. As seen in Fig. 6, cos θ increased linearly to the increase of α. Figure 7 represents the vertical density profiles in the contact layer for the two most probable surface densities in the smallest system, with l x = 3.13 nm. For the largest surface density, the contact layer is slightly farther away from the wall, consistent with its lower commensurability with the surface. Figure 8 represents the probability distribution of surface density for the smallest system (l x = 3.13 nm) and for two different wetting coefficients α. For α = 0.774, the probability is asymmetrical and highly peaked, as discussed in the main text. In contrast, for α = 0.387, the probability assumes a symmetrical bell curve, reminiscent of the probability distributions observed at much larger box sizes for α = 0.774. Consequently, one can expect that finite size effects on friction will be negligible for α = 0.387, as observed in the main text.