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 Δ

**corrects the incomplete lrc Δ**

*F*_{1}

**of Lotfi**

*F**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 Δ

**of Janeček and a cut-off radius r**

*F*_{c}= 3

*σ*agree with the results of Mecke

*et al.*[J. Chem. Phys.

**107**, 9264 (1997)] obtained with the lrc Δ

_{1}

**of Lotfi**

*F**et al.*and

*r*

_{c}= 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

*T*

_{c}= 1.317 66 and is in very good agreement with

*T*

_{c,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)].

## I. INTRODUCTION

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^{7} 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^{8} 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 Δ_{1}** F** = ∇Δ

*u*was obtained by differentiation of Δ

*u*. By using the force correction Δ

_{1}

**and enlarging**

*F**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.

^{8}In a subsequent article by Mecke, Winkelmann, and Fischer,

^{9}MD simulations were made for the LJ liquid-vapor interface using Δ

_{1}

**in order to obtain improved orthobaric densities and also the surface tensions**

*F**γ*including a tail correction

*γ*

_{tail}. In that article,

^{9}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 Δ

_{1}

**. 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**

*F**γ*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.*

^{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 Fischer

^{11}also performed MD simulations for liquid-vapor interfaces of the LJ mixture argon + methane. In view of the weak effect of the lrc Δ

_{1}

**in Ref. 9 such a correction was not used anymore in the mixture simulations**

*F*^{11}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^{12} 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^{12} 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 Δ

**, Janeček pointed out that the lrc for the force Δ**

*F*_{1}

**= ∇Δ**

*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^{8} and by Janeček^{12} are discussed in detail. In particular, we explore the mathematical reason for the incomplete lrc Δ_{1}** 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

^{9}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

*T*

_{c}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*$ = *r*/*σ*, spatial coordinate perpendicular to the interface $z*$ = *z*/*σ,* energy $u*$ = *u*/*ε*, force component *F*_{z}* = *F*_{z}/(*ε*/*σ*), temperature $T*$ = *k*_{B}*T*/*ε,* density *ρ** = *ρσ*^{3}, and surface tension *γ** = *γσ*^{2}/*ε*. 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^{8} 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*_{1} as

where *u*(r_{12}) is the intermolecular potential between particles at *r*_{1} = (*x*_{1}, *y*_{1}, *z*_{1}) and **r**_{2} = (*x*_{2}, *y*_{2}, *z*_{2}) with *r*_{12} = *r*_{2} − *r*_{1,}*r*_{12} = |*r*_{12}| and *z*_{12} = *z*_{2} − *z*_{1}. For the special case of the LJ potential, Eq. (2) yielded^{8} for the lrc Δ*u*(*z*_{1}) at a position *z*_{1},

Therefrom, these authors^{8} 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 Δ_{1}*F*_{z}(*z*_{1}) as

Several years later Janeček^{12} 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^{12} 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^{12} 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^{9} 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 *F*_{z}(*z*_{1}) at a position *z*_{1} 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 Δ_{1}*F*_{z}(*z*_{1}) given by Eq. (4). The second term Δ_{2}*F*_{z}(*z*_{1}) was given by Janeček^{12} 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 Δ_{2}*F*_{z}(*z*_{1}) which accounts for the lrc contribution at *z*_{1} from those particles which are outside the cut-off radius *r*_{c} and have a *z*-coordinate *z*_{2} with –*r*_{c} ≤ *z*_{12} ≤ *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 Δ_{1}*F*_{z}(*z*_{1}) from Ref. 8 and Δ*F*_{z}(*z*_{1}) 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 *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\sigma $**= n***dσ =* (*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*_{12}), then *G*(**x**, *t*) is *u*(*r*_{12})*ρ*(**r**_{2}) with *r*_{12} = |**r**_{2} − **r**_{1}|, **r**_{1} = (*x*_{1} = 0, *y*_{1} = 0, *z*_{1}) and **r**_{2} = (*x*_{2}, *y*_{2}, *z*_{2}), *t* = *z*_{1}, **x** = **r**_{2}, the domain *D*_{t} is given by *r*_{12} > *r*_{c}, and the velocity **v** is the unit vector in the z-direction **v** = **e**_{z}. Therewith **v**·$d\sigma $ = *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 −Δ_{1}*F*_{z}(*z*_{1}) from Ref. 8 given in the present Eq. (4), whilst the first integral on the rhs corresponds to the negative lrc −Δ*F*_{z}(*z*_{1}) 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*_{12}/*r*_{c}. Therewith one obtains

Finally, changing the integration over the surface area *dσ* to the integration over *z*_{2} by using *dσ* = 2π *r*_{c}*dz*_{2} and replacing *dz*_{2} by *dz*_{12} yields

which corresponds to the term Δ_{2}*F*_{z}(*z*_{1}) obtained by Janeček^{12} 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*_{1}) in integral form as

## III. SURFACE TENSION EQUATION FOR THE LENNARD-JONES FLUID

After the disagreement between Δ_{1}*F*_{z}(*z*_{1}) from Ref. 8 and Δ*F*_{z}(*z*_{1}) 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^{9} from which we take the results obtained with 2048 particles, the lrc Δ_{1}*F*_{z}(*z*_{1}), and r_{c} = 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 Δ*F*_{z}(*z*_{1}), 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.

T
. | γ_{S1}
. | γ_{S2}
. | Δγ_{S} (%)
. | γ_{C1}
. | γ_{C2}
. | γ_{C3}
. |
---|---|---|---|---|---|---|

0.70 | 1.1452 (194) | 1.150(4) | −0.42 | 1.1451 | 1.1515 | 1.1534 |

0.80 | 0.930(10) | 0.9187 | 0.9249 | 0.9248 | ||

0.85 | 0.8096(162) | 0.8096 | 0.8152 | 0.8144 | ||

0.90 | 0.707(8) | 0.7034 | 0.7082 | 0.7070 | ||

1.00 | 0.502(5) | 0.5012 | 0.5034 | 0.5021 | ||

1.10 | 0.3150 (122) | 0.310 (4) | 1.61 | 0.3150 | 0.3133 | 0.3129 |

1.20 | 0.144(8) | 0.1500 | 0.1434 | 0.1450 | ||

1.25 | 0.075(4) | 0.0787 | 0.0696 | 0.0726 |

T
. | γ_{S1}
. | γ_{S2}
. | Δγ_{S} (%)
. | γ_{C1}
. | γ_{C2}
. | γ_{C3}
. |
---|---|---|---|---|---|---|

0.70 | 1.1452 (194) | 1.150(4) | −0.42 | 1.1451 | 1.1515 | 1.1534 |

0.80 | 0.930(10) | 0.9187 | 0.9249 | 0.9248 | ||

0.85 | 0.8096(162) | 0.8096 | 0.8152 | 0.8144 | ||

0.90 | 0.707(8) | 0.7034 | 0.7082 | 0.7070 | ||

1.00 | 0.502(5) | 0.5012 | 0.5034 | 0.5021 | ||

1.10 | 0.3150 (122) | 0.310 (4) | 1.61 | 0.3150 | 0.3133 | 0.3129 |

1.20 | 0.144(8) | 0.1500 | 0.1434 | 0.1450 | ||

1.25 | 0.075(4) | 0.0787 | 0.0696 | 0.0726 |

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^{18}

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

## IV. CONCLUSIONS

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.

## ACKNOWLEDGMENTS

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.