Unexpected finite size effects in interfacial systems: Why bigger is not always better—Increase in uncertainty of surface tension with bulk phase width

Air-liquid interfaces are very sensitive and peculiar systems that show complex behavior, different from bulk materials, and continue to surprise both experimental and theoretical communities. In particular, the air-water interface is of fundamental importance to a huge range of environmental, biological, and industrial chemistry and possesses many anomalous properties.^1–5 

Many other previous investigations^6–13 have modelled the thermodynamic properties of water slabs in order to test a range of water forcefields of varying complexity. Typically these forcefields are designed to be used as explicit solvents for biochemical simulations of drugs, proteins, or lipid membranes. Consequently, their accuracy is generally assessed by being able to reproduce dynamic properties for single molecules and solvent cages in the bulk phase, as well as the average thermodynamic properties of single phases. However, it is also very important to be able to fully assess the physical properties that can be computed by a forcefield for any type of simulated system since they provide the ties that bind the microscopic and macroscopic worlds. Measuring the thermodynamic properties of an interface provides an understanding of not only how well atomistic simulations can reproduce experimental behavior but also how the time average of particle micro-states can represent a full ensemble average of a system that has an inhomogeneous density distribution. Although the bridge between statistical mechanics and thermodynamics using pairwise potentials has been well established in molecular mechanics simulations, methods to test distribution assumptions rely on studying systems of a size and accuracy that is just beginning to become computationally feasible.^12,14,15

In order to assess their suitability, forcefields for molecular mechanics simulations are generally parametrized to accurately recreate the homogeneous and isotropic mechanical, electronic, and thermodynamic properties of the materials and molecular species that they aim to emulate. The properties of inhomogeneous, anisotropic regions are assumed to arise from emergent behavior between atoms and molecules located nearby. In addition, it is necessary for simulation programs using periodic boundary conditions to employ a cutoff radius r_c for non-bonded forces between particles that is typically no greater than half the length of the shortest cell dimension.¹⁶ Consequently, extending r_c to provide a more accurate description of long-range interactions requires expansion of the simulation cell itself and therefore a significant increase in computational resources. Methods to avoid this demand include Ewald summation techniques that converge long-range contributions on the complex plane and have become a standard practice in molecular simulation codes for truncated electrostatic interactions. Recently these techniques have also been extended to dispersion forces modelled by a Lennard-Jones (LJ) type potential.^9,17–19 However, typically these missing pairwise interactions are not dealt with explicitly since they are considered to have a minimal effect on the system dynamics. The impact of these two assumptions on interfacial and surface regions of atomistic and molecular systems can be significant. Consequently, corrections are required for any property of interest and have traditionally been formulated using a continuum model during post-processing of data from a simulation on an ad hoc basis. For water models, these corrections can account for up to 25% of the absolute value of the computed properties for a typical LJ radial cutoff of 8 Å. It has been recently suggested that van der Waals forces are key to the formation of the hydrogen bonding network in water that gives rise to many of its unique and unusual properties.²⁰ Therefore, more suitable methods for including long-range LJ interactions in systems containing inhomogeneous density distributions are desirable.²¹ 

Considering the sensitivities of interface regions, it may appear beneficial to employ higher theory ab initio methods when studying them, especially for fluids with unique properties such as water. However, such techniques are still limited by the computational cost for large scale systems, which are necessary to produce converged surface behavior without artefacts.^22–27 Therefore for this investigation, we have chosen to focus on maximising the number of particles in our simulations at the expense of limiting model complexity and so employ a suitable non-polarisable water forcefield with fixed bond lengths and angles. In particular, we use the four-site transferable intermolecular potential, 2005 variant (TIP4P/2005)²⁸ forcefield, which has been shown to reproduce surface tension values across a range of temperatures between 273 and 600 K that are accurate to within 1–3 mJ m⁻² of experimental values.^7–13 In order to highlight the importance of dispersion forces and as a comparison between forcefields of varying complexity, we also choose to include simulation data from the surface thermodynamics of a LJ fluid. We then use parameters derived from experimental measurements of the cohesive forces in liquid argon²⁹ to compare with the experimental surface thermodynamic properties.

Despite all these considerations, simulated surface tension values of the air-water interface accurate to within 0.2 mJ m⁻² of experimental data remain from an investigation performed over 20 years ago, using systems containing only 1000 molecules of the fixed charge three-site extended simple point charge (SPC/E) forcefield.⁶ It is clear, therefore, that the intricacies of modeling air-liquid interfaces are still far from understood and any investigation that would push the boundaries of this field can be considered a vital contribution to the current knowledge.

