Communication: Long range corrections in liquid-vapor interface simulations

Early studies^1–6 of the planar liquid-vapor interface of the Lennard-Jones (LJ) fluid were concerned with the structure of the interface and also presented results^1,2,4–6 for the surface tension γ. In the simulation studies,^2–5 the long range parts of the energies or forces were neglected beyond a cut-off radius r_c as is usually done in simulations of homogeneous systems; in Refs. 4 and 5, the cut-off radius is explicitly given as r_c = 2.5σ with σ being the LJ size parameter. It was, however, pointed out in a Monte Carlo study of adsorption by Rowley, Nicholson, and Parsonage⁷ that for inhomogeneous systems, the effects of the long range interactions do not cancel out but can be included by appropriate corrections. For the LJ liquid-vapor interface, it was found by Lotfi, Vrabec, and Fischer⁸ that the dew densities from interface simulations with r_c = 2.5σ are mostly too high by a factor of 3 in comparison with those obtained via bulk fluid simulations.

Hence, in Ref. 8 a long range correction (lrc) for the energy Δu was derived for the liquid-vapor interface and therefrom a lrc for the force Δ₁F = ∇Δu was obtained by differentiation of Δu. By using the force correction Δ₁F and enlarging r_c up to 5.0σ, direct molecular dynamics (MD) simulations of the LJ liquid-vapor interface with 1372 particles brought the orthobaric densities close to those obtained from bulk fluid simulations.⁸ In a subsequent article by Mecke, Winkelmann, and Fischer,⁹ MD simulations were made for the LJ liquid-vapor interface using Δ₁F in order to obtain improved orthobaric densities and also the surface tensions γ including a tail correction γ_tail. In that article,⁹ three setups with different particle numbers N and cut-off radii r_c were used: (a) N = 1372, r_c = 2.5σ, (b) N = 1372, r_c = 5.0σ, and (c) N = 2048, r_c = 6.5σ. It was found that for the low temperature T = 0.7 (reduced by k_B/ε with k_B being the Boltzmann constant and ε the LJ energy parameter) the surface tension γ increased by 7.6% in going from (a) to (b) and by 0.5% in going from (b) to (c). For the high temperature T = 1.10, the surface tension increased by 78% in going from (a) to (b) and decreased by 1.35% in going from (b) to (c). The strong variation of the results in going from r_c = 2.5σ to r_c = 5.0σ is somewhat surprising and indicates only a weak effect of the lrc Δ₁F. Hence, in Ref. 9, it was concluded, “In order to obtain reliable values for the surface tension, cut-off radii of at least 5 molecular diameters supplemented by a tail correction are required.” It will be shown below that the results for γ obtained in Ref. 9 with N = 2048 and r_c = 6.5σ including γ_tail agree within −0.4% to +1.6% with recent MD results obtained by Werth et al.¹⁰ which we believe to be presently the most reliable simulation results for the LJ liquid-vapor interface. It should still be mentioned that Mecke, Winkelmann, and Fischer¹¹ also performed MD simulations for liquid-vapor interfaces of the LJ mixture argon + methane. In view of the weak effect of the lrc Δ₁F in Ref. 9 such a correction was not used anymore in the mixture simulations¹¹ but a rather large cut-off radius r_c = 7.0 σ_Ar = 6.38 σ_CH4 was used for all interactions.

Several years after Refs. 8 and 9 had appeared, Janeček published an interesting article¹² on lrcs for the energy and the force in inhomogeneous simulations. The lrc for the energy given in Ref. 12 is the same as was already given in Ref. 8 which was overlooked by Janeček and subsequent authors. The merit of Janeček¹² is that he derived a lrc for the force ΔF by directly averaging over the forces from outside the cut-off sphere. Based on his result for ΔF, Janeček pointed out that the lrc for the force Δ₁F = ∇Δu derived in Ref. 8 and used in Refs. 8 and 9 is incomplete. From the physical point of view, this becomes immediately evident by looking at the upper part of Fig. 1 in Ref. 12.

