Reducing uncertainty in simulation estimates of the surface tension through a two-scale finite-size analysis: thicker is better

Recent simulation studies of the surface tension γ, and other properties of thin free-standing films, have revealed unexpected finite size effects in which the variance of the properties vary monotonically with the in-plane width of the films, complicating the extrapolation of estimates of film properties to the thermodynamic limit. We carried out molecular dynamics simulations to determine the origin of this phenomenon, and to address the practical problem of developing a more reliable methodology for estimating γ in the thermodynamic limit. We find that there are two distinct finite size effects that must be addressed in a finite size analysis of γ in thin films. The first finite size scale is the in-plane width of the films and the second scale is the simulation cell size in the transverse direction. Increasing the first scale enhances fluctuations in γ, measured by the standard deviation of their distribution, while increasing the second reduces γ fluctuations due to a corresponding increased ‘freedom’ of the film to fluctuate out of plane. We find that using progressively large simulation cells in the transverse direction, while keeping the film width fixed to an extent in which the full bulk liquid zone is developed, allows us to obtain a smooth extrapolation to the thermodynamic limit, enabling a reduction of the γ uncertainty to a magnitude on the order of 1% for systems having a reasonably large size, i.e., O (1 μm).


Introduction
components of the pressure in the interfacial region. Even so, the pressure profiles of Harasima have been used in several surface-tension studies because it does not matter whether the interaction forces are distributed uniformly in several slabs, or in only two; when it is integrated to obtain the profile of the surface tension, the same result is obtained. 10, [20][21][22][23] As the contributions to the Harasima pressure profiles are located in the regions that originated the interactions, the density inhomogeneities at the interfaces are highlighted in the Harasima profiles. The Harasima normal pressure profiles are mechanically stable because the outermost regions show large negative pressures (attractive) holding together the whole liquid layer. Before the attractive zones, there are small peaks with positive pressure (repulsive), which, in comparison to the attractive ones, are smaller, and are not long and wide enough to burst the liquid layer. Here, the reduced surface tension * was calculated through its mechanical definition using the Harasima pressure profiles. 4,10 In this way, the whole system (two interfaces) is used to calculate de surface tension.
In order to obtain more realistic profiles, we allowed the system to move in the inhomogeneous direction, and calculated the density and pressures profiles every 100 steps. At the end, we averaged these by correcting the positions of the profiles according to their positions in the inhomogeneous direction of the center of the layer, which was calculated as the midpoint between the two positions of the Gibbs' dividing surface of each interface, K * , with the thickness of the layer, L * , corresponding to the separation between the positions of the two dividing surfaces.
This allowed us to obtain more realistic profiles of the calculated properties. To obtain the positions of K * , each of the density profiles, calculated every 100 steps, was adjusted to the hyperbolic tangent expression commonly used in vapor/liquid phase equilibrium studies: 11,19 * ( ) =

Results
Molecular dynamics are commonly carried out in finite simulation cells, and it is important to consider those finite effects of the simulation cell on the property calculations. 24-26 Two finitesize variables that can be studied are the length of the layer and the interfacial area. Lengths of the layer below a critical thickness induce a layer break-up, while at the critical thickness, they produce metastable systems with no real bulk liquid phase. 4 Layers that are wide enough correspond to coexisting phases that are thermodynamically stable, and it is expected that further increments in the thickness will not change the interfacial and coexisting average properties. The size effects in the interfacial area are not intuitive due to the periodic boundary conditions employed, which try to reproduce the dynamics of systems with infinite and continuous interfacial areas; however, the periodic images create artificial cohesive forces at the interfaces, affecting certain thermophysical and interfacial properties, while the surface tension remains almost insensitive. 4 The cohesive forces are stronger in simulation cells having smaller interfacial areas, and result of the localization of the Lennard-Jones atoms, and thus enhancing the interatomic cohesive interactions, in the interfacial region. Due to the cutoff of the Lennard-Jones potential used, we performed the simulations using a reduced interfacial length, Y * , from 16 to 150.4, and thicknesses from the critical length to a reduced length of 16.8. Also, for the system with the smallest interfacial area used, we simulated systems with reduced lengths of up to 67.8.
The artificial cohesiveness forces induce a more structured behavior of the systems in the bulk phases and at the interfaces. Figure 2 shows the density profiles averaged over 100 timesteps for systems with L * = 16.8 (far away from the critical length 4 ), T* = 0.72 and Y * = 16 and 150.4, which represent the smallest and largest surface areas studied. Even these profiles can be more well defined using more timesteps to average the profiles, the average over 100 timesteps allow us to understand the instantaneous density profile estimates. The system with a Y * of 16 was more In a similar way, the average g is also insensitive to the value of Y * employed in the simulation, as previously reported, 4 but also it is insensitive to the length of the layer, as shown in average values of the two sets at different Y * , for each L * , we found very small differences that were less than the size of the symbols. On the other hand, we found a similar behavior between the distributions of the bulk liquid density and the distributions of the surface tension, in which the distributions of narrowed as Y * increased, which manifested as lower standard deviations.
We studied the origins of the insensitivity of the behavior of the average value of * with the size of Y * at constant L * through analyzing the components of the average pressure profiles. In In laboratory units: where K is a characteristic size of the simulation cell in the transverse direction. At very short Y sizes (close to K ), the interfacial thickness is equal to their intrinsic value. A study on water simulations have employed Equation 3 to predict the surface tension at a specific thickness of the layer, 27 but no comparisons have been made at different thicknesses of the layer, and none have explored the limits of Equation 3 at very narrow thicknesses. Figure 9 shows (∆ * ) 7 as a function of ( Y * ) for the systems with L * of 8.0 and 16.8. Linear regressions of the two sets of data produce almost parallel lines, which also reflect the insensitive behavior of * not only with Y * , but also with L * .
The standard deviation of * , k * , is plotted as a function of L * in Figure 10 for the simulated systems using sizes of Y * of 16, 40, 120 and 150.4. We calculated the profiles of k * for each Y * from the first stable layers at its critical length, 4 at long Y * values. Figure 12 shows k * as a function of Y * for two thin layers with L * of 8.0 and 16.8. The two sets of data show an asymptotic behavior of k * with Y * , in agreement with the expression developed by Longford et al. 1 , which can be modeled as an inversely asymptotic function of the form: As wider lengths of the liquid layer are employed, the interfaces become freer, and there is no coordination between these, making them free to develop their own dynamics.
The size effect due to the finite size of the simulation cell in the interfacial area is clearly an artefact of the methodology employed to simulate phase equilibria under periodic conditions in the interfacial surface, but for systems confined to these interfacial-length scales, probably this size effect has a physical meaning and future work should elucidate its influence. The size effect due to the thickness of the layer is only important for narrow layers in the scale of the nanometer, beyond that its influence decreases until they disappear at larger thicknesses layers. In the future, we also plan to assess the effects of additives and multicomponent equilibria to study the dependence of finite size effects on the distributions of their interfacial properties.

Conflicts of interest
There