We show here that the variance of simulated surface tension values, derived from the statistical mechanical description of interfacial stress by Kirkwood and Buff³⁰ (KB), is dependent on the bulk phase width-to-surface area ratio. The KB method is the standard approach to calculate pressures and surface tensions of each system configuration generated during a molecular simulation. Alarmingly, the linear dependence of this finite size effect on the width of the bulk phase appears to contradict the law of large numbers, from which we would expect bigger simulations to yield ensemble averaged properties with a greater numerical precision. Therefore the identification and explanation of this previously over-looked finite size effect are of significance for the design of simulations either employing a constant interfacial tension or investigating the behavior and reproducibility of surface tension values by particular forcefields.

We go on to demonstrate this significance by investigating the surface tension, surface energy, and surface entropy of a range of air-liquid simulations with differing liquid phase width-to-surface area ratios. We show that due to our reported finite size effect there is an advantage in statistical accuracy as well as an implicit theoretical advantage in using the configurational energy to estimate surface energetics, rather than calculating all three thermodynamic values from studying the KB surface tension alone. We finally attempt to provide some general guidelines for cell dimensions of interfacial simulations for which the variance of KB surface tension values is important.

There are three thermodynamic surface properties that can be abstracted from an MD simulation to give an insight into the behavior of molecules at an interface. These can also be compared against both experimental and simulated literature values. The surface tension γ of a system can be described thermodynamically in terms of energy U, temperature T, and entropy S as the change in Helmholtz free energy ϝ = U − TS, upon the change in surface area A, whilst maintaining a constant temperature and volume V,³¹ 

Likewise, the change in energy and entropy upon the change in the surface area can be defined in a similar way, yielding expressions for the surface energy U_s and surface entropy S_s,³¹ 

These surface properties are linked through the following equation and can be used to assess the ability of simulations to reproduce the behavior of particles at interfaces:

The most suitable forcefield, simulation procedure, and theoretical methodology should produce a balanced set of γ, U_s, and S_s values in comparison to experimental results. This is not guaranteed, as only 2 of these 3 variables are independent. There are multiple independent routes to estimate surface thermodynamic properties using molecular simulations, although they can be abstracted from the behavior of γ alone. Simulated surface tension values are commonly defined by Harasima,³² based on the KB derivation,³⁰ as the difference between the total pressure normal P_N and tangential P_T to the surface plane. For a periodic slab-like system with two interfaces in the xy plane and cell length L_z along the axis normal to the plane (as shown in Fig. 1), the surface tension across the whole system can be calculated using the following definition:

The global pressure of a molecular simulation can be defined at each time step in terms of the net force exerted by N particles upon the walls of the system. For a 3-dimensional many body system, this can be represented as a 3 × 3 Cauchy stress tensor P, where each component, P_αβ, represents the force in the β direction acting on the surface plane that has a normal vector in the α direction. The pressure components include the kinetic energy of the system and the virial sum of the distance r_ij and force F_ij vectors of each pairwise interaction between particles i and j [Eq. (5), where m_i is the mass of particle i, $viα$ is the α component of the particle’s Cartesian velocity v_i, and $rijα$ and $Fijα$ are the α components of the Cartesian distance and force vectors, respectively, between particles i and j],

Combining Eqs. (4) and (5) gives us a representation of the instantaneous surface tension for any system configuration in terms of pressure tensor components,

Values of S_s can then be inferred from γ and T using Eq. (1) and the relation S ≡ −(∂ϝ/∂T)_V,A,³¹ 

Subsequent estimates of U_s are provided via the rearrangement of (3). The expectation value of γ from a single simulation can be represented as an ensemble average $γ$ and calculated as the mean value in a time series γ_k of length τ from a system at equilibrium,

However, we have discovered that there is a significant dependence on the variance of surface tension values recorded during a single simulation, Var[γ], on the cell-dimensional ratio L_l/A, which can have a significant impact on the assessment of a forcefield’s ability to reproduce interfacial behavior via Eqs. (6) and (7). We provide the following theoretical explanation for this unexpected finite size dependency.

Our aim is to formulate the relationship between Var[γ] to L_l/A. Each element γ_k in Eq. (8) can be represented by combining Eqs. (4) and (5) for a system containing N pairwise interacting particles with 2 interfaces in the xy plane so that A = 2L_xL_y. It is assumed that at equilibrium, the kinetic terms of the pressure tensor components are equivalent to ρk_BT and therefore cancel out in Eq. (6), leaving only the virial sum. The pairwise force is also assumed to be proportional to the radial distance between particles i and $jrij=rijx2+rijy2+rijz2$⁠. Consequently, the expression for γ_k becomes

Without knowing the expected value of the summation, the nature of the forces between particles, or the distribution of radial distances, we separate out the virial contributions χ_ij from each pairwise interaction

