JOURNAL OF CHEMICAL PHYSICS VOLUME 111, NUMBER 18 8 NOVEMBER 1999
Computer simulations of liquid/vapor interface in Lennard-Jones fluids:
Some questions and answers
Andrij Trokhymchuka)
Instituto de Quı́mica de la UNAM, Circuito Exterior, Coyoacán 04510, México D.F., Mexico
José Alejandreb)
Departamento de Quı́mica, Universidad Autónoma Metropolitana, Iztapalapa, Apdo. Postal 55-534, 09340,
México D.F., Mexico
共Received 23 February 1999; accepted 11 August 1999兲
Canonical molecular dynamics 共MD兲 and Monte Carlo 共MC兲 simulations for liquid/vapor
equilibrium in truncated Lennard-Jones fluid have been carried out. Different results for coexistence
properties 共orthobaric densities, normal and tangential pressure profiles, and surface tension兲 have
been reported in each method. These differences are attributed in literature to different set up
conditions, e.g., size of simulation cell, number of particles, cut-off radius, time of simulations, etc.,
applied by different authors. In the present study we show that observed disagreement between
simulation results is due to the fact that different authors inadvertently simulated different model
fluids. The origin of the problem lies in details of truncation procedure used in simulation studies.
Care has to be exercised in doing the comparison between both methods because in MC calculations
one deals with the truncated potential, while in MD calculations one uses the truncated forces, i.e.,
derivative of the potential. The truncated force does not uniquely define the primordial potential. It
results in MD and MC simulations being performed for different potential models. No differences
in the coexistence properties obtained from MD and MC simulations are found when the same
potential model is used. An additional force due to the discontinuity of the truncated potential at
cut-off distance becomes crucial for inhomogeneous fluids and has to be included into the virial
calculations in MC and MD, and into the computation of trajectories in MD simulations. The normal
pressure profile for the truncated potential is constant through the interface and both vapor and
liquid regions only when this contribution is taken into account, and ignoring it results in incorrect
value of surface tension. © 1999 American Institute of Physics. 关S0021-9606共99兲52441-0兴
I. INTRODUCTION that in MC calculations one deals directly with the pair po-
tential, while in MD calculations one uses the pair forces,
Computer simulations are one of the most powerful tools i.e., derivative of the potential. Truncation of the potential is
of the modern statistical mechanical theory of condensed performed in MC and truncation of the forces in MD algo-
matter. Generally, it is assumed that computer simulations rithm. MC creates configurations according to energy criteria
produce exact data for a given potential model. Two main while MD uses force route. The truncated force does not
factors which can affect the accuracy of simulation data are define uniquely the primordial potential. It results in the MD
caused by computational limitations, i.e., system size and and MC simulations being performed for different potential
truncation of interactions. Both of them have been discussed models: MC for spherically truncated 共ST兲 potential while
in literature and for the case of single-phase homogeneous MD for spherically truncated and shifted 共STS兲 potential.
systems, reasonable criteria have been established.1 In the The truncation of interactions has different consequences
case of two-phase systems, particularly when densities are dependent on the physical nature of the system under mod-
varying through the simulation cell, as in the case of liquid/ eling, i.e., whether it is simple 共nonpolar兲 or complex 共polar,
liquid or liquid/vapor phase coexistence, it becomes more ionic, etc.兲 fluids. Electrostatic interactions 共Coulombic, di-
complicated.2 polar, and higher multipole兲 determine the peculiar proper-
Two methodologies are usually used in direct computer ties of a such systems 共conductivity, dielectric permittivity,
simulations 共coexisting phases are in physical contact and etc.兲 caused by charge and polarization fluctuations. Nonad-
interface region is presented兲 of phase coexistence in equate 共not large enough兲 truncation of full interactions
fluids:3–43 conventional canonical Monte Carlo 共MC兲 and might change the physics of the original system. Examples
molecular dynamics 共MD兲. Chapela et al.8 have pointed out are water and aqueous solutions where application of ad-
equate truncation procedure for electrostatic forces in inho-
a兲
Permanent address: Institute for Condensed Matter Physics, National mogeneous simulations becomes crucial to preserve a physi-
Academy of Sciences of the Ukraine, Lviv 11, Ukraine. cally correct microscopic model, especially when
b兲
Also at Departmento de Simulación Molecular, Instituto Mexicana del
Petróleo, Eje central Lazaro C’ardenas 152, Apdo. Postal 14-805, 07730, electrostatic information is to be obtained. In particular, it
Mexico D.F., Mexico. has been shown by Spohr44 that the use of truncated interac-
0021-9606/99/111(18)/8510/14/$15.00 8510 © 1999 American Institute of Physics
Downloaded 02 Jun 2006 to 131.151.76.4. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
, J. Chem. Phys., Vol. 111, No. 18, 8 November 1999 Liquid/vapor interface in Lennard-Jones fluids 8511
tions might lead to unphysical results for electrostatic poten- different potential models have not been distinguished, and
tials. In contrast, for simple fluids, where the long-range in- 共ii兲 the impulsive contribution from the derivative of the
teractions between molecules are dominated by dispersion truncated potential at cut-off distance to the pressure and
forces, truncation does not change the physics of original forces has not been taken into account. As a consequence,
system, but slightly modifies original model due to the the results obtained for liquid/vapor coexistence in LJ fluid
changes in potential energy. Additionally, for the spherically are confusing and some questionable conclusions have been
truncated potential, an ‘‘impulsive’’ contribution to the pres- drawn.
sure, due to the discontinuous change of the potential at cut- The differences between ST and STS models in simula-
off distance, appears. In the modeling of single-phase simple tions of liquid/vapor coexistence in complex fluids,28–43
fluids the changes in original model, introduced by trunca- where LJ interactions are part of total Hamiltonian, have not
tion procedure, have been discussed in literature 共see, e.g., been addressed as well. Results of prior works on nonpolar
Refs. 1,2兲 and found as not of great importance for the ma- molecular systems show that LJ-type contribution to the co-
jority of bulk properties. Is it true for a two-phase system? existence properties might be significant,32–34,36,39 and there-
The review of the literature sources on the subject of the fore a correct handling of truncated LJ interactions is impor-
phase coexistence in simple fluids seems to indicate that it tant. In the simulations of liquid/vapor coexistence in highly
may not be valid any longer. polar fluids, like water, the contribution of electrostatic
To obtain a conclusive answer on the above question is forces prevails31 though care has to be taken since besides
indeed the main objective of the present article. To accom- the LJ interactions, there is contribution of so-called ‘‘real’’
plish this objective we have chosen a classical pure atomic part of electrostatic interactions separated from the full elec-
Lennard-Jones 共LJ兲 12–6 fluid. This model itself is of im- trostatic interactions if Ewald summation technique is ap-
portance and is by far the best studied continuous potential plied. The real part has a complementary function form and
because it provides moderately well a description of liquid/ its contribution depends on the parameters employed in
vapor coexistence properties of nonpolar real fluids of Ewald scheme.
spherical 共argon, krypton, xenon兲 and homonuclear diatomic Additional questions arise regarding simulation of
共oxygen, nitrogen兲 molecules.3–22,45–47 LJ model also is liquid/vapor interface in LJ fluids. Is the LJ potential an ad-
widely employed in modeling of phase coexistence in mix- equate model to describe the liquid/vapor coexistence and
tures of simple fluids23–27 as well as for more complex sys- interfacial properties in a real nonpolar fluid?4,6,46,47 To an-
tems, like liquid/vapor coexistence in nonpolar 共chlorine, swer this question the simulations of the untruncated LJ po-
hexane, alkanes兲 and polar 共water, methanol兲 molecular tential are required. For coexisting densities and pressure,
fluids,28–36 liquid/liquid interfaces in aqueous solutions 共eth- simulations of full LJ potential have been done using the
anol, octane, hydrocarbon, benzene, dichloroethane兲,37–41 bi- GEMC method.50,51 For interfacial properties 共e.g., surface
layers and monolayers in water,42,43 etc., where the disper- tension兲 it can be done either in the way similar to the single-
sion or atom–atom interactions are described by LJ potential. phase case, applying the long-range corrections 共LRC兲 to the
Therefore, it becomes important to have correct and unam- ‘‘reference’’ results obtained in direct canonical simulations
biguous computer modeling of phase coexistence in LJ fluid. with truncated potential,5,8,15,21 or including the LRC scheme
Despite a number of articles on liquid/vapor coexistence explicitly in the direct canonical simulations. The first way is
in LJ fluid published in the past, the situation is far from clearer and simpler but results obtained depend on the prop-
clear. In particular, Smit and Frenkel48,49 using Gibbs en- erties of the reference system. In particular, it has been
semble MC 共GEMC兲 simulations 共in the context of the prob- shown by Blokhuis et al.21 that, corrected in such manner,
lems discussed in the present study, GEMC simulations can surface tension is different for reference calculations per-
be considered as homogeneous or bulk, since coexisting formed with different cut-off radius. To apply LRC explicitly
phases are not in physical contact and no interface is pre- in inhomogeneous simulations is not an easy task. Different
sented兲 have shown that in the case of LJ fluid different schemes to include LRC during the direct canonical simula-
temperature–density coexisting envelopes are produced de- tions of the interface have been proposed recently17,18 and
pending on how the potential is truncated, i.e., whether it is distinct results have been obtained. Guo and Lu17 have simu-
ST or STS. In direct canonical simulations of phase equilib- lated the liquid/vapor interface in LJ fluid using canonical
ria, where the phases coexist through the interface, only the MC and performing LRC of the potential energy, pressure
effect of the value of cut-off distance on the properties of tensor, and surface tension. These authors have concluded
liquid/vapor interface has been analyzed so far;12,15,17–19 the that without taking LRC into account, constant normal pres-
way the potential is truncated 共ST or STS兲 has not been taken sure profile through the interface and the vapor and liquid
into account. However, different truncation procedures bulk regions 共that is a necessary condition of the mechanical
change the phase diagram and can by no means be ignored equilibrium in canonical simulations兲 could not be reached in
when coexisting densities and properties of interfacial region simulations with a cut-off radius smaller than 3.1. We dis-
are to be obtained, and when different studies are to be com- agree with this conclusion. The MD simulations results re-
pared. Nevertheless, this issue continues not to be suffi- ported in the literature indicate 共see, e.g., Refs. 7,10兲 that
ciently addressed by many authors 共see, e.g., Refs. 12,15– normal pressure profile is constant for a cut-off radius as
17兲 and a high scattering of the results on coexisting small as 2.5.
densities and surface tension in LJ fluid may be found be- It is our primary interest to clarify these subtle points in
cause 共i兲 applications of canonical MC or MD simulations to the comparison of canonical MD and MC simulations of
Downloaded 02 Jun 2006 to 131.151.76.4. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
Computer simulations of liquid/vapor interface in Lennard-Jones fluids:
Some questions and answers
Andrij Trokhymchuka)
Instituto de Quı́mica de la UNAM, Circuito Exterior, Coyoacán 04510, México D.F., Mexico
José Alejandreb)
Departamento de Quı́mica, Universidad Autónoma Metropolitana, Iztapalapa, Apdo. Postal 55-534, 09340,
México D.F., Mexico
共Received 23 February 1999; accepted 11 August 1999兲
Canonical molecular dynamics 共MD兲 and Monte Carlo 共MC兲 simulations for liquid/vapor
equilibrium in truncated Lennard-Jones fluid have been carried out. Different results for coexistence
properties 共orthobaric densities, normal and tangential pressure profiles, and surface tension兲 have
been reported in each method. These differences are attributed in literature to different set up
conditions, e.g., size of simulation cell, number of particles, cut-off radius, time of simulations, etc.,
applied by different authors. In the present study we show that observed disagreement between
simulation results is due to the fact that different authors inadvertently simulated different model
fluids. The origin of the problem lies in details of truncation procedure used in simulation studies.
Care has to be exercised in doing the comparison between both methods because in MC calculations
one deals with the truncated potential, while in MD calculations one uses the truncated forces, i.e.,
derivative of the potential. The truncated force does not uniquely define the primordial potential. It
results in MD and MC simulations being performed for different potential models. No differences
in the coexistence properties obtained from MD and MC simulations are found when the same
potential model is used. An additional force due to the discontinuity of the truncated potential at
cut-off distance becomes crucial for inhomogeneous fluids and has to be included into the virial
calculations in MC and MD, and into the computation of trajectories in MD simulations. The normal
pressure profile for the truncated potential is constant through the interface and both vapor and
liquid regions only when this contribution is taken into account, and ignoring it results in incorrect
value of surface tension. © 1999 American Institute of Physics. 关S0021-9606共99兲52441-0兴
I. INTRODUCTION that in MC calculations one deals directly with the pair po-
tential, while in MD calculations one uses the pair forces,
Computer simulations are one of the most powerful tools i.e., derivative of the potential. Truncation of the potential is
of the modern statistical mechanical theory of condensed performed in MC and truncation of the forces in MD algo-
matter. Generally, it is assumed that computer simulations rithm. MC creates configurations according to energy criteria
produce exact data for a given potential model. Two main while MD uses force route. The truncated force does not
factors which can affect the accuracy of simulation data are define uniquely the primordial potential. It results in the MD
caused by computational limitations, i.e., system size and and MC simulations being performed for different potential
truncation of interactions. Both of them have been discussed models: MC for spherically truncated 共ST兲 potential while
in literature and for the case of single-phase homogeneous MD for spherically truncated and shifted 共STS兲 potential.
systems, reasonable criteria have been established.1 In the The truncation of interactions has different consequences
case of two-phase systems, particularly when densities are dependent on the physical nature of the system under mod-
varying through the simulation cell, as in the case of liquid/ eling, i.e., whether it is simple 共nonpolar兲 or complex 共polar,
liquid or liquid/vapor phase coexistence, it becomes more ionic, etc.兲 fluids. Electrostatic interactions 共Coulombic, di-
complicated.2 polar, and higher multipole兲 determine the peculiar proper-
Two methodologies are usually used in direct computer ties of a such systems 共conductivity, dielectric permittivity,
simulations 共coexisting phases are in physical contact and etc.兲 caused by charge and polarization fluctuations. Nonad-
interface region is presented兲 of phase coexistence in equate 共not large enough兲 truncation of full interactions
fluids:3–43 conventional canonical Monte Carlo 共MC兲 and might change the physics of the original system. Examples
molecular dynamics 共MD兲. Chapela et al.8 have pointed out are water and aqueous solutions where application of ad-
equate truncation procedure for electrostatic forces in inho-
a兲
Permanent address: Institute for Condensed Matter Physics, National mogeneous simulations becomes crucial to preserve a physi-
Academy of Sciences of the Ukraine, Lviv 11, Ukraine. cally correct microscopic model, especially when
b兲
Also at Departmento de Simulación Molecular, Instituto Mexicana del
Petróleo, Eje central Lazaro C’ardenas 152, Apdo. Postal 14-805, 07730, electrostatic information is to be obtained. In particular, it
Mexico D.F., Mexico. has been shown by Spohr44 that the use of truncated interac-
0021-9606/99/111(18)/8510/14/$15.00 8510 © 1999 American Institute of Physics
Downloaded 02 Jun 2006 to 131.151.76.4. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
, J. Chem. Phys., Vol. 111, No. 18, 8 November 1999 Liquid/vapor interface in Lennard-Jones fluids 8511
tions might lead to unphysical results for electrostatic poten- different potential models have not been distinguished, and
tials. In contrast, for simple fluids, where the long-range in- 共ii兲 the impulsive contribution from the derivative of the
teractions between molecules are dominated by dispersion truncated potential at cut-off distance to the pressure and
forces, truncation does not change the physics of original forces has not been taken into account. As a consequence,
system, but slightly modifies original model due to the the results obtained for liquid/vapor coexistence in LJ fluid
changes in potential energy. Additionally, for the spherically are confusing and some questionable conclusions have been
truncated potential, an ‘‘impulsive’’ contribution to the pres- drawn.
sure, due to the discontinuous change of the potential at cut- The differences between ST and STS models in simula-
off distance, appears. In the modeling of single-phase simple tions of liquid/vapor coexistence in complex fluids,28–43
fluids the changes in original model, introduced by trunca- where LJ interactions are part of total Hamiltonian, have not
tion procedure, have been discussed in literature 共see, e.g., been addressed as well. Results of prior works on nonpolar
Refs. 1,2兲 and found as not of great importance for the ma- molecular systems show that LJ-type contribution to the co-
jority of bulk properties. Is it true for a two-phase system? existence properties might be significant,32–34,36,39 and there-
The review of the literature sources on the subject of the fore a correct handling of truncated LJ interactions is impor-
phase coexistence in simple fluids seems to indicate that it tant. In the simulations of liquid/vapor coexistence in highly
may not be valid any longer. polar fluids, like water, the contribution of electrostatic
To obtain a conclusive answer on the above question is forces prevails31 though care has to be taken since besides
indeed the main objective of the present article. To accom- the LJ interactions, there is contribution of so-called ‘‘real’’
plish this objective we have chosen a classical pure atomic part of electrostatic interactions separated from the full elec-
Lennard-Jones 共LJ兲 12–6 fluid. This model itself is of im- trostatic interactions if Ewald summation technique is ap-
portance and is by far the best studied continuous potential plied. The real part has a complementary function form and
because it provides moderately well a description of liquid/ its contribution depends on the parameters employed in
vapor coexistence properties of nonpolar real fluids of Ewald scheme.
spherical 共argon, krypton, xenon兲 and homonuclear diatomic Additional questions arise regarding simulation of
共oxygen, nitrogen兲 molecules.3–22,45–47 LJ model also is liquid/vapor interface in LJ fluids. Is the LJ potential an ad-
widely employed in modeling of phase coexistence in mix- equate model to describe the liquid/vapor coexistence and
tures of simple fluids23–27 as well as for more complex sys- interfacial properties in a real nonpolar fluid?4,6,46,47 To an-
tems, like liquid/vapor coexistence in nonpolar 共chlorine, swer this question the simulations of the untruncated LJ po-
hexane, alkanes兲 and polar 共water, methanol兲 molecular tential are required. For coexisting densities and pressure,
fluids,28–36 liquid/liquid interfaces in aqueous solutions 共eth- simulations of full LJ potential have been done using the
anol, octane, hydrocarbon, benzene, dichloroethane兲,37–41 bi- GEMC method.50,51 For interfacial properties 共e.g., surface
layers and monolayers in water,42,43 etc., where the disper- tension兲 it can be done either in the way similar to the single-
sion or atom–atom interactions are described by LJ potential. phase case, applying the long-range corrections 共LRC兲 to the
Therefore, it becomes important to have correct and unam- ‘‘reference’’ results obtained in direct canonical simulations
biguous computer modeling of phase coexistence in LJ fluid. with truncated potential,5,8,15,21 or including the LRC scheme
Despite a number of articles on liquid/vapor coexistence explicitly in the direct canonical simulations. The first way is
in LJ fluid published in the past, the situation is far from clearer and simpler but results obtained depend on the prop-
clear. In particular, Smit and Frenkel48,49 using Gibbs en- erties of the reference system. In particular, it has been
semble MC 共GEMC兲 simulations 共in the context of the prob- shown by Blokhuis et al.21 that, corrected in such manner,
lems discussed in the present study, GEMC simulations can surface tension is different for reference calculations per-
be considered as homogeneous or bulk, since coexisting formed with different cut-off radius. To apply LRC explicitly
phases are not in physical contact and no interface is pre- in inhomogeneous simulations is not an easy task. Different
sented兲 have shown that in the case of LJ fluid different schemes to include LRC during the direct canonical simula-
temperature–density coexisting envelopes are produced de- tions of the interface have been proposed recently17,18 and
pending on how the potential is truncated, i.e., whether it is distinct results have been obtained. Guo and Lu17 have simu-
ST or STS. In direct canonical simulations of phase equilib- lated the liquid/vapor interface in LJ fluid using canonical
ria, where the phases coexist through the interface, only the MC and performing LRC of the potential energy, pressure
effect of the value of cut-off distance on the properties of tensor, and surface tension. These authors have concluded
liquid/vapor interface has been analyzed so far;12,15,17–19 the that without taking LRC into account, constant normal pres-
way the potential is truncated 共ST or STS兲 has not been taken sure profile through the interface and the vapor and liquid
into account. However, different truncation procedures bulk regions 共that is a necessary condition of the mechanical
change the phase diagram and can by no means be ignored equilibrium in canonical simulations兲 could not be reached in
when coexisting densities and properties of interfacial region simulations with a cut-off radius smaller than 3.1. We dis-
are to be obtained, and when different studies are to be com- agree with this conclusion. The MD simulations results re-
pared. Nevertheless, this issue continues not to be suffi- ported in the literature indicate 共see, e.g., Refs. 7,10兲 that
ciently addressed by many authors 共see, e.g., Refs. 12,15– normal pressure profile is constant for a cut-off radius as
17兲 and a high scattering of the results on coexisting small as 2.5.
densities and surface tension in LJ fluid may be found be- It is our primary interest to clarify these subtle points in
cause 共i兲 applications of canonical MC or MD simulations to the comparison of canonical MD and MC simulations of
Downloaded 02 Jun 2006 to 131.151.76.4. Redistribution subject to AIP license or copyright, see http://jcp.aip.org/jcp/copyright.jsp