The present paper is organized such that in Sec. II, the lrcs for the energy and for the force as given by Lotfi, Vrabec, and Fischer⁸ and by Janeček¹² are discussed in detail. In particular, we explore the mathematical reason for the incomplete lrc Δ₁F given in Ref. 8 by using an extended version of the Leibniz rule for a three-dimensional integral with a parameter.^13–16 In Sec. III, we compare the MD results for the surface tension γ obtained by Mecke, Winkelmann, and Fischer⁹ with the recent MD results obtained by Werth et al.¹⁰ Moreover, we suggest for the LJ fluid, an improved correlation for the surface tension based only on the MD data of Ref. 10 and compare the critical temperature T_c with the value from the recent reference equation of state of the LJ fluid.²⁰ 

We consider only spherically symmetric intermolecular potentials and a planar liquid-vapor interface. If not stated otherwise, we use the LJ potential u(r) with the energy parameter ε, the size parameter σ, and r being the intermolecular distance,

Henceforth the following reduced quantities are used: intermolecular distance $r*$ = r/σ, spatial coordinate perpendicular to the interface $z*$ = z/σ, energy $u*$ = u/ε, force component F_z* = F_z/(ε/σ), temperature $T*$ = k_BT/ε, density ρ* = ρσ³, and surface tension γ* = γσ²/ε. For convenience, the stars are omitted where no confusion can occur.

In molecular simulations, a cut-off radius r_c has to be introduced beyond which the potential or the force are set equal to zero. Hence, in order to obtain the interface properties of the liquid-vapor interface of the LJ fluid with the full potential given in Eq. (1), Lotfi, Vrabec, and Fischer⁸ derived lrcs for the potential energy and the force. Assuming a cylindrical coordinate system with the z-axis being perpendicular to the interface, they obtained the lrc Δu to the potential energy at a point r₁ as

where u(r₁₂) is the intermolecular potential between particles at r₁ = (x₁, y₁, z₁) and r₂ = (x₂, y₂, z₂) with r₁₂ = r₂ − r_1,r₁₂ = |r₁₂| and z₁₂ = z₂ − z₁. For the special case of the LJ potential, Eq. (2) yielded⁸ for the lrc Δu(z₁) at a position z₁,

Therefrom, these authors⁸ obtained by differentiation of Eq. (3) on the basis of the usual one-dimensional Leibniz rule for the lrc of the z-component of the force Δ₁F_z(z₁) as

which was used in Refs. 8 and 9.

Several years later Janeček¹² also derived lrcs for the potential energy and the force. Let us first consider his lrc for the potential energy given in his Eqs. (11) and (12). The only difference between his formulas and lrc for the potential energy of Ref. 8, given above as Eq. (3), is that Eq. (3) presents the lrc as one-dimensional integrals, whilst Janeček¹² uses summations over strips. But this formal difference should not have a relevant impact on the results. A point of concern, however, is that Janeček¹² did not mention the correct lrc for the energy from Ref. 8, the above Eq. (3), but has given a lrc for the potential energy in his Eq. (16) in which the upper line, w(ξ) = 0 for ξ ≤ r_c, is wrong. He claims that this lrc for the energy in Monte Carlo simulations is equivalent to the lrc for the force in Ref. 9. Unfortunately, this statement together with the caption of Fig. 2 in Ref. 12 can be misunderstood in the sense that the wrong Eq. (16) in Ref. 12 is due to Mecke, Winkelmann, and Fischer⁹ which definitely is not the case.

The merit of Janeček,¹² as already mentioned, is that he also gave a lrc for the force. We remind that in Ref. 8, the lrc for the force was derived from Δu by interchanging differentiation and integration following the usual Leibniz rule for a one-dimensional integral with a parameter. Janeček, however, used a different route. He first derived from the LJ potential u(r) the full force in the z-direction F_z(z₁) at a position z₁ resulting in the algebraic expression

Thereafter, he averaged the long range contributions of the forces similar as it is done in the above Eq. (2) for the potential energy so that he avoided the interchange of differentiation and integration. Thus he obtained for the lrc for the force,