and split the summation over all pairwise interactions into contributions where both particles i and j are in the bulk region and where at least one is in the surface. By doing so, we group surface-surface and bulk-surface interactions together, as it will be shown later on that only the bulk-bulk interactions are the source of the finite size effect under investigation. We expect the bulk summation to contain n_B elements and the surface summation to contain n_S elements,

Subsequently, combining Eqs. (8) and (11) yields

The expectation value for the contribution from the bulk phase $χB=0$⁠. Therefore we can deduce that increasing the number of particles in the system by increasing the width of the bulk liquid phase L_l without a subsequent increase in A will not affect $γ=χS/A$⁠. However, we show below that estimating Var[γ] in an equivalent way leads to a non-zero contribution from the bulk region. We have assumed no covariance between each element χ_ij so that the variance contributions from each region can be treated as a linear summation,

Combining Eqs. (8) and (13) shows that the ensemble average of $χB2$ from a single molecular dynamics simulation is normalised over time by the number of system configurations τ, but not over the number of contributing bulk pairwise interactions n_B. Therefore $χB2$ grows with n_B, and since $χS2$ and $γ2$ are independent of the bulk phase, Var[γ] will also grow in proportion to $χB2$⁠,

Consequently, increasing n_B by expanding the width of the liquid phase L_l will lead to a subsequent increase in Var[γ]. We believe this to be the cause of the observed relationship between Var[γ] and L_l/A. We shall now attempt to derive an analytical expression from Eq. (13) that shows the linear dependence of Var[γ] on L_l/A.

We first represent $χB2$ in terms of a pairwise average Ω_B [Eq. (11)],

We expect Ω_B to be an intensive property for a system in equilibrium,

This allows us to demonstrate using Eq. (14) that the variance in surface tension values Var[γ] is linearly proportional to n_B,

The number of pairwise interactions n_B present in a simulation typically depends on the volume, liquid phase density ρ_l, and whether a cutoff radius r_c has been used to truncate force interactions. Upholding the minimum image convention when calculating the virial sum of a periodic system¹⁶ is typically achieved by employing a spherical cutoff radius of $rc≤minLx,Ly,Lz/2$⁠. Therefore we now adopt a cylindrical coordinate system r(z, r, θ) and define a simple step function representing the liquid phase particle density along the z axis normal to the interface plane, with respect to the centre of mass z = 0,

where Θ(z) is the Heaviside function,

An estimation of the total number of pairwise interactions present in a system described by Eq. (18) can be given by using a spherical shell continuum model,

Rearranging and performing integration by substitution where s = cos θ brings the density product ρ(z)ρ(z + rs) into a single integral over z,

Considering the form of ρ(z), we geometrically estimate the integral of ρ(z)ρ(z + rs) over all space as the product of the square of the bulk density $ρl2$ and the overlapping distance between the two rectangular density functions L_l − |rs| (Fig. 2),

The absolute value of rs is used to account for the overlap’s conservation of parity. Incorporating this into Eq. (21) yields

which can be solved analytically to give n_B,

Subsequently, combining Eq. (24) with (17), we form an expression for the linear dependence of Var[γ] on L_l/A,

Therefore, we expect the variance of calculations of γ via the KB methodology for pairwise interactions described by Eqs. (5) and (6) to be linearly proportional to the system dimensional ratio L_l/A. This expression assumes the use of a radial cutoff for the direct calculation of pairwise force components and does not take into account the effect of techniques to deal with long-range interactions and periodicity, such as the Ewald summation or Particle-Mesh Ewald (PME) summation. It is therefore most applicable to molecular systems dominated by dispersion forces, for which the LJ potential has traditionally been applied via direct summation. We now present an example of how this finite size dependency can have a significant effect in the calculation of simulated surface thermodynamic values.

Using Eq. (7), we can estimate U_s and S_s from the behavior of γ with respect to T. We therefore perform a series of interfacial simulations across a range of temperatures for TIP4P/2005 water and LJ argon to investigate this behavior. We then use these data to parameterize empirical equations describing the relation of γ with T, sourced from the work of Vargaftik, Volkov, and Voljak³³ for water and from the work of Lemmon and Penoncello³⁴ for argon. The same procedure is performed on various experimental and simulated data sets reported in the literature.

Vargaftik, Volkov, and Voljak³³ suggested the following equation as an amendment to Guggenheim³⁵ to fit reference surface tension data for water, where T_c is the critical temperature at which γ = 0:

The derivative of (26) with respect to T can then be used in order to produce an estimate of S_s,

Consequently, U_s values can then also be derived from the rearrangement of (3),

Lemmon and Penoncello³⁴ suggested the following equation to fit reference surface tension data for argon to (29), where T_c is the critical temperature at which γ = 0:

The derivative of (29) can be used in order to produce an estimate of S_s,

Consequently, U_s values can then also be derived from the rearrangement of (3),

