Surface thermodynamics of planar, cylindrical, and spherical vapour-liquid interfaces of water

Surface thermodynamics of planar, cylindrical, and spherical vapour-liquid interfaces of water Gabriel V. Lau,1 Ian J. Ford,2 Patricia A. Hunt,3 Erich A. Müller,1 and George Jackson1 1Department of Chemical Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ, United Kingdom 2Department of Physics and Astronomy and London Centre for Nanotechnology, University College London, Gower Street, London WC1E 6BT, United Kingdom 3Department of Chemistry, Imperial College London, South Kensington Campus, London SW7 2AZ, United Kingdom


I. INTRODUCTION
Interfaces are omnipresent in real materials and understanding how systems behave at boundaries is key to technological control and manipulation at the smallest scales. Theories of the mechanics and thermodynamics of interfaces have been under development for over two centuries, 1 but only since the advent of direct molecular simulation in the late 1950s has it been possible to investigate interfaces from the microscopic perspective. In our current paper, we examine the interfacial properties of vapour-liquid interfaces by molecular dynamics simulation and examine the effect of different surface geometries (planar, cylindrical, and spherical) on the vapour-liquid interfacial tension and related thermodynamic properties such as the energy and entropy associated with these interfaces. The focus is on the interfacial properties of water as the most ubiquitous of fluids.
One of the key phenomena where the surface tension of water plays an important role is the process of atmospheric particle nucleation, which involves the formation of condensed-phase molecular clusters from the supersaturated vapour phase. [2][3][4][5][6] Within the framework of classical nucleation theory (CNT), the nucleation rate, which represents a measure of how quickly these clusters form, is very sensitive to the value of the vapour-liquid interfacial tension. Calculating the correct interfacial tension (or more precisely the surface free energy) for water clusters is therefore crucial to developing successful theories of nucleation, although one may rightly question the validity of a macroscopic description of this type because concepts such as the surface tension are not expected to apply at small length scales. A further complication is that the surfaces of nanoscale water clusters are clearly not planar, and as a result one also has to consider the effect of curvature on the interfacial tension.
Methods for the determination of the interfacial tension of fluids in molecular simulation (Monte Carlo or molecular dynamics) at the microscopic level fall into three general classes. 7 In the first, and by far the most common, approach, one takes a mechanical route where the surface tension is computed by evaluating the appropriate components of the pressure tensor. [8][9][10][11][12][13] For a planar configuration, the interfacial tension γ can be obtained from the condition of mechanical equilibrium as an integral in the difference of the normal P n = P bulk and tangential P t components of the pressure tensor across the interface The second class of approaches follows a thermodynamic perspective and involves the computation of the difference in free energy for systems in macrostates with different interfacial areas. 7,[14][15][16] The third and final class of methods is based on the concepts of finite-size scaling (FSS) within the grand canonical ensemble. [17][18][19] We employ a specific type of thermodynamic approach in our current work, where a perturbation in the mean interfacial area of the system is performed to determine the corresponding change in free energy. In this so-called test-area (TA) method, 7 the difference in free energy between a reference system of interest and a perturbed system of different mean interfacial areas can be expressed as a ratio of their partition functions, which reduces to a ratio of configurational phasespace integrals for isothermal perturbations. A thermodynamic expression for the interfacial tension can then be written as the following ensemble average: where ∆A is the difference in mean interfacial area (suitably defined) and ∆U is the difference in configurational energy between the perturbed and reference systems, T is the temperature, and k B is the Boltzmann constant. The angled brackets denote a canonical average in the ensemble associated with the reference system. The surface tension of the planar vapour-liquid interface of water has been studied extensively by computer simulation (e.g., see Refs. [20][21][22][23][24][25][26]. The procedure typically involves placing a film of liquid in contact with vapour and computing the surface tension using the mechanical (pressure tensor) route. The value of the vapour-liquid interfacial tension reported in the literature can vary significantly due to several factors such as the system size and lengths of the runs (finite-size effects become increasingly important when the interfacial area decreases), and details of the intermolecular potential model such as the cutoff range can significantly affect the computed properties.
The TA method has also been employed to obtain the interfacial tension of the planar vapour-liquid interface of water. Of particular mention is the thorough study by Vega and de Miguel 25 of the interfacial tension for the most popular classical non-polarizable potential models of water using both the traditional mechanical route and the thermodynamic TA approach. They found that the values of the interfacial tension computed from the mechanical and thermodynamic routes are consistent with each other. The TIP4P/2005 model 27 has been shown to provide the most accurate overall description of the vapourliquid interfacial tension for temperatures ranging from the triple to the critical point, the main reason for the choice of this particular model in our current work: at ambient temperature (300 K), the value of the tension obtained for the TIP4P/2005 model is γ ∞ = 69.3 mN m −1 which compares favourably with the experimental value of 71.7 mN m −1 for water.
The theoretical treatment of the thermodynamics of curved interfaces poses a greater challenge than that of planar interfaces. 28,29 Its history dates back to the seminal macroscopic mechanical description of Young 30 and Laplace 31 at the turn of the nineteenth century, followed by the formal thermodynamic framework of Gibbs 32 and Tolman. 33 Young and Laplace noted that the pressure difference between the inside and outside of a spherical liquid drop in equilibrium is proportional to the tension and inversely proportional to its radius where P l and P v are the pressures of the liquid inside and the vapour outside the drop, respectively, and R is the radius of the drop. Kelvin 34 subsequently employed the Young-Laplace description to relate the saturation properties associated with a liquid drop to those of the planar interface where P v (R) and P sat v are the saturated vapour pressures over a drop of radius R and a planar surface, respectively, γ ∞ is the interfacial tension associated with the planar vapour-liquid interface, M is the molecular mass, and ρ l is the mass density of the liquid. When the dimension of the drop is decreased to nanoscopic length scales, macroscopic concepts involving the bulk properties of the liquid such as the tension, pressure, and density become increasingly inappropriate. The macroscopic mechanical treatment of Young, Laplace, and Kelvin is nevertheless still commonplace in the description of curved liquid interfaces.
Recognising the importance of the effect of curvature on the interfacial tension of fluids, Tolman 33 followed the purely thermodynamic treatment introduced by Gibbs 32 to derive an expression for the dependence of the tension γ on the drop size where R s is the radius of the spherical "surface of tension" (defined as the surface where the Young-Laplace Eq. (3) holds exactly 1,35 ), and δ ∞ is the so-called Tolman length (in the planar limit). The common implementation of the Tolman relation is as a series expansion in powers of the curvature (inverse radius), though the physical significant of terms beyond leading order have been brought into question; 1,29 furthermore, the existence of non-analytical contributions have also be suggested. 36,37 In any case, the magnitude of the Tolman length has been shown to be small and difficult to determine even for the simplest of fluids (see Refs. [38][39][40]. There has been much debate and controversy on the nature of the curvature dependence of the surface tension of drops (the reader is directed to Ref. 29 for a relatively recent account). Gibbs argued that curvature corrections are only significant when drops approach molecular length scales and hypothesized that the tension decreases monotonically with decreasing radius. 41 On the other hand, Thomson and Thomson 42 suggested a nonmonotonic behaviour as a possibility, while Bakker 43 insisted that the tension should be independent of the radius.
Literature on the determination of the surface tension of drops and other curved fluid interfaces using molecular simulation is much less extensive than for planar interfaces. One of the earliest molecular investigations of the vapourliquid interfacial tension of drops was by Rusanov and Brodskaya 44 who studied clusters of spherically truncated and shifted Lennard-Jones (STS-LJ) particles and extracted the normal and tangential components of the Irving-Kirkwood 45 pressure tensor and also used the Young-Laplace equation to obtain the tension. The interfacial tension was found to decrease with decreasing radius of the drop for the cluster sizes studied (R = 3.75σ to 8σ, where σ is the particle diameter). A number of subsequent studies also employed a mechanical route involving an integration of the gradient of the normal component of the pressure tensor from the centre of the drop to the bulk vapour phase. 13,38,39,46 In particular, the work of Vrabec et al. 13 offers a very comprehensive study of the properties of planar and spherical interfaces of STS-LJ systems. They studied spherical droplets for radii between R = 5σ and 16σ and obtained the vapour-liquid interfacial tension using the mechanical route. The tension was again found to increase monotonically with the equimolar radius but only very slowly, not reaching the planar limit even for the largest drops considered. An important point to note is that the standard mechanical route to the surface tension becomes problematic for spherical interfaces due to the non-uniqueness of the components of the pressure tensor arising from the choice of contour joining two interacting particles, which results in an ambiguity in its definition. As pointed out early on by Schofield and Henderson, 47 the non-uniqueness associated with the local pressure tensor and the definition of the internal pressure for drops of high curvature will lead to a marked discrepancy in the free energies of formation determined with the mechanical and thermodynamic routes as has been observed in simulation. 48 The cause of the discrepancy between the values of the vapour-liquid interfacial tension obtained using the mechanical and thermodynamic routes has been examined recently with a TA simulation study of spherical drops of STS-LJ fluids. 49 In accord with the findings of the thermodynamic analysis in the grand canonical ensemble by Binder and coworkers, 18,19 the vapour-liquid interfacial tension was found to increase very slightly above the planar value on increasing the radius of the drop to R ∼ 8σ, followed by a slow decay to the planar limit, which would corresponds to a small negative Tolman length. 49 For a planar vapour-liquid interface, the change in free energy associated with a deformation in the interfacial area was shown to be fully characterized by the average of the change in configurational energy; this leadingorder contribution to the surface tension is found to be equivalent to that obtained with the mechanical relation (cf. Eq. (1)). However, in the case of liquid drops, a large second-order contribution to the tension (of the same order of magnitude as the leading-order term) was found to be associated with fluctuations in energy, invalidating the use of the mechanical relation for spherical geometries. The effect of fluctuations in the thermodynamics of drops was attributed to an additional entropic contribution. In our current paper, we will provide an alternative interpretation of the contributions to the free energy determined by the TA method, explaining the nature of the fluctuations and their fundamental connection with the geometry of the fluid structures being studied.
Investigations into the vapour-liquid interfacial tension of water droplets by molecular simulation are worth a particular mention as this is the specific subject of our current study. Samsonov et al. 50 used the Stockmayer potential parameterized to represent water and equated the change in free energy of formation of a drop to the change in internal energy accompanying its excision from a bulk parent phase. The surface tension determined in this way at 300 K was found to be a monotonically increasing function of the equimolar radius, starting from γ ∼ 20 mN m −1 for R ∼ 2 Å and reaching the asymptotic planar value for relatively small droplets corresponding to R < 10 Å. Ghoufi and Malfreyt 51 have performed mesoscale Monte Carlo simulations of water nanodrops with a DPD (dissipative particle dynamics) 52 model reparameterized to represent the surface tension at ambient conditions 26 using the local components of the pressure tensor; a monotonic behaviour for the curvature dependence of the surface tension was again found. Following an alternative thermodynamic approach, Joswiak et al. 53 employed the mitosis method to evaluate the free energy associated with separating a liquid drop into a pair of smaller drops, thus providing an estimate of the vapour-liquid interfacial tension of drops of water for the TIP4P/2005 potential model. Surprisingly, these authors found that the vapour-liquid interfacial tension at ambient conditions (300 K) now increased continuously from the planar limit of γ ∞ = 65.5 mN m −1 (for the TIP4P/2005 model not including long-range corrections 25 ) to γ ∼ 77 mN m −1 on decreasing the size of the nanodrop to a radius of R ∼ 6 Å. A similar increase in the interfacial tension of the TIP4P model by more than 10 mN m −1 from the planar limit on decreasing the radius to R ∼ 7 Å was reported recently by Homman et al. 54 using the Laplace relation with a novel approach for the determination of the surface of tension R s ; no discussion was offered as to why the trend is the opposite to that shown in earlier work by the group with the DPD model. 51 On the other hand, in their analysis of the vapour pressure of drops of water represented with the mW coarse-grained model, Factorovich et al. 55 demonstrated that the macroscopic Kelvin relation (cf. Eq. (4)) holds down to radii of R ∼ 6 Å (i.e., less than two molecular diameters) which suggests that the corresponding vapour-liquid interfacial tension is essentially insensitive to the curvature even at these small dimensions. Given the preceding discussion, it is clearly apparent that there is an inconsistency in the findings for the interfacial tension of water nanodrops worthy of further consideration.
In Sec. II, we discuss the basis of the test-area method and its connection with transformations of the spatial metric of the computational cell. We develop expressions for the surface tension, and its entropic and energetic components, in terms of the underlying perturbation of the configurational energy. For spherical geometries, we show that the quantities of interest depend on the scale of the distortion to second order, while in the case of planar and cylindrical geometries, the firstorder contributions are sufficient to characterize the interfacial thermodynamics. In Secs. III-V, we describe our numerical studies for the planar, cylindrical, and spherical geometries of liquid water structures of various sizes described with the TIP4P/2005 model to investigate the curvature dependence of the interfacial properties. The conclusions of our work are given in Sec. VI.

II. TEST-AREA METHOD
As has already been mentioned in Sec. I, the test-area method is a technique for calculating the change in free energy associated with a distortion of the system. In the simplest implementation, this distortion involves a volume-preserving modification to the dimensions of a periodically repeated cell, and the corresponding change in free energy can be associated with changes in the geometry of a condensed-phase system that extends through the periodic boundaries. We can employ this interpretation in studies of films of condensed matter that take the geometries of a planar slab or a cylindrical rod. These structures are mechanically stretched or compressed by the confining cell boundaries as result of the distortion, and the associated quasi-static work can be related to the difference in free energy brought about by the change in surface area of the structure. The test-area method can therefore be used to calculate the interfacial free energy, or surface tension.
As a simple example, one can consider a planar slab oriented in the x− y plane. A volume-preserving distortion of the lengths L i of each side of the cell may be defined by the transformations L x → L x (1 + ε) 1/2 , L y → L y (1 + ε) 1/2 , and L z → L z (1 + ε) −1 where ε is a dimensionless variable that characterizes the magnitude of the transformation. For a fluid system in a thermodynamic state within the binodal envelope, the structure adopted by N particles in an appropriate rectangular simulation prism is characterized by two vapour-liquid interfaces in the x− y plane with normals parallel to the z axis, as illustrated in Figure 1. The total interfacial area of the fluid in the original cell is therefore A(ε = 0) = 2L x L y and that of the distorted cell is A(ε) = 2L x L y (1 + ε), so that the corresponding change in area induced by the distortion is ∆A = 2L x L y ε. For this incremental change in area, the change in free energy corresponds completely to dF = γdA since the system is isothermal, and the number of particles and the volume of the box are also constant. In this case the vapourliquid interfacial tension may be expressed as where F(ε) is the free energy of the structure corresponding to the distortion ε.
A less intuitive interpretation of the TA method is necessary when it is applied to non-periodic structures, such as quasi-spherical droplets free to move in a rarefied vapour. Here, a volume-preserving distortion of the cell walls would not in itself be expected to influence the free energy of the system, and so the interpretation made for planar slab geometries would not be appropriate. Instead, we can take the distortion to be a volume-preserving change in the metric of the physical space, and therefore effectively in the strength of the interactions between the particles. We shall show that the corresponding change in free energy can be associated with a nominal change in the mean surface area of the droplet under the distortion, at least in the perturbative limit, and hence, the method again provides a route to the surface tension.
In Secs. II A and II B, we develop the second of these interpretations, in readiness for applications of the test-area method to slabs, cylinders, and droplets of condensed fluid in Secs. III-V, respectively. Expressions for the vapour-liquid interfacial tension are derived as well as the associated relations for the interfacial energies and entropies per unit area. The feasibility of evaluating these quantities numerically on the basis of the derived expressions is considered for the various geometries.

A. Interfacial free energy
We begin by calculating the partition function of a system of N particles in the original and distorted canonical cells, confining configurations to those that correspond to a slab in the x− y plane, i.e., those that pass through the periodic boundaries in the x−z and y−z planes, as illustrated in Fig. 2. The partition function (or more precisely the configurational integral) of the original system is Z =  dr exp(− βU(r)) over a phase space of particle coordinates r ≡ {x 1 , y 1 , z 1 , · · ·, z N } lying within the undistorted cell, where U is the potential energy of a configuration and β = 1/(k B T). The particle coordinates in the distorted cell can be written as a transformation of the coordinates in the original cell: for our current deformation G is a transformation matrix that depends on ε. With r and r ′ defined as vectors of dimension 3N, G is a 3N × 3N matrix that takes a block diagonal form consisting of N 3 × 3 matrices g defined by The partition function of the distorted system is Z ′ =  dr ′ exp(− βU(r ′ )) with the particle coordinates running over the distorted cell. Since the Jacobian of the transformation from r ′ to r is unity, we can express the partition function in terms of the coordinates of the undistorted cell as where the phase space is seen to be the same for the two partition functions. The potential energy of the distorted system can be defined as The difference in the free energies of the distorted and undistorted systems can then be written as where F = − 1 β ln Z, F ′ = − 1 β ln Z ′ , and the ensemble average, ⟨· · · ⟩, is taken over the undistorted reference system described by U.
The change in the free energy is clearly related to a change in interaction energy and we emphasize that this change can be associated with a modification to the metric of the space represented by the transformation matrix g. When evaluating Z ′ , the energy of a given configuration is to be determined using effective particle separations represented by ∆x ′ = g x x ∆x + g x y ∆ y + g x z ∆z, etc., where ∆x and ∆x ′ are particle separations in the x direction according to the original and distorted metrics.
Following the formalism of Ref. 49, one can develop an expansion in ∆U to give where ∆F i corresponds to the terms involving β i−1 . The moments of ∆U can be evaluated numerically in the limit ε → 0 to obtain the interfacial tension as γ = lim ε→ 0 (∆F/∆A), where ∆A is the generic change in area of the interface brought about by the distortion of the cell dimensions characterized by ε.
To better understand the nature of the TA methodology as applied to curved interfaces, one can express the change in the interaction energy of a molecular configuration r as a Taylor expansion in the perturbation parameter which satisfies ∆U = 0 when ε = 0 as required. The coefficients a and b are the first and second derivatives of ∆U with respect to ϵ, respectively (a = ∂∆U/∂ε and b = 1 2 ∂ 2 ∆U/∂ε 2 ). The dependence of ∆A on ε can then be considered. For the planar slab, we have already established that ∆A = 2L x L y ε, which implies that we need only consider the contributions to order ε in the expression for the change in free energy ∆F. From an inspection of Eqs. (11) and (12), it is apparent that the lowest-order terms in ∆F 2 and ∆F 3 are of order ε 2 and ε 3 , respectively, and therefore do not contribute to the interfacial tension for the planar geometry. The only relevant part of the change in configurational energy ∆U in Eq. (12) is the term that is linear in ε, and as a consequence and we conclude that the tension can be expressed simply as where Similarly, the change in surface area of a cylindrical structure aligned along the z axis brought about by the same where R c is the radius of the cylinder; this means that for small deformations the change in the exposed surface area of the cylindrical structure is again governed by contributions which are linear in ε. Following the same analysis as for the planar geometry, we find that once again the relevant contribution to the change in free energy is the first term ∆F 1 and the surface tension is where c 2 = −πR c L z . In Sec. IV, we scrutinize the significance of the radius R c in this context. It is important to note that if the cylinder is not aligned along the z axis, the said transformation will give rise to contributions in ε 2 as one finds for spherical structures.
The analysis is more complicated in the case of quasispherical geometries because the leading-order change in surface area of the structure is quadratic in ε: this is because distortions of a sphere into oblate (ε > 0) or prolate (ε < 0) spheroids both lead to an increase in the surface area. It can be shown that for the same spatial transformation as for the planar and cylindrical cases, the leading-order change in surface area s p /5 and R s p is the radius of the sphere, as will be shown later in Sec. V. Terms in ∆F up to second order in ε now become relevant. Indeed, there is no term proportional to ε due to the symmetrical nature of the transformation, which requires that the mean change in the energy of the system should be minimized at ε = 0, or equivalently that ⟨∂∆U/∂ε⟩ = ⟨a⟩ = 0 at ε = 0 in this system. Hence, ⟨∆U⟩ = ⟨b⟩ε 2 + O(ε 3 ) and we can write and while ∆F 3 remains third order in ε, and so The tension for a spherical geometry is then obtained by combining the first-and second-order terms as The dependence on second-order contributions (⟨b⟩ = 1/2⟨∂ 2 ∆U/∂ε 2 ⟩ and ⟨a 2 ⟩ = ⟨(∂∆U/∂ε) 2 ⟩) to the perturbation in the configurational energy suggests that the evaluation of the vapour-liquid interfacial tension of spherical liquid droplets will present a greater numerical challenge than the evaluation of the tension for planar and cylindrical interfaces, which only depends on the first-order energetic contribution ⟨a⟩.

B. Interfacial energies and entropies
The analysis can be extended to determine the separate energetic and entropic contributions to the interfacial tension for the three geometries considered. The entropic contribution is an important quantity to study in the context of the Gibbs adsorption relation where the surface excess entropy gives information about the temperature dependence of the tension.
. (20) On expanding the change in free energy as a series in energetic contributions ∆U and differentiating with respect to T, the corresponding changes in the entropy and energy can be expressed as and We are now in a position to extract the energetic and entropic contributions to the interfacial tension, e = lim ε→ 0 (∆E/∆A) and s = lim ε→ 0 (∆S/∆A) from Eqs. (21) and (22), respectively. For planar and cylindrical geometries, we need only evaluate the terms in Eqs. (21) and (22) that are proportional to ε, leading to and with i = 1 and 2 for planar slabs and cylinders, respectively. These leading order contributions were also derived early on for planar interfaces by Lekner and Henderson. 56 It is also useful to realise that even though the mechanical work ⟨∆U⟩ fully characterizes the tension for planar and cylindrical interfaces, this does not preclude the existence of a higher-order entropic contribution. The dependence on a difference between ⟨Ua⟩ and ⟨U⟩⟨a⟩ suggests that the computation of the separate energetic and entropic contributions to the surface tension will require good statistics in order to obtain reliable results. For a spherical interface, the expressions are once again different because the change in area ∆A does not have a linear term in ε. The terms proportional to ε in Eqs. (21) and (22) are expected to vanish due to the symmetrical nature of the transformation, and therefore we need to identify second-order contributions in ε 2 . We find that and Since the expressions involve second-order terms in ε, we expect the numerical evaluation of the energetic and entropic contributions to the surface tension to be much more challenging for a spherical interface than for planar or cylindrical geometries. The numerical scheme has to provide statistics that are good enough to distinguish ⟨Ub⟩ from ⟨U⟩⟨b⟩ (as well as ⟨Ua 2 ⟩ from ⟨U⟩⟨a 2 ⟩), which is considerably more difficult than the resolution of ⟨Ua⟩ and ⟨U⟩⟨a⟩ required for planar or cylindrical systems, cf. Eqs. (23) and (24), as b represents a higher-order perturbation to the system energy than a.
A positive energetic contribution to the interfacial tension is a reflection of the lower coordination of molecules at the interface compared with the bulk. The larger entropic contribution quantifies the increased freedom of movement of molecules at an interface compared with the bulk liquid; from the Gibbs adsorption equation, the entropy change per unit area can also be seen to be related to the temperature dependence of the surface tension, s = −∂γ/∂T, suggesting that the surface tension decreases with temperature. Furthermore, since the temperature dependence of the surface tension of liquids is quite often found to be nearly linear, 57,58 such that ∂γ/∂T ∼ −γ/T, the value of T s is expected to be of the same order as γ so that e would also be similar in magnitude.

III. PLANAR INTERFACE
We first investigate a planar interface of water as a base case, which has been studied extensively in the past and for which we expect to have the fewest problems both in terms of obtaining good statistics and in unambiguously defining changes in interfacial area. The atomistically detailed TIP4P/2005 model of water 27 is used to represent water in all of our simulations; TIP4P/2005 is a four-site model (with three point charges placed in a fixed planar configuration within a Lennard-Jones core) that has proved to be one of the most successful classical nonpolarizable force fields for the simulation of condensed water including the vapour-liquid interfacial tension. 25

A. Methodology
A liquid slab consisting of N = 1000 water molecules is placed between two empty regions in the centre of a periodic orthorhombic cell with dimensions L x = L y = 31.1 Å and L z = 100 Å (cf. Fig. 2). Molecular dynamics simulations are carried out in the canonical NVT ensemble at a temperature of T = 293 K, maintained using a Nosé-Hoover thermostat with a relaxation time of 0.1 ps and using a molecular dynamics timestep of 0.5 fs. 59 Interactions are cut off at 15 Å (roughly half the shortest box dimension) and the Ewald summation method is used to compute the long-range electrostatic interactions; no other long-range corrections are made for the computed properties. A notable point is the lack of particles in the vapour at the ambient conditions of the simulation. This is partially due to a limitations of the TIP4P/2005 water model, which has an electric dipole moment parameterized to bulk water rather an isolated gas-phase water molecule, leading to an underprediction of the saturated vapour pressure and vapour density. 60 Following an equilibration period to stabilize the mean total energy, running averages of the properties are collected over a 10 ns period. Changes in configurational energy, ∆U, are computed by applying virtual perturbations every 100 timesteps equivalent to the transformation L x → L x (1 + ε) 1/2 , L y → L y (1 + ε) 1/2 , L z → L z (1 + ε) −1 for separate values of ε from −0.005 to 0.005 in steps of 0.001 for all the simulations. It is important to reiterate that the transformation of the spatial metric involves not only a change in the box dimensions but that the relative positions between the centres of mass of the water molecules are scaled according to the affine transformation. The corresponding changes in free energy, ∆F(ε), are calculated directly using Eq. (10) and the changes in area are ∆A(ε) = 2L x L y ε. Since the leading order term in ∆F in this case is expected to be proportional to ε, the interfacial tension can be obtained from the difference [∆F(ε) − ∆F(−ε)]/(2c 1 ε) in the limit as ε → 0; this is equivalent to the form, 1 2 [∆F(ε)/∆A(ε) + ∆F(−ε)/∆A(−ε)]. Similar constructions employing ∆E(ε) and T∆S(ε) defined in Eqs. (22) and (21), respectively, can be used to extract the energetic and entropic contributions to the surface tension. The errors in the computed properties are estimated from the standard deviations of the means computed from 10 separate block averages of 1 ns.

B. Results
The TA methodology for the determination of the surface tension γ from area perturbations is illustrated in Figure 3: both positive and negative perturbations in the area are undertaken, and averages of the associated changes in free energy per change in area are shown. As expected, the contributions to the free energy are seen to converge in the limit ε → 0. A similar construction is followed in Figure 4 Table I where the energetic part γ energetic = e is computed according to γ energetic = γ + γ entropic . Our estimate of γ = 68.4 ± 0.5 mN m −1 at T = 293 K is consistent with the value of γ = 70.4 ± 1.3 mN m −1 obtained from the correlation proposed by Vega and de Miguel 25 for the tension of the TIP4P/2005 model determined with data from both the mechanical and thermodynamic TA routes (noting that these authors also employ long-range corrections); this is to be compared with the experimental value of 72.8 mN m −1 for water at 293 K. Furthermore, it is found that the first-order contribution ⟨∆U⟩ to the change in free energy fully characterizes the surface tension and the second-order contribution is essentially zero as consistent with the analysis in Sec. II A and with previous work on Lennard-Jones systems. 49 The value of the entropic contribution reported in Table I can be seen to be the same order of magnitude as the value of the tension, in line with our expectations.

IV. CYLINDRICAL INTERFACE
The curvature dependence of the vapour-liquid interfacial tension of cylindrical fluid interfaces is examined next by studying cylinders of TIP4P/2005 liquid water of varying radii. As with planar interfaces, the tension can be obtained from 1 2 [∆F(ε)/∆A(ε) + ∆F(−ε)/∆A(−ε)] in the limit as ε → 0 and similarly for the energetic γ energetic and entropic γ entropic contributions.

A. Methodology
Molecular dynamics simulations of cylindrical liquid structures consisting of between N = 64 and N = 1000 water molecules are carried out in a periodic cell of dimensions L x = L y = 100 Å and L z = 30 Å; this choice of rectangular prism facilitates the formation of a cylindrical structure with a principal axis along the z axis. As for the planar interface, interactions are cut off at 15 Å and the Ewald summation technique is used to compute the long-range electrostatic interactions. Initial configurations of cuboidal lattice configurations are allowed to relax and equilibrate into periodic cylinders, as illustrated in Fig. 5, until the mean total energy becomes time independent. This is followed by a 10 ns production run to collect statistics with a timestep of 0.2 fs for all simulations. A smaller timestep is chosen than for the case of the planar interface as cylindrical interfaces are more inherently unstable and have larger natural fluctuations in the thermodynamics properties. The smaller timestep also ensures a better overall conservation of energy for this system. The temperature is maintained at T = 293 K with the Nosé-Hoover thermostat, consistent with the conditions studied for the planar interface. In some circumstances (small system sizes), the cylinders cannot be stabilized for the entire 10 ns period and the structures relax to form drops; in these cases, independent simulations of shorter length (timescales where the cylinders remain mechanically stable) are carried out in parallel and the results averaged. Changes in configurational energy and the corresponding free energies are computed by applying the TA distortion transformation and scaling the particle coordinates in the same fashion as for the planar geometry.
The calculation of the associated changes in the surface area of the cylinders requires a knowledge of their radii which can be determined from an appropriate analysis of the density profiles: the mass-density profiles ρ(r) are obtained for each cylinder by counting the number of water molecules lying in cylindrical bins from the centre of mass and averaging over the entire trajectory. The data are correlated with a simple profile, 46 where ρ l is the liquid density, ρ v is the vapour density, and D is the interfacial width. The density profile is used to locate the equimolar (Gibbs) dividing surface, which corresponds to a vanishing adsorption in a one-component system, ensuring that the surface free energy per unit area so defined corresponds to the surface tension. 1 The radius R e of the equimolar dividing surface can then be calculated using numerical quadrature, 46 Though Eq. (28) for the equimolar radius was specifically developed for spherical geometries (which are treated separately in Section V), we employ this in our analysis of the cylinders for convenience as the differences are found to be very small; the resulting differences in the values of surface tension are within the statistical uncertainties of the simulation. The change in area corresponding to a distortion parameter ε is then obtained from ∆A(ε) = 2πR e L z , where the cylinder radius R c introduced in Sec. II A is now represented by the equimolar radius R e . The errors in the tension are estimated from standard deviations of the means computed from 10 separate block averages of 1 ns.

B. Results
The vapour-liquid interfacial tension of the cylindrical liquid structures formed by TIP4P/2005 water and the corresponding entropic contribution are depicted in Figure 6 as a function of cylinder radius R e ; the corresponding data are summarized in Table II. A weak curvature dependence in the tension is observed with a broad faint peak at a cylinder radius of R e ∼ 8 Å (though it should be recognized that the thinnest cylinders are also the most unstable); the values are seen to be very close to (and slightly above) the planar limit after taking into account statistical error and the uncertainty in defining the equimolar radius. The relative insensitivity of the vapourliquid interfacial tension to curvature for equilibrium cylindrical geometries of the LJ fluid has already been noted by El Bardouni et al. 61 As in the case of the planar interface, the change in free energy is found to be equal to the average of the change in configurational energy (∆F = ⟨∆U⟩), consistent with the analysis in Sec. II A.
No clear evidence of a curvature dependence of the entropic and energetic contributions to the tension for cylindrical geometries can be seen in Figure 6. Furthermore, it is apparent that the error in γ entropic is significantly larger than the error in the overall γ, which is negligible in comparison. A closer investigation of the entropic component of the interfacial tension of fluid cylinders formed by N = 216 TIP4P/2005 water molecules for 1 ns independent simulation blocks indicates that the mean values for each block fluctuate significantly, in contrast to the highly reproducible values for the full interfacial tension (see Figure 7). The fluctuations observed for this system are to be expected in view of the analysis carried out in Sec. II B, cf. Eq. (24), which demonstrates that the entropic contribution depends on quantities that are likely to be more susceptible to less extensive statistical sampling. This behaviour has been observed before, e.g., Fleischman and Brooks, 62 who found that the errors in the energetic and entropic components of the Helmholtz free energies of hydration for alkanes and alcohols were at least an order of magnitude greater than the error of the free energy itself. One can nevertheless conclude that the entropic and energetic contributions are of the same magnitude as γ, and much longer simulations would be required to obtain more reliable values with smaller statistical errors.

V. SPHERICAL INTERFACE
We now come to the central theme of our current paper: an in-depth analysis of the curvature dependence of the interfacial tension of quasi-spherical condensed water clusters of varying size. As was mentioned in the Introduction such an understanding is important in the modelling of nucleation; in classical nucleation theory, quantities such as the size of the critical cluster and its free energy of formation are characterized by a sensitive dependence on the surface tension. The capillarity approximation used in CNT, where the surface tension of small clusters is assumed to be that of a macroscopic planar interface, is most likely a contributor to the discrepancy between theoretical predictions and experimental nucleation rates which are often seen to be different by several orders of magnitude. 63

A. Methodology
Molecular dynamics simulations of TIP4P/2005 water droplets in equilibrium with vapour are performed for clusters consisting of between N = 16 and N = 1000 molecules in a cubic cell of dimensions L x = L y = L z = 270 Å. The initial configurations are constructed from cubic lattices in the centre of the cell that are allowed to equilibrate into spherical clusters until the mean of the total energy is essentially time independent. One observes immediately that the assumption that the drops are perfectly spherical is inaccurate, especially for the smaller clusters, as can be seen in Fig. 8. This intrinsic non-sphericity is problematic when considering statistics over short time scales because the "surface area" of the drop is not that of a sphere, the fundamental assumption made when computing the changes in area. The problem can be overcome by running sufficiently long simulations so that the assumption of a spherical cluster geometry holds true on average.
For the same reason as in the case of cylindrical interfaces, a relatively small timestep of 0.2 fs is used for all the simulations to ensure the stability of the droplets and suitable conservation of the total energy. The temperature of the system is again maintained at T = 293 K using the Nosé-Hoover thermostat. Pairwise interactions between all atomic sites are calculated by choosing a suitably large cutoff of 135 Å, precluding the need for periodic boundary conditions and Ewald summations. The configurational sampling is split into 50 independent runs (each 2 ns in length) executed in parallel corresponding to a total 100 ns worth of data. The velocity components of any water molecules (vapour phase) that reach a position which is 10 Å beyond the equimolar radius of the cluster are reversed. The velocity reversal and the choice of a very large cutoff effectively allow one to study isolated drops during the simulation. The neglect of the interaction between the condensed phase and the vapour beyond the reversal point is an appropriate approximation because the saturated-vapour density of the TIP4P/2005 model is very low at 293 K. The centre of mass of the drops is also repositioned to the centre of the simulation box every 1000 timesteps to minimize boundary effects. The procedure does not affect the energy of the system as the potential and kinetic energy of the target molecule are both conserved.
The changes in configurational energy are obtained using the TA distortion transformation specified for the planar interface, but with additional equivalent transformations in the other axes made possible by the (average) spherical symmetry of the system: L x → L x (1 + ε) −1 , L y → L y (1 + ε) 1/2 , L z → L z (1 + ε) 1/2 ; and L x → L x (1 + ε) 1/2 , L y → L y (1 + ε) −1 , L z → L z (1 + ε) 1/2 . A positive value of the perturbation parameter ε > 0 results in a deformation of the spherical drop into an oblate ellipsoid after the scaling of the particle coordinates, while a negative value of ε < 0 corresponds to a prolate ellipsoid.
Radial density profiles for each cluster are computed by counting the number of water molecules in spherical bins centred at the centre of mass and correlated to Eq. (27); Eq. (28) is then used to extract the equimolar radii R e . The differences in area between the perturbed oblate/prolate ellipsoids and the reference sphere can be calculated analytically, 64 ∆A oblate (ε) = 2πR 2 e (ε − 1) + πR 2 e ξ(1 + ε) 2 ln and where ξ is the ellipticity of the oblate ellipsoid defined as and the parameter α = arccos(1 − ε) 3 2 is convenient to characterize the change in area for the deformation into a prolate ellipsoid. One can show that both expressions for the change in area tend towards the form ∆A = c 3 ε 2 for small ε, with c 3 = 8πR 2 e /5. In the following analysis, we use the equimolar radius R e to represent the radius R s p of the sphere introduced in Sec. II A.
The leading-order terms in the change in free energy ∆F accompanying the TA distortion of a quasi-spherical geometry turn out to be quadratic in ε and so the interfacial tension is obtained from 1 2 [∆F oblate (ε)/∆A oblate (ε) + ∆F prolate (−ε)/ ∆A prolate (−ε)] or (for small ε) from the central difference [∆F(ε) + ∆F(−ε) − 2∆F(0)]/(2c 3 ε 2 ) in the limit as ε → 0. Similar expressions can be used for the energetic γ energetic and entropic γ entropic contributions to the interfacial tension. Standard deviations of the means computed for 50 separate block averages of 2 ns are employed to estimate the errors in the properties.

B. Results
The radial mass-density profiles determined for water clusters formed by the TIP4P/2005 model are presented in Figure 9. More substantial density fluctuations can be observed as the number of particles in the cluster is decreased. An oscillatory nature of the density is observed on the liquid side of the cluster as one approaches the centre of mass of the drop (particularly for the smaller radii). This could be partly due to a lack of statistical sampling in this region, but is expected for smaller clusters as a real consequence of the nanoscopic confinement of the system; cf. the non-local density functional theory calculations of liquid drops of the Lennard-Jones system. 29 Although Eq. (27) can be used to correlate the profiles (see Table III for the values of the coefficients), parameters such as ρ l do not necessarily coincide with the "bulk" liquid  III. Liquid and vapour densities, interfacial widths and equimolar radii obtained for the spherical clusters by correlating the density data with Eq. (27). density for such small structures. The average equimolar radii are also indicated on the figure for the various cluster sizes; R e is used to compute the change in area under the ellipsoidal distortion as described in Sec. V A.
The dependence of the vapour-liquid surface tension on the equimolar radius of the water drop is depicted in Figure  10, and the corresponding data are reported in Table IV. The tension computed with the TA deformations of TIP4P/2005 water clusters at T = 293 K is found to be a monotonically increasing function of the equimolar radius, tending to the planar limit γ → γ ∞ as R e → ∞. At the lower end of the scale, the tension is apparently found to vanish for clusters of N = 16 water molecules corresponding to R e ≈ 4.5 Å. In previous studies with Lennard-Jones systems, the vapour-liquid interfacial tension was found to exhibit a small maximum followed by a weak decay towards the planar value as R e was increased. 49 A similar maximum is not apparent for the TIP4P/2005 water drops, at least within the statistical uncertainty of our calculations. Although not explicitly presented here, the change in free energy for spherical water interfaces is found to comprise both first-order and second-order contributions in the configurational energy (∆F = ⟨∆U⟩ − β ), consistent with the analysis in Sec. II A.  Table I). As is apparent from Figure 10, our findings for the curvature dependence of the interfacial tension TIP4P/2005 water at 293 K are in general qualitative agreement with those of other studies for Stockmayer (dipolar Lennard-Jones) 50 and DPD 51 water nanodroplets. A direct comparison cannot easily be made because these authors employ different intermolecular potential models, run their simulations at the slightly higher temperatures of T = 298 K 51 and 300 K, 50 report a different planar limit, and do not quote error bars. The asymptotic planar value is reached at a significantly larger value of R e ∼ 15 Å (within error bars) based on our TA calculations for the TIP4P/2005 model than based on the calculations of the energy of excision for the Stockmayer model or the TA values for the DPD model. Our finding that the tension essentially vanishes, albeit with a large error, for a small system size of N = 16 water molecules (corresponding to R e ≈ 4.5 Å) is however inconsistent with the stable drops reported by Samsonov et al. 50 and by Ghoufi and Malfreyt 51 for systems of similar and even smaller sizes: Samsonov et al. 50 find a value of the tension which is very close to the planar limit at a radius of ∼4.5 Å while Ghoufi and Malfreyt 51 find a value which is ∼65% of the planar limit at the same radius. It is possible that this difference can be attributed to the different force fields considered.
The general overall trend obtained in our and the aforementioned 50,51 studies that the vapour-liquid interfacial tension of nanoscopic drops of water decreases quite sharply with a decrease in the drop radius after a given size is reached is in stark contrast with the recent findings of Joswiak et al. 53 for TIP4P/2005 water (the same potential model as we use) and of Homman et al. 54 for TIP4P water; the data for the TIP4P model have not been included in Figure 10 because the force field leads to an underestimate of the planar value of the tension of real water of almost 20%. 25 Both of these groups report a significant increase in the interfacial tension of water drops by more than 10 mN m −1 from the planar limit as the radius is decreased to R e ∼ 7.0 Å at ambient temperature (300 K). Joswiak et al. 53 computed the free energy associated with separating a liquid drop into a pair of smaller drops to estimate the interfacial tension, while Homman et al. 54 determine the tension from the surface of tension R s and the Laplace relation (with R = R s in Eq. (3)). It is difficult to uncover the reason for the discrepancy with our findings for the curvature dependence of water nanodrops, though it is possible that the use of a local pressure of the liquid at the interior of the drop by Homman et al., 54 rather than the pressure of an equivalent bulk liquid phase with the same chemical potential, could be problematic; the reader is directed to Ref. 29 for a discussion of the issues associated with the use of the Laplace relation for small systems. One should also recall at this stage that Factorovich et al. 55 do not find a significant curvature dependence for the vapour-liquid interfacial tension of the mW coarsegrained model of water determined from the Kelvin relation (cf. Eq. (4)) even for small drops with radii down to R e ∼ 6.0 Å, though the validity of the macroscopic Kelvin description can also be questioned at these microscopic length scales.
It is important to reflect on the statistical uncertainties in the reported interfacial tensions for the spherical structures. In Sec. II A, we demonstrate that the tension for spherical drops depends on second-order contributions to the configurational energy perturbation. Specifically, the change in free energy due to the distortion is derived from a mean and a fluctuation term in the statistics of the configurational energy change, both of which turn out to be of the same order of magnitude (and orders of magnitude larger than the value of the shift in free energy) but of opposite sign. The requirement of evaluating a difference between two large numbers means that extremely long runs are needed to provide sufficient statistics to extract an accurate value for the tension: typically, 100 ns simulation runs are required for each spherical structure in contrast with 10 ns for a planar slab. We illustrate this point in Figure 11 where fluctuations between estimates of the tension from separate 2 ns simulation runs are seen to range from negative values to values more than twice the average for the N = 216 cluster. The uncertainty in extracting a reliable average from this "noise" is represented by the error bars associated with the reported vapour-liquid interfacial tensions in Figure 10.
The entropic and energetic contributions to the surface tension discussed in Sec. II B for spherical structures involve a difference of second-order configurational energy correlations, and we find that our simulations are too limited in extent to allow for a reliable computation of these properties.

VI. CONCLUSIONS
We provide an interpretation of the test-area method in terms of a perturbation in the metric of the coordinate space and hence of the configurational energy, reminiscent of a distortion of the cell boundaries but more far-reaching. The change in energy ∆U(r) of a molecular configuration under a shift in the metric can be written as a Taylor expansion in a dimensionless parameter ε that characterizes the extent of the distortion. Expressions are then derived for the interfacial tension of condensed-phase structures for three different geometries: planar, cylindrical, and spherical. It is shown that when the change in free energy ∆F associated with the distortion is expanded in powers of ∆U, only the first moment of ∆U contributes to the tension of planar and cylindrical interfaces because the leading-order change in the surface area is linear in ε. As a consequence, the only part in ∆U that is relevant for the surface thermodynamics of planar and cylindrical interfaces is the term that is linear in ε. By contrast, the distortional change in area is quadratic in ε for spherical geometries, and terms in ∆F up to second order in ε are now relevant implying that the second moment of ∆U contributes to the tension.
The energetic and entropic contributions to the interfacial tension are also derived. For planar and cylindrical geometries, these quantities are shown to depend on first-order correlations between the configurational energy U and its change due to the distortion ∆U. In the case of spherical geometries, the contributions depend on both first-and second-order correlations, involving second moments of ∆U. Our analysis indicates that the calculation of the surface tension of spherical structures with the test-area perturbative approach will be more challenging than for planar and cylindrical cases. Furthermore, the extraction of entropic and energetic contributions to the tension will be more difficult than the task of calculating the surface tension itself for all geometries.
As an important application, we investigate different geometries of condensed structures formed by water represented with the TIP4P/2005 model at ambient temperature using molecular dynamics simulation. For the planar interface, our calculations of the vapour-liquid interfacial tension are consistent with literature values obtained with both the mechanical (pressure-tensor) route to the tension and the test-area method. Additionally, we find that the magnitude of the entropic contribution is similar to that of the tension, in line with our analysis. The surface entropy term has received little attention in previous studies. In the case of cylindrical interfaces formed by TIP4P/2005 water, the tension is computed for a range of cylindrical radii: no clear curvature dependence is observed, at least within the range of radii studied and the statistical uncertainty of the method. There is also no clear curvature dependence for the entropic and energetic contributions, but the magnitude of these contributions is again consistent with our analysis. On the other hand, we find a significant curvature dependence of the vapour-liquid interfacial tension for spherical clusters formed by TIP4P/2005 water: the tension is found to be a monotonically increasing function of the drop radius, taking a value of around zero at an average equimolar radius of ∼ 4.5 Å (corresponding to just N = 16 water molecules), and then approaching the planar limit for a radius of about ∼ 15 Å (corresponding to ∼5 molecular diameters). This is consistent with previous test-area simulations of Lennard-Jones drops 49 where the tension retains a value close to that of the planar limit for drops with equimolar radii down to ∼ 6 molecular diameters; in the case of the LJ system, however, a weak maximum in the tension is found which is not seen for water within the statistical uncertainty of our results. The curvature dependence that we obtain for the TIP4P/2005 model is qualitatively in line with some other studies on water nanodrops, 50,51 though in the latter work, the tension was found to approach the planar limit at much smaller radii (corresponding to 1-2 molecular diameters). The general trend uncovered here does not however support the findings reported in some recent work on water nanodrops 53,54 where the interfacial tension was found to increase from the planar limit with decreasing drop radius.
Statistical uncertainties are found to be much larger for the surface tension of spherical interfaces than for those of planar and cylindrical interfaces, as expected from the theoretical analysis. The importance of second-order contributions means that much longer simulations (by an order of magnitude) are required in the simulations of droplet to extract surface tensions, and the statistical error associated with our results is higher than that obtained from the shorter simulations of the planar and cylindrical cases. The entropic and energetic contributions for planar and cylindrical interfaces are also noted to be susceptible to limited statistical sampling. In the case of spherical interfaces, a computational effort far beyond that undertaken here would be needed to resolve the surface tension into these contributions with acceptable statistical accuracy.