with Δ₁F_z(z₁) given by Eq. (4). The second term Δ₂F_z(z₁) was given by Janeček¹² in his Eq. (13) and in the first line of Eq. (14) (for the case ξ ≤ r_c) and can be written in integral form as

where u(r_c) is the LJ potential at r = r_c.

The difference between the result of Ref. 12 and the result of Ref. 8 is in the additional term Δ₂F_z(z₁) which accounts for the lrc contribution at z₁ from those particles which are outside the cut-off radius r_c and have a z-coordinate z₂ with –r_c ≤ z₁₂ ≤ r_c (see the upper part of Fig. 1 in Ref. 12). As these particles have to contribute to the lrc for the force, the statement of Janeček is definitely correct that the lrc for the force used in Ref. 9 (taken from Ref. 8) is incomplete.

When it became clear that the lrc for the force derived in Ref. 8 is incomplete, the challenge arose to understand the difference between Δ₁F_z(z₁) from Ref. 8 and ΔF_z(z₁) from Ref. 12. For that purpose, we use the generalization of the Leibniz rule to the case of a three-dimensional space in the form given by Flanders.¹³ He considers a fluid flowing through a region of space. The Euler description gives the velocity v(x, t) at time t at position x. Suppose now a domain D_t that moves with the flow and a function G(x, t) on the region of flow. For that case, Flanders^13,14 proved the following formula which was already known in continuum mechanics:^15,16

Here $dσ$= ndσ = (n_x, n_y,n_z)dσ is the outward directed vectorial area element on the closed surface ∂D_t of the domain D_t. Graphical representations of the considered situation are given in Ref. 13 and in Ref. 14. Moreover we mention that in these references, the function G(x, t) is called F(x, t) which we changed in order to avoid confusion with the force terms in the present paper.

If we apply now Eq. (8) to our case with planar geometry and any spherically symmetric interaction u(r₁₂), then G(x, t) is u(r₁₂)ρ(r₂) with r₁₂ = |r₂ − r₁|, r₁ = (x₁ = 0, y₁ = 0, z₁) and r₂ = (x₂, y₂, z₂), t = z₁, x = r₂, the domain D_t is given by r₁₂ > r_c, and the velocity v is the unit vector in the z-direction v = e_z. Therewith v·$dσ$ = n_z dσ and in these notations the Leibniz rule takes the form

Let us first discuss the meaning of this equation for the LJ fluid. The expression on the lhs corresponds to the negative lrc −Δ₁F_z(z₁) from Ref. 8 given in the present Eq. (4), whilst the first integral on the rhs corresponds to the negative lrc −ΔF_z(z₁) from Ref. 12 given in the present Eq. (6). The difference between the two lrcs obtained in Ref. 8 and Ref. 12 is represented by the second integral on the rhs of Eq. (9). This can still be rewritten in simpler form where we have to keep in mind that the integration is made outside the cut-off sphere and hence the vector n is directed inside the cut-off sphere yielding n_z = −z₁₂/r_c. Therewith one obtains

Finally, changing the integration over the surface area dσ to the integration over z₂ by using dσ = 2π r_cdz₂ and replacing dz₂ by dz₁₂ yields

which corresponds to the term Δ₂F_z(z₁) obtained by Janeček¹² for the LJ-fluid and is given above in Eq. (7). Summarizing, we obtain from Eq. (9) after rearranging the sequence of the terms

We emphasize that Eq. (12) is valid for all spherically symmetric interactions u(r) and planar interfaces. It means that the lrc for the force can be obtained as a negative derivative of the lrc for the potential energy plus a term with an integral over the cut-off sphere.

Application of Eq. (12) to the LJ fluid yields the lrc for the force ΔF_z(z₁) in integral form as

which is in agreement with Eq. (6) in combination with Eqs. (4) and (7).

Equation (12) can, e.g., also be used to get the lrc for the force in integral form for the Mie n-m potentials which were considered recently in Ref. 17.