We propose a methodology to calculate U_s values independent of γ by using (2) to compare the internal energies of a range of simulated systems that possess different surface area to volume ratios (i.e., different slab “thicknesses”). We scale the energies and surface area of our system by the number of molecules N and assume a negligible expansion of the slab during the simulation. Therefore the central phase of each slab replicates a bulk environment of homogeneous density. The total energy of an interfacial system U is given as a sum of the bulk liquid U₀ and surface phases U_sA, leading to the following expression:

In this way, we define the surface energy by the following equation, which is gained by linear regression (LR) of the variables U/N and A/N:

To calculate U_s by this method, we propose to run simulations of liquid-vapour systems across a range of A/N dimensional proportions and perform a linear regression with the time-average total configurational energy $U/N$ for each system. An estimation of U_s can be used along with the average γ across all our systems in a rearrangement of (3) to estimate the surface entropy S_s. It is assumed that simulated interfacial properties are free from finite size effects above a certain slab dimension, typically L_x = L_y ≳ 10σ across the surface plane^23,24,36,37 and L_l ≳ 15σ along the normal axis²⁷ for particles using a LJ radius of σ. In addition, the use of trajectory corrections to obtain converged liquid densities further minimises the influence of system size artefacts. Therefore we do not expect $γ$ to be dependent on A/N, although considering that (A/N)⁻¹ ∝ L_l/A, we expect Var[γ] to differ between simulation dimensions. The average $γ$ across all system dimensions will therefore need to be calculated as an error weighted mean. It must be noted that there are more elegant solutions available in the literature to calculate γ, U_s, and S_s, such as the Test-Area (TA)¹⁴ and Two-Phase Thermodynamic (2PT)^38,39 methods, which have been used to estimate all three properties from a single simulation^12,40 without susceptibility to the artefact we present. We provide our linear-regression calculation to highlight an example of where the finite size dependence of Var[γ] and L_l/A in the KB methodology becomes significant. It is possible that alternative routes to calculate surface thermodynamic properties will be superior, especially when investigating systems that require large bulk regions. However, calculations using the virial remain the standard estimation of pressure related properties for molecular simulations.

It has been ascertained that for LJ fluids, such as argon, an accurate representation of long-range dispersion forces in an interfacial simulation requires cutoff radii on the order of r_c ≥ 5σ ≈ 17 Å.^41,42 The same investigations demonstrated that Janeček’s long-range corrections²¹ can effectively replace these lost interactions whilst running at significantly lower cutoff radii of r_c ≥ 3σ ≈ 10 Å. Therefore we implement the energy, force, and surface tension corrections for systems with inhomogeneous densities, as developed by Janeček.²¹ The methodology employs a planar surface geometry rather than a spherical shell geometry, allowing radial and angular potential terms to be solved analytically. The resultant global integral that ranges over the surface normal axis can be estimated by a summation using a discrete density histogram ρ(z), leading to a significant increase in the speed of calculation. Consequently this allows force corrections to be applied to particle trajectories during the simulation. Energy and surface tension corrections are then applied during post-processing using a time average discrete density distribution $ρ(z)$⁠.

Ensemble averages of internal energies $U$ and surface tension $γ$ from a single simulation are calculated from discrete time series that have a tendency to be highly correlated in time. The variance of any measurement A taken from a time series of length τ can be shown to have an inverse dependence on the number of uncorrelated measurements $τAC/τ$⁠, where $τAC$ is the correlation time or minimum time length between independent samples. If measurements are uncorrelated, then $τAC=1$ and the variance is directly proportional to τ⁻¹,

Therefore, the standard error in the mean of both properties needs to be amended by $τAC$⁠. We estimate this parameter by studying the behavior of block averages. A time series containing N measurements is split into “block” lengths of time τ_B, for which there is a mean value Ā_B. The variance of n = N/τ_B blocks is given by

so that according to Eq. (34),

Therefore if we can find the value of τ_B at which $f(τB)=τBVarĀB/Var[A]$ converges, we can compute the true standard error in the mean σ_A,

In practice, this is easier to estimate if we assume that there exists a linear relationship between $τB−1$ and $f(τB)−1$ so that the linear regression of both would yield an intercept of $τAC$⁠.

We outline a generic procedure to form an equilibrated liquid-vapour system with periodic boundary conditions, as illustrated in Fig. 1. Systems containing electrostatic interactions employ the AMBER 12 PMEMD code.⁴³ This includes a B-spline Smoothed PME (SPME)⁴⁴ integration routine with a full 3D Fast Fourier Transform (3DFFT) charge density mesh spacing of 0.5 Å and a time step of Δt = 2 fs. A radial cutoff of r_c = 10 Å for the direct summation of non-bonded forces is employed. Molecular geometries with covalent interactions use bond constraints via the SHAKE algorithm. Steps involving a constant temperature used a Langevin thermostat with a collision frequency of 10.0 ps⁻¹, whereas steps involving constant pressure used a Berendsen barostat. Our simulation procedure is then as follows:

