Long range corrections (lrc) for the potential energy and for the force in planar liquid-vapor interface simulations are considered for spherically symmetric interactions. First, it is stated that for the Lennard-Jones (LJ) fluid the lrc for the energy Δu of Janeček [J. Phys. Chem. B 110, 6264 (2006)] is the same as that of Lotfi et al. [Mol. Simul. 5, 233 (1990)]. Second, we present the lrc for the force ΔF for any spherically symmetric interaction as a derivative of Δu plus a surface integral over the cut-off sphere by using the extended Leibniz rule of Flanders [Am. Math. Monthly 80, 615 (1973)]. This ΔF corrects the incomplete lrc Δ1F of Lotfi et al. and agrees with the result of Janeček obtained by direct averaging of the forces. Third, we show that the molecular dynamics (MD) results for the surface tension γ of the LJ fluid with size parameter σ obtained by Werth et al. [Physica A 392, 2359 (2013)] with the lrc ΔF of Janeček and a cut-off radius rc = 3σ agree with the results of Mecke et al. [J. Chem. Phys. 107, 9264 (1997)] obtained with the lrc Δ1F of Lotfi et al. and rc = 6.5σ within −0.4% to +1.6%. Moreover, using only the MD results for γ of Werth et al., we obtain for the LJ fluid a new surface tension correlation which also represents the γ-values of Mecke et al. within ±0.7%. The critical temperature resulting from the correlation is Tc = 1.317 66 and is in very good agreement with Tc,ref = 1.32 of the reference equation of state for the LJ fluid given by Thol et al. [J. Phys. Chem. Ref. Data 45, 023101 (2016)].
Early studies1–6 of the planar liquid-vapor interface of the Lennard-Jones (LJ) fluid were concerned with the structure of the interface and also presented results1,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 rc as is usually done in simulations of homogeneous systems; in Refs. 4 and 5, the cut-off radius is explicitly given as rc = 2.5σ with σ being the LJ size parameter. It was, however, pointed out in a Monte Carlo study of adsorption by Rowley, Nicholson, and Parsonage7 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 Fischer8 that the dew densities from interface simulations with rc = 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 Δ1F = ∇Δu was obtained by differentiation of Δu. By using the force correction Δ1F and enlarging rc 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.8 In a subsequent article by Mecke, Winkelmann, and Fischer,9 MD simulations were made for the LJ liquid-vapor interface using Δ1F in order to obtain improved orthobaric densities and also the surface tensions γ including a tail correction γtail. In that article,9 three setups with different particle numbers N and cut-off radii rc were used: (a) N = 1372, rc = 2.5σ, (b) N = 1372, rc = 5.0σ, and (c) N = 2048, rc = 6.5σ. It was found that for the low temperature T = 0.7 (reduced by kB/ε with kB 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 rc = 2.5σ to rc = 5.0σ is somewhat surprising and indicates only a weak effect of the lrc Δ1F. 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 rc = 6.5σ including γtail agree within −0.4% to +1.6% with recent MD results obtained by Werth et al.10 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 Fischer11 also performed MD simulations for liquid-vapor interfaces of the LJ mixture argon + methane. In view of the weak effect of the lrc Δ1F in Ref. 9 such a correction was not used anymore in the mixture simulations11 but a rather large cut-off radius rc = 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 article12 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ček12 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 Δ1F = ∇Δ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 Fischer8 and by Janeček12 are discussed in detail. In particular, we explore the mathematical reason for the incomplete lrc Δ1F 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 Fischer9 with the recent MD results obtained by Werth et al.10 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 Tc with the value from the recent reference equation of state of the LJ fluid.20
II. LONG RANGE CORRECTIONS
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/σ, spatial coordinate perpendicular to the interface = z/σ, energy = u/ε, force component Fz* = Fz/(ε/σ), temperature = kBT/ε, density ρ* = ρσ3, and surface tension γ* = γσ2/ε. For convenience, the stars are omitted where no confusion can occur.
In molecular simulations, a cut-off radius rc 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 Fischer8 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 r1 as
where u(r12) is the intermolecular potential between particles at r1 = (x1, y1, z1) and r2 = (x2, y2, z2) with r12 = r2 − r1,r12 = |r12| and z12 = z2 − z1. For the special case of the LJ potential, Eq. (2) yielded8 for the lrc Δu(z1) at a position z1,
Several years later Janeček12 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ček12 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ček12 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 ξ ≤ rc, 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 Fischer9 which definitely is not the case.
The merit of Janeček,12 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 Fz(z1) at a position z1 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,
where u(rc) is the LJ potential at r = rc.
The difference between the result of Ref. 12 and the result of Ref. 8 is in the additional term Δ2Fz(z1) which accounts for the lrc contribution at z1 from those particles which are outside the cut-off radius rc and have a z-coordinate z2 with –rc ≤ z12 ≤ rc (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 Δ1Fz(z1) from Ref. 8 and ΔFz(z1) 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.13 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 Dt that moves with the flow and a function G(x, t) on the region of flow. For that case, Flanders13,14 proved the following formula which was already known in continuum mechanics:15,16
Here = ndσ = (nx, ny,nz)dσ is the outward directed vectorial area element on the closed surface ∂Dt of the domain Dt. 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(r12), then G(x, t) is u(r12)ρ(r2) with r12 = |r2 − r1|, r1 = (x1 = 0, y1 = 0, z1) and r2 = (x2, y2, z2), t = z1, x = r2, the domain Dt is given by r12 > rc, and the velocity v is the unit vector in the z-direction v = ez. Therewith v· = nz 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 −Δ1Fz(z1) from Ref. 8 given in the present Eq. (4), whilst the first integral on the rhs corresponds to the negative lrc −ΔFz(z1) 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 nz = −z12/rc. Therewith one obtains
Finally, changing the integration over the surface area dσ to the integration over z2 by using dσ = 2π rcdz2 and replacing dz2 by dz12 yields
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 ΔFz(z1) in integral form as
III. SURFACE TENSION EQUATION FOR THE LENNARD-JONES FLUID
After the disagreement between Δ1Fz(z1) from Ref. 8 and ΔFz(z1) 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 Fischer9 from which we take the results obtained with 2048 particles, the lrc Δ1Fz(z1), and rc = 6.5σ. The other source is the publication of Werth et al.10 from which we take the results obtained with 300 000 particles, the lrc ΔFz(z1), and rc = 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.
|T .||γS1 .||γS2 .||ΔγS (%) .||γC1 .||γC2 .||γC3 .|
|1.10||0.3150 (122)||0.310 (4)||1.61||0.3150||0.3133||0.3129|
|T .||γS1 .||γS2 .||ΔγS (%) .||γC1 .||γC2 .||γC3 .|
|1.10||0.3150 (122)||0.310 (4)||1.61||0.3150||0.3133||0.3129|
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 Kohnstamm18
In Ref. 9, the parameters A, Tc, and b were obtained by a fit to the available three simulation results which gave A = 2.960 19, Tc = 1.325 21, and b = 1.264 15. In Ref. 10, the parameter Tc was taken from an external source19 as Tc = 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 Tc = 1.317 66 is in very good agreement with the critical temperature Tc,ref = 1.32 in the reference equation of state for the Lennard-Jones fluid.20
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,13 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.