After the disagreement between Δ₁F_z(z₁) from Ref. 8 and ΔF_z(z₁) from Ref. 12 has been clarified in Sec. II, it seems interesting to compare representative MD results for the LJ fluid obtained with either lrc for the force. One source is the publication of Mecke, Winkelmann, and Fischer⁹ from which we take the results obtained with 2048 particles, the lrc Δ₁F_z(z₁), and r_c = 6.5σ. The other source is the publication of Werth et al.¹⁰ from which we take the results obtained with 300 000 particles, the lrc ΔF_z(z₁), and r_c = 3.0σ.

The simulation results are compiled in Table I. Column 2 shows the simulation results γ_S1 from Ref. 9 and column 3 shows the simulation results γ_S2 from Ref. 10. Direct comparison can be made for the temperatures T = 0.70 and T = 1.10. The relative differences Δγ_S = (γ_S1/γ_S2 − 1) × 100 are shown in column 4 and amount to −0.4% at the low temperature and 1.6% at the high temperature which is quite satisfactory in view of the uncertainties of other simulation results.

Another point of interest is the correlation equation for the surface tension for which we use the equation given in the book of van der Waals and Kohnstamm¹⁸ 

In Ref. 9, the parameters A, T_c, and b were obtained by a fit to the available three simulation results which gave A = 2.960 19, T_c = 1.325 21, and b = 1.264 15. In Ref. 10, the parameter T_c was taken from an external source¹⁹ as T_c = 1.3126 and A and b were obtained by a fit to the simulation results which gave A = 2.94 and b = 1.23. The surface tensions obtained from these two correlation equations are also contained in Table I as γ_C1 (Ref. 9) and γ_C2 (Ref. 10). Therefrom we learn that at the highest temperature T = 1.25, the correlation γ_C2 yields a considerably too low value for the surface tension which can be attributed to the chosen critical temperature. Surprisingly, at that temperature the correlation γ_C1 is closer to the simulation value γ_S2 than γ_C2. As a consequence, we made a least square fit using only all simulation data of Ref. 10, Tables 1 and 2, for 300 000 particles and obtained as a correlation equation,

The surface tensions obtained from this correlation equation are also contained in Table I as γ_C3. First, we see that this correlation reproduces the simulation results γ_S2 from Ref. 10 very well and also represents the γ_S1 values of Ref. 9 within ±0.7%. Moreover, we note that the obtained critical temperature T_c = 1.317 66 is in very good agreement with the critical temperature T_c,ref = 1.32 in the reference equation of state for the Lennard-Jones fluid.²⁰ 

First, it was stated that for the LJ fluid the lrc for the energy in Ref. 12 is the same as that given earlier in Ref. 8. Second, by using the Leibniz rule for the three-dimensional integral given by Flanders,¹³ we detected the reason for the incomplete lrc for the force derived from the lrc for the energy in Ref. 8. Third, we found agreement between the surface tension results from Ref. 9 obtained with the incomplete force-lrc and a cut-off radius of 6.5σ and those from Ref. 10 with the correct force-lrc within −0.4% to +1.6%. Finally, we obtained a correlation equation for the surface tension of the LJ fluid by taking only MD data from Ref. 10. This correlation reproduces also the MD results of Ref. 9 for the surface tension within ±0.7%. Moreover, the critical temperature resulting from the correlation is in very good agreement with that from the recent reference equation of state for the LJ fluid in Ref. 20.

The authors gratefully acknowledge their participation in the ESI/CECAM Workshop on “Physics and Chemistry at Fluid/Fluid Interfaces” held at the Erwin Schrödinger International Institute for Mathematics and Physics (ESI), Vienna, on December 11 to 13, 2017, which stimulated the present Communication. S.L. also acknowledges support from ESI. Moreover, the authors thank Professor Jadran Vrabec, TU Berlin, for fruitful discussions, and J.F. thanks Professor Wolfgang Wagner, Ruhr-Universität Bochum, for the “Lehrbuch der Thermodynamik” by van der Waals and Kohnstamm.