1. Fill a cuboid of dimensions L_x = L_y, L_l with a homogeneous mixture of desired molecules using leap, AMBER’s general utilities software.

2. Heat the cuboid system from 0 to T K over a period of 0.1 ns.

3. Equilibrate the liquid system under NPT conditions for 0.5 ns at T K with P = 1 atm.

4. Add a vacuum region along the z axis, adhering to the dimensions in Fig. 1, where L_v = 2L_l ≥ 100 Å.

5. Equilibrate the liquid-vapour system under NVT at T K conditions for 0.5 ns with Janeček’s trajectory corrections.

6. A further simulation of the cell under NVT conditions at T K lasting for 8 ns is performed with Janeček’s trajectory corrections. System configurations are recorded every 1000 Δt.

Systems of TIP4P/2005 water and LJ argon (forcefield parameters shown in Table I) using various combinations of dimensions 15σ ≤ L_x, L_y, L_l, L_v ≤ 35σ were simulated across a range of temperatures between 273–478 K and 85–135 K for water and argon, respectively, in order to predict surface thermodynamic properties (see the supplementary material). The length of the vapour region L_v = 2L_l ≥ 100 Å was chosen in an attempt to ensure a balance of computational efficiency and minimal influence of electrostatic interactions between periodic slab surfaces. Estimations of the $ρ(z)$ were calculated from 4000 system configurations, recorded over the simulation described in step 6. Energy and surface tension long-range corrections were calculated for both models, as described by Janeček,²¹ whilst a block averaging method (as described in Sec. II G) was used to estimate both the errors in the mean internal energy σ_U and surface tension σ_γ. Anisotropic variations in cell dimensions were investigated for LJ argon, following Fig. 1.

We first present data illustrating the finite size dependence of Var[γ] on L_l/A for surface tensions calculated via the KB methodology.

An initial insight into the behavior of surface tension with system dimensions is provided in Figs. 3 and 4, where we have compared the ensemble average $γ$ for each cell dimensional ratio L_l/A. We see that there is a slight oscillation and a general increase in the uncertainty of $γ$ with L_l/A, both of which are most apparent for argon (Fig. 4). An oscillation in $γ$ has been previously reported for both KB^24,36,37 and TA¹⁵ estimates, though it has been assumed that this finite size effect was predominately proportional to A alone and was a result of the surface area being insufficiently larger than the molecular hard sphere radial area σ.² Our relatively small oscillations, lying within the range of one standard error, are therefore to be expected since, as mentioned previously, our simulations have been designed to avoid the dimensional range in the xy plane, in which area finite size effects are significant. From our LJ argon data in Fig. 4, we also report no additional impact of cell anisotropy on $γ$⁠, which can be problematic for systems with too small dimensions.²³ 

From Figs. 5 and 6, we can clearly see the linear relationship predicted by Eq. (25) for both water and argon simulations. We find that this relationship holds for all systems of varying combinations of L_l and A investigated. Additionally, it appears that the increase in surface tension variance is significant for both liquids, suggesting that systems dominated by dispersion interactions only using direct force summations within a radial cutoff r_c, and those including electrostatic interactions employing Ewald techniques for long-range interactions outside this range, are both susceptible to this finite size effect. However, further testing of simulations using the Ewald treatment of dispersion LJ forces is required to confirm such an observation. Nevertheless, we expect the Var[γ] of other simulated liquids to be also sensitive to L_l/A, as long as the surface tension has been calculated via the KB method for pairwise interactions, and the simulation has employed a radial cutoff for the direct sum of non-bonded forces.

We shall now go on to explore estimations of γ, U_s, and S_s, as predicted by our two methodologies described in Secs. II D and II E to examine possible consequences of this finite size effect when calculating thermodynamic properties.

Tables II and III include average surface thermodynamic values γ, U_s, and S_s of water and argon, calculated via a variety of methods including this investigation’s linear regression approach. KB values are calculated from the parametrisation of Eq. (26) for water and Eq. (29) for argon, as described in Sec. II D using surface tension data across a range of temperatures. All thermodynamic properties are then estimated via (26), (28), and (27) for water and (29), (31), and (30) for argon. KB/LR (Kirkwood and Buff/Linear Regression) values are estimated using a series of simulations with differing A/N ratios at 298 K for water and 85 K for argon, as described in Sec. II E. The surface tension γ is calculated as an error-weighted mean of all simulation averages presented in Figs. 3 and 4, whereas U_s is calculated via linear regression of $U/N$ and A/N, and S_s inferred from Eq. (3). Full results for the parametrisation of Eqs. (26) and (29), as well as the linear regression calculation of surface thermodynamic properties, can be found in the supplementary material.

