We present an unexpected finite size effect affecting interfacial molecular simulations that is proportional to the width-to-surface-area ratio of the bulk phase *L*_{l}/*A*. This finite size effect has a significant impact on the variance of surface tension values calculated using the virial summation method. A theoretical derivation of the origin of the effect is proposed, giving a new insight into the importance of optimising system dimensions in interfacial simulations. We demonstrate the consequences of this finite size effect via a new way to estimate the surface energetic and entropic properties of simulated air-liquid interfaces. Our method is based on macroscopic thermodynamic theory and involves comparing the internal energies of systems with varying dimensions. We present the testing of these methods using simulations of the TIP4P/2005 water forcefield and a Lennard-Jones fluid model of argon. Finally, we provide suggestions of additional situations, in which this finite size effect is expected to be significant, as well as possible ways to avoid its impact.

## I. INTRODUCTION

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.^{16} 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.^{20} Therefore, more suitable methods for including long-range LJ interactions in systems containing inhomogeneous density distributions are desirable.^{21}

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)^{28} 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^{−2} 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^{29} 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^{−2} 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.^{6} 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^{30} (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.

## II. BACKGROUND

### A. Thermodynamic surface properties

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*,^{31}

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},^{31}

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,^{32} based on the KB derivation,^{30} 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\alpha $ is the *α* component of the particle’s Cartesian velocity **v**_{i}, and $rij\alpha $ and $Fij\alpha $ 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},^{31}

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 $\gamma $ 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.

### B. Finite size effects of Var[*γ*]

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* = 2*L*_{x}*L*_{y}. It is assumed that at equilibrium, the kinetic terms of the pressure tensor components are equivalent to *ρk*_{B}*T* 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 $j\u2009rij=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,

The expectation value for the contribution from the bulk phase $\chi 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 $\gamma =\chi 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 $\chi 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 $\chi B2$ grows with *n*_{B}, and since $\chi S2$ and $\gamma 2$ are independent of the bulk phase, Var[*γ*] will also grow in proportion to $\chi 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*.

### C. Var[*γ*] dependence on *L*_{l}/*A*

We first represent $\chi 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^{16} is typically achieved by employing a spherical cutoff radius of $rc\u2264minLx,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 $\rho 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.

### D. Surface thermodynamics via physical inference

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^{33} for water and from the work of Lemmon and Penoncello^{34} for argon. The same procedure is performed on various experimental and simulated data sets reported in the literature.

#### 1. Water

Vargaftik, Volkov, and Voljak^{33} suggested the following equation as an amendment to Guggenheim^{35} 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),

#### 2. Argon

Lemmon and Penoncello^{34} 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),

### E. Surface thermodynamics via linear regression

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*_{0} and surface phases *U*_{s}*A*, 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^{27} 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 $\gamma $ to be dependent on *A*/*N*, although considering that (*A*/*N*)^{−1} ∝ *L*_{l}/*A*, we expect Var[*γ*] to differ between simulation dimensions. The average $\gamma $ 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)^{14} 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.

### F. Long-range corrections

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^{21} 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.^{21} 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 $\rho (z)$.

### G. Estimation of uncertainty

Ensemble averages of internal energies $U$ and surface tension $\gamma $ 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 $\tau AC/\tau $, where $\tau AC$ is the correlation time or minimum time length between independent samples. If measurements are uncorrelated, then $\tau AC=1$ and the variance is directly proportional to *τ*^{−1},

Therefore, the standard error in the mean of both properties needs to be amended by $\tau 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(\tau B)=\tau BVar\u0100B/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 $\tau B\u22121$ and $f(\tau B)\u22121$ so that the linear regression of both would yield an intercept of $\tau AC$.

## III. SIMULATION METHODOLOGY

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.^{43} This includes a B-spline Smoothed PME (SPME)^{44} 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^{−1}, whereas steps involving constant pressure used a Berendsen barostat. Our simulation procedure is then as follows:

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.Heat the cuboid system from 0 to T K over a period of 0.1 ns.

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

Add a vacuum region along the

*z*axis, adhering to the dimensions in Fig. 1, where*L*_{v}= 2*L*_{l}≥ 100 Å.Equilibrate the liquid-vapour system under NVT at T K conditions for 0.5 ns with Janeček’s trajectory corrections.

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} = 2*L*_{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 $\rho (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,^{21} 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.

. | LJ . | TIP4P/2005 . |
---|---|---|

. | Argon . | Water . |

M_{m} (g mol^{−1}) | 39.948 | 18.016 |

ϵ (kJ mol^{−1}) | 0.9937 | 0.7749 |

σ (Å) | 3.40 | 3.1589 |

r_{OH} (Å) | 0.9572 | |

r_{Oq} (Å) | 0.1546 | |

q_{H} (e) | 0.5564 | |

q_{q} (e) | −1.1128 | |

θ (deg) | 104.52 | |

φ (deg) | 52.26 |

. | LJ . | TIP4P/2005 . |
---|---|---|

. | Argon . | Water . |

M_{m} (g mol^{−1}) | 39.948 | 18.016 |

ϵ (kJ mol^{−1}) | 0.9937 | 0.7749 |

σ (Å) | 3.40 | 3.1589 |

r_{OH} (Å) | 0.9572 | |

r_{Oq} (Å) | 0.1546 | |

q_{H} (e) | 0.5564 | |

q_{q} (e) | −1.1128 | |

θ (deg) | 104.52 | |

φ (deg) | 52.26 |

## IV. RESULTS

We first present data illustrating the finite size dependence of Var[*γ*] on *L*_{l}/*A* for surface tensions calculated via the KB methodology.

### A. Var[*γ*] vs *L*_{l}/*A*

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 $\gamma $ for each cell dimensional ratio *L*_{l}/*A*. We see that there is a slight oscillation and a general increase in the uncertainty of $\gamma $ with *L*_{l}/*A*, both of which are most apparent for argon (Fig. 4). An oscillation in $\gamma $ has been previously reported for both KB^{24,36,37} and TA^{15} 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 *σ*.^{2} 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 $\gamma $, which can be problematic for systems with too small dimensions.^{23}

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.

### B. Surface thermodynamics

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.

Water 298 K . | Method . | r_{c} (Å)
. | γ (mJ m^{−2})
. | U_{s} (mJ m^{−2})
. | S_{s} (mJ m^{−2} K^{−1})
. |
---|---|---|---|---|---|

Experimental | |||||

Vargaftik, 1983 | … | … | 72.0 (0.1) | 118.0 (0.1) | 0.154 (0.001) |

TIP4P/2005 | |||||

This work | KB | 10.0 | 67.4 (1.8) | 108.3 (4.0) | 0.137 (0.008) |

KB/LR | 10.0 | 67.9 (0.6) | 109.7 (0.3) | 0.140 (0.002) | |

Ghoufri, 2008 | KB | 9.8 | 70.5 (8.9) | 111.4 (22.1) | 0.137 (0.046) |

Mountain, 2009 | KB | 12.6 | 68.3 (2.2) | 104.6 (5.0) | 0.122 (0.010) |

Alejandre, 2010 | KB | 9.5 | 69.1 (2.6) | 111.6 (6.1) | 0.143 (0.012) |

Sakamaki, 2011 | KB | 9.5 | 69.0 (2.2) | 112.1 (5.1) | 0.145 (0.011) |

Hu, 2015 | KB | ∞ | 68.5 (1.0) | 110.7 (2.4) | 0.142 (0.005) |

Lau, 2015 (293 K) | TA | 15.0 | 68.4 (0.5) | 134.7 (21.7) | 0.226 (0.074) |

Pascal, 2014 | 2PT | 12.0 | 64.4 (0.7) | 107.5 (2.5) | 0.145 (0.003) |

Water 298 K . | Method . | r_{c} (Å)
. | γ (mJ m^{−2})
. | U_{s} (mJ m^{−2})
. | S_{s} (mJ m^{−2} K^{−1})
. |
---|---|---|---|---|---|

Experimental | |||||

Vargaftik, 1983 | … | … | 72.0 (0.1) | 118.0 (0.1) | 0.154 (0.001) |

TIP4P/2005 | |||||

This work | KB | 10.0 | 67.4 (1.8) | 108.3 (4.0) | 0.137 (0.008) |

KB/LR | 10.0 | 67.9 (0.6) | 109.7 (0.3) | 0.140 (0.002) | |

Ghoufri, 2008 | KB | 9.8 | 70.5 (8.9) | 111.4 (22.1) | 0.137 (0.046) |

Mountain, 2009 | KB | 12.6 | 68.3 (2.2) | 104.6 (5.0) | 0.122 (0.010) |

Alejandre, 2010 | KB | 9.5 | 69.1 (2.6) | 111.6 (6.1) | 0.143 (0.012) |

Sakamaki, 2011 | KB | 9.5 | 69.0 (2.2) | 112.1 (5.1) | 0.145 (0.011) |

Hu, 2015 | KB | ∞ | 68.5 (1.0) | 110.7 (2.4) | 0.142 (0.005) |

Lau, 2015 (293 K) | TA | 15.0 | 68.4 (0.5) | 134.7 (21.7) | 0.226 (0.074) |

Pascal, 2014 | 2PT | 12.0 | 64.4 (0.7) | 107.5 (2.5) | 0.145 (0.003) |

Argon 85 K . | Method . | r_{c} (Å)
. | γ (mJ m^{−2})
. | U_{s} (mJ m^{−2})
. | S_{s} (mJ m^{−2} K^{−1})
. |
---|---|---|---|---|---|

Experimental | |||||

Various authors | … | … | 13.1 (0.1) | 34.8 (0.2) | 0.255 (0.001) |

LJ argon | |||||

This work | KB | 10.0 | 16.3 (1.0) | 40.1 (1.3) | 0.280 (0.010) |

KB/LR | 10.0 | 16.2 (0.2) | 38.2 (0.2) | 0.259 (0.003) | |

Trokhymchuk, 1999 | KB | 18.7 | 15.1 (2.0) | 36.7 (2.5) | 0.253 (0.020) |

Baidakov, 2000 | KB | 23.1 | 13.5 (0.7) | 36.4 (0.9) | 0.269 (0.008) |

Gloor, 2005 | KB | 18.7 | 15.2 (0.4) | 38.4 (0.5) | 0.273 (0.004) |

Goujon, 2014 | KB | 18.0 | 15.4 (0.7) | 39.1 (0.9) | 0.280 (0.008) |

Argon 85 K . | Method . | r_{c} (Å)
. | γ (mJ m^{−2})
. | U_{s} (mJ m^{−2})
. | S_{s} (mJ m^{−2} K^{−1})
. |
---|---|---|---|---|---|

Experimental | |||||

Various authors | … | … | 13.1 (0.1) | 34.8 (0.2) | 0.255 (0.001) |

LJ argon | |||||

This work | KB | 10.0 | 16.3 (1.0) | 40.1 (1.3) | 0.280 (0.010) |

KB/LR | 10.0 | 16.2 (0.2) | 38.2 (0.2) | 0.259 (0.003) | |

Trokhymchuk, 1999 | KB | 18.7 | 15.1 (2.0) | 36.7 (2.5) | 0.253 (0.020) |

Baidakov, 2000 | KB | 23.1 | 13.5 (0.7) | 36.4 (0.9) | 0.269 (0.008) |

Gloor, 2005 | KB | 18.7 | 15.2 (0.4) | 38.4 (0.5) | 0.273 (0.004) |

Goujon, 2014 | KB | 18.0 | 15.4 (0.7) | 39.1 (0.9) | 0.280 (0.008) |

#### 1. Surface tension

For TIP4P/2005 water, we report KB/LR *γ* = 67.9 ± 0.6 mJ m^{−2} in comparison with KB *γ* = 67.4 ± 1.8 mJ m^{−2} 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 $\gamma $ estimates at values of *T* < 400 K are consistently lower than those reported from other investigations by up to 2 mJ m^{−2}. 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^{−2} in comparison with KB *γ* = 16.1 ± 1.0 mJ m^{−2} 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^{21} 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,^{51} 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.^{41} 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 $\gamma $ 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 $\gamma $ or Var[*γ*] beyond a mesh displacement of 1 Å, in agreement with previous investigations.^{56} It is worth noting that an equally weighted mean of $\gamma $ 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.

#### 2. Surface energy

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^{2} 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^{−2} in comparison with KB *U*_{s} = 108.3 ± 4.0 mJ m^{−2} at 298 K. The difference between both estimates (+1.4 mJ m^{−2}) is significantly greater than that of *γ* (+0.5 mJ m^{−2}) 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^{−2} lower than the experimental value of *U*_{s} = 118.0 ± 0.1 mJ m^{−2}, 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,^{12} 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^{−2} in comparison with KB *U*_{s} = 40.1 ± 1.3 mJ m^{−2} at 85 K. Again the difference between these measurements (+1.9 mJ m^{−2}, approximately 3/2 the standard error in KB *U*_{s}) is greater than that between *γ* (+0.1 mJ m^{−2}). 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^{−2}, 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^{21} 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.

#### 3. Surface entropy

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^{−2} K^{−1} in comparison with KB *S*_{s} = 0.137 ± 0.008 mJ m^{−2} at 298 K. Again, despite a significant difference of around −0.015 mJ m^{−2} K^{−1} 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^{−2} K^{−1} in comparison with KB *S*_{s} = 0.280 ± 0.010 mJ m^{−2} K^{−1} at 85 K. We notice a large discrepancy here between both methodologies of 0.021 mJ m^{−2} K^{−1}. We believe this to be caused by the influence of Janeček’s^{21} 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^{−2} K^{−1}, 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.

## V. DISCUSSION

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.^{58} 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 $\gamma $ 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.^{27} 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*σ*^{−1}), 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.*^{64} 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.

## VI. CONCLUSION

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.

## SUPPLEMENTARY MATERIAL

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.

## ACKNOWLEDGMENTS

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.^{65}