For TIP4P/2005 water, we report KB/LR γ = 67.9 ± 0.6 mJ m⁻² in comparison with KB γ = 67.4 ± 1.8 mJ m⁻² at 298 K. We notice a higher level of precision for our KB/LR estimate, which we attribute to lower variations in γ between system dimensions than across temperature ranges, illustrated by the variations in Figs. 3 and 7. TA and 2PT calculations from the literature also produce estimations of γ with a greater precision than KB values, which could be explained by the expected insensitivity of these two methods to L_l/A. Despite possessing a relatively high level of precision in general, our simulated $γ$ estimates at values of T < 400 K are consistently lower than those reported from other investigations by up to 2 mJ m⁻². Consequently, we generate KB estimations of thermodynamic properties using parameterized equations (26)–(28) that are ≈5% lower than experimental values at 298 K. Our predicted critical temperature T_c = 645.0 ± 3.1 K of TIP4P/2005 is also lower than the experimental value of water T_c = 646.9 ± <0.1 K. However, it is not clear that this should imply a systematic inaccuracy in our methodology since we have sampled over a much longer equilibration time and use a significantly finer PME mesh grid size (0.5 Å) than any previous study. In addition, there is a large ≈10 K variation between estimated T_c values from reported TIP4P/2005 surface tension data in the literature (see the supplementary material).

For argon, we report KB/LR γ = 16.2 ± 0.2 mJ m⁻² in comparison with KB γ = 16.1 ± 1.0 mJ m⁻² at 85 K. Again, we notice similar disparities in variation between both KB and KB/LR estimates, illustrated in Figs. 4 and 8. From Fig. 8, we also see that our simulated γ of LJ argon across nearly all temperatures is consistently above that of experimental and simulated values in the literature. Janeček²¹ commented on the tendency for trajectory corrected simulations to overestimate values of γ that would be expected for r_c = ∞, which may explain the difference between our predictions and other studies. Our predicted critical temperature T_c = 155.9 ± 1.9 K of LJ argon is similar to that estimated from the surface tension data of Gloor 2005 T_c = 156.0 ± 0.6,⁵¹ both of which are significantly higher than the experimental value of argon T_c = 150.7 ± 0.01 K. It has been ascertained that explicitly including three-body terms into the KB definition of γ for argon leads to a reduction in predicted surface tension values resulting in estimates that are much more comparable to experimental data.⁴¹ Therefore, it is possible that the overestimations simply derive from the limitation of the pairwise LJ potential to reproduce many-body interactions between noble gases.^52,53

The size of the error bars of $γ$ in Figs. 3 and 4 also gives an insight into the significance of Var[γ] and the requirement to run for long simulation times in order to achieve a converged ensemble average. It has been suggested^6,54–56 that the way in which long-range electrostatic forces are dealt in the Ewald summation has a significant effect on the accuracy of the virial sum. However, we use a PME summation to calculate periodic long-range electrostatics and have found no correlation between $γ$ or Var[γ] beyond a mesh displacement of 1 Å, in agreement with previous investigations.⁵⁶ It is worth noting that an equally weighted mean of $γ$ for each system L_l/A produces estimates of γ for TIP4P/2005 water and LJ argon that are equivalent to the error-weighted calculation to within one standard deviation. However, this is achieved by using a large number of different simulation cell dimensions than would normally occur. It will be shown next that the impact on the estimation of U_s and S_s via the inference of (∂γ/∂T)_V can become more significant.

Figures 9 and 10 show the linear relationship between $U/N$ against A/N for TIP4P/2005 water and LJ argon, respectively, as predicted by Eq. (32). The R² values of each fit are above 0.999, indicating that our KB/LR methodology produces estimates of the surface energy with a low error in fitting precision.

For TIP4P/2005 water, we report KB/LR U_s = 109.7 ± 0.3 mJ m⁻² in comparison with KB U_s = 108.3 ± 4.0 mJ m⁻² at 298 K. The difference between both estimates (+1.4 mJ m⁻²) is significantly greater than that of γ (+0.5 mJ m⁻²) and there is also a clear decrease in the uncertainty of KB/LR U_s in comparison to KB U_s. Despite being 10 mJ m⁻² lower than the experimental value of U_s = 118.0 ± 0.1 mJ m⁻², we see from Table II that both these estimates also lie within the range of simulated values available in the literature. Interestingly, the lower KB estimate appears most similar to simulated values using larger cutoff radii of r_c ≥ 12 Å (Mountain 2009, Pascal 2014), whereas the KB/LR value is more comparable to studies using r_c ≈ 9.5. The 2PT U_s estimate taken from the literature appears to lie within the range of reported KB values and also possesses an equivalent precision. However, the TA estimate taken from the literature is of significantly lower precision and accuracy in comparison to experimental data. Considering this was abstracted from S_s in the original paper,¹² inaccuracies in the TA U_s term were attributed to the statistical difficulties of obtaining reliable entropic values from MD simulations.^12,57

For argon, we report KB/LR U_s = 38.2 ± 0.2 mJ m⁻² in comparison with KB U_s = 40.1 ± 1.3 mJ m⁻² at 85 K. Again the difference between these measurements (+1.9 mJ m⁻², approximately 3/2 the standard error in KB U_s) is greater than that between γ (+0.1 mJ m⁻²). We also note that, similar to the surface tension, both estimates of U_s for argon give values that overestimate the experimental surface energy U_s = 34.8 ± 0.2 mJ m⁻², although they lie within the range of simulated predictions in the literature. We have already commented on the tendency for Janeček’s force corrections²¹ and the absence of many body terms in the virial sum to overestimate the surface free energy (γ) in LJ fluids such as argon, so it is not wholly surprising that the same appears to occur with U_s. Nevertheless, our KB/LR estimate is consistently closer to experimental U_s than KB, which may suggest that calculating surface energetic terms independently of γ can alleviate some of these issues.

Calculating the surface entropy via Eq. (7) provides an example of an estimation that is sensitive to Var[γ]. Considering the increase in the precision of our KB/LR U_s estimates, we therefore also expect our KB/LR S_s calculations to show greater precision than KB estimations. However, it is not clear whether they will also be more comparable to experimental data since this will depend on the accuracy of both γ and U_s, whereas KB estimates rely on the behavior of γ alone.

We report TIP4P/2005 KB/LR S_s = 0.140 ± 0.003 mJ m⁻² K⁻¹ in comparison with KB S_s = 0.137 ± 0.008 mJ m⁻² at 298 K. Again, despite a significant difference of around −0.015 mJ m⁻² K⁻¹ between our estimates and the experimental measurement S_s = 0.154 ± 0.001, both values lie within the range found in the simulation literature using the TIP4P/2005 forcefield. As expected, the precision of our KB/LR estimate is also significantly greater than that of our KB estimate. Similar to our surface energy results, the lower KB S_s estimate is more similar to literature values with a larger r_c, which is to be expected, since stronger cohesive forces would reduce the impact of greater thermal motion on the surface structure and therefore the magnitude of (∂γ/∂T). Similar to the surface energy, the 2PT S_s estimate taken from the literature appears to be of equivalent precision to our KB/LR estimate and also lies within the range of KB literature values. However, as previously discussed, the TA S_s estimate taken from the literature is of significantly lower precision and accuracy in comparison to experimental data due to the considerably larger variation of simulated entropic terms.^12,57

For argon, we report KB/LR S_s = 0.259 ± 0.005 mJ m⁻² K⁻¹ in comparison with KB S_s = 0.280 ± 0.010 mJ m⁻² K⁻¹ at 85 K. We notice a large discrepancy here between both methodologies of 0.021 mJ m⁻² K⁻¹. We believe this to be caused by the influence of Janeček’s²¹ trajectory corrections (as previously mentioned) and also the behavior of our simulation as it approaches T_c ≈ 150 K, at which insufficient cohesive forces can cause the system to prematurely form a single phase, causing γ → 0. Inaccuracies at this region can significantly affect the functional fit of Eq. (29) and therefore KB predictions of γ, U_s, and S_s. We note from Table III that applying the same technique to simulation data from other investigation can have a similar outcome. As a consequence, we have limited our temperature range to ≤135 K. Due to these difficulties, in combination with our KB/LR S_s estimate being significantly more accurate to the experimental measurement S_s = 0.255 ± 0.001 mJ m⁻² K⁻¹, it appears beneficial to estimate surface thermodynamic values of LJ argon by measuring γ across a range of cell dimensions at a single temperature, rather than across a range of temperatures at a single set of cell dimensions.

Generally, there is consistency in all three surface properties produced via our KB/LR methodology for the TIP4P/2005 water and LJ argon models, which are also highly comparable to experimental and other simulated estimates. However, due to a newly identified finite size effect affecting Var[γ] and difficulties encountered whilst measuring simulated surface tension values at high temperatures, it appears that estimating surface energies via linear regression of $U/N$ with A/N provides a more robust method of calculating surface thermodynamic properties than studying γ alone.

Including an explicitly determined value of U_s in Eq. (3) reduces the uncertainty of S_s since there is a much lower variance in configurational energies Var[U] throughout the molecular dynamics simulation compared to that of the differences in components of the pressure tensor. We also expect Var[U] to be independent of L_l/A, unlike Var[γ], removing the dimensional restriction on larger simulations that can decrease the uncertainty in U_s further. As described, the dependence of Var[γ] on L_l/A causes expansion of the slab thickness L_l to lead to a greater uncertainty in calculated surface tension values, whereas expansion of the surface cross section A reduces the uncertainty. This result has strong implications for investigations either assessing the accuracy of surface tension estimates for a new forcefield or using a constrained interfacial tension, commonly employed in lipid-bilayer simulations.⁵⁸ In both cases, it is expected that using systems with a large dimensional ratio of L_l/A will require longer simulation times to either yield converged estimates of $γ$ or lead to a constrained interface at equilibrium.

A similar dependence of γ on the system size has been observed by Schmitz et al.,^59,60 who theorised that the typical length over which the interface position fluctuates is proportional to the dimensional fraction $Lz/A$⁠. This “domain breathing” effect is also described as being independent of the system Hamiltonian, instead being a phenomenological theory shown to exist for both Ising models and LJ fluids, caused by density fluctuations in the bulk. However, we have not recorded any evidence at the moment to suggest a link between the variance in surface tension values derived by the virial summation and any extra mobility in the interface region, as predicted by domain breathing. It is also likely that these two effects are independent since in each of our systems the density profile and time-averaged thickness of the interfaces remained constant, indicating that the behavior of the interface region was relatively unaffected by the change in $Lz/A$ compared to the size of the variations in γ.

Considering that the increase in Var[γ] originates from contributions to the virial sum from particles in the bulk liquid, if the positions of any interfaces are known prior to the simulation, simply excluding any pairwise interactions outside this region could remove this size dependency. Typically these positions are not known unless a constraint is being placed on the system or the density distribution along the axis normal to the interface is being recorded. Simulations employing the Janeček method for long-range corrections to particle trajectories would possess such a density distribution, and so an amendment to the pressure tensor could be made for these cases. However, care must be taken when analyzing instantaneous particle distributions since the same density fluctuations can cause the exact position of interfaces to be unclear. In addition, methodologies to map the structure of the intrinsic surface^61–63 at each time step would be too computationally intensive to perform during the simulation and therefore reserved for post-processing analysis only. Unfortunately we are unable to directly minimise Var[γ] with respect to L_l/A in Eq. (25) since we have not resolved all possible terms contributing to Var[γ] and simply reducing L_l is likely to introduce additional finite size effects.²⁷ The simplest solution may simply be to maintain a low-dimensional L_l/A ratio beyond the range of any additional system size artefacts when designing simulations of interfaces (L_x = L_y ≳ 10σ, L_l ≈ 15σ, therefore L_l/A ≲ 0.15σ⁻¹), though it is important to remember that other methodologies used to calculate interfacial tension values, such as the TA and 2PT methods, are not expected to be susceptible to the same size dependency. Alternatively, calculating the surface tension via a local pressure profile may allow for contributions far away from the average position of the interface to be discarded. This solution was recently proposed by Yoon et al.⁶⁴ whilst performing post-processing analysis on nano-films of water and graphene. Finally, we have not assessed the impact of the ratio L_l/A on contributions to the surface tension that lie outside the direct summation of the virial. Therefore Ewald methodologies that include contributions on the imaginary plane may not show the same linear dependency on the width of the bulk phase. It is also expected that this may be the case for many body contributions since the expression for the pressure tensor components given in Eq. (5) only considers pairwise interactions.

We have proven the existence of a previously unexplained finite size effect on the variance of the Kirkwood and Buff description of γ for pairwise interactions that is linear to L_l/A. This is an unexpected and highly intriguing result since it infers that increasing the size of any simulated interfacial system may not always lead to ensemble average properties that possess a lower uncertainty. We derive an explanation of the origin of this dependency, as well as suggesting ways to combat or avoid it. In order to demonstrate possible implications of this finite size effect, we have developed a novel way to estimate the surface energy U_s of air-liquid interfaces using simulations of slabs with varying dimensional parameters. Using this method, we have calculated three thermodynamical properties γ, U_s, and S_s of the air-water and air-argon interfaces using the TIP4P/2005 forcefield and a parameterized LJ fluid. These values compare well to other experimental and simulated estimates from studies found in the literature and show lower susceptibility to the influence of L_l/A in comparison with calculations of all three via the thermodynamic inference of γ alone.

See supplementary material for data tables of simulated surface thermodynamic values generated by this investigation and parameters used in equations (26)–(31) for all studies presented in Tables II and III.

This work was supported by the EPSRC DTC Grant No. EP/G03690X/1. Simulations were performed using the IRIDIS High Performance Computing Facility at the University of Southampton, UK. We cite the developer group of AMBER 12 in Ref. 66.

We acknowledge the use of the MDTraj python library for analysis of simulation data.⁶⁵