RESEARCH ARTICLE | FEBRUARY 01 1995

Computer simulations of vapor–liquid phase equilibria of n‐

alkanes 
Berend Smit; Sami Karaborni; J. Ilja Siepmann

J. Chem. Phys. 102, 2126–2140 (1995)

https://doi.org/10.1063/1.469563

CrossMark

 
View Export
Online Citation

Articles You May Be Interested In

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

Detection of chemical compounds and its antioxidant activity of aniseeds extract
AIP Conference Proceedings (October 2022)

Uncertainty quantification and propagation of errors of the Lennard-Jones 12-6 parameters for n-alkanes
J. Chem. Phys. (May 2017)

On the use of F0 variations in automatic speech recognition

J Acoust Soc Am (August 2005)
Computer simulations of vapor–liquid phase equilibria of n -alkanes
Berend Smit,a) Sami Karaborni, and J. Ilja Siepmannb)
Shell Research B.V., Koninklijke/Shell-Laboratorium, Amsterdam, P.O. Box 38000, 1030 BN Amsterdam, The
Netherlands
~Received 1 April 1994; accepted 25 October 1994!
For petrochemical applications knowledge of the critical properties of the n-alkanes is of interest
even at temperatures where these molecules are thermally unstable. Computer simulations can
determine the vapor–liquid coexistence curve of a large number of n-alkanes ranging from pentane
~C5! through octatetracontane ~C48!. We have compared the predicted phase diagrams of various
models with experimental data. Models which give nearly identical properties of liquid alkanes at
standard conditions may have critical temperatures that differ by more than 100 K. A new n-alkane
model has been developed by us that gives a good description of the phase behavior over a large
temperature range. For modeling vapor–liquid coexistence a relatively simple united atom model
was sufficient to obtain a very good agreement with experimental data; thus it appears not necessary
to take the hydrogen atoms explicitly into account. The model developed in this work has been used
to determine the critical properties of the long-chain alkanes for which experiments turned out to be
difficult and contradictory. We found that for the long-chain alkanes ~C8 –C48! the critical density
decreases as a function of the carbon number. These simulations were made possible by the use of

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

a recently developed simulation technique, which is a combination of the Gibbs-ensemble technique
and the configurational-bias Monte Carlo method. Compared with the conventional Gibbs-ensemble
technique, this method is several orders of magnitude more efficient for pentane and up to a hundred
orders of magnitude for octatetracontane. This recent development makes it possible to perform
routinely phase equilibrium calculations of complex molecules. © 1995 American Institute of
Physics.

I. INTRODUCTION One of the Monte Carlo steps in the Gibbs-ensemble

Anyone who has tried to construct an igloo in the desert technique is the transfer of molecules between the liquid
can testify to the importance of elementary knowledge of phase and gas phase. For chain molecules, this step results in
phase behavior. Besides being of practical interest, phase a prohibitively low acceptance of transfers from the gas
equilibria have been the topic of many fundamental studies phase into the liquid phase. Therefore the Gibbs ensemble
since the seminal work of van der Waals. Since knowledge of used to be limited to systems containing atoms or small mol-
the phase behavior is essential in many practical applica- ecules. Recently, the Gibbs-ensemble technique has been
tions, there have been significant experimental efforts to- combined with the configurational-bias Monte Carlo
wards the determination of unknown or partially known method.6 – 8 Instead of a random insertion, in the
phase diagrams. configurational-bias Monte Carlo scheme molecules are
It is interesting to note that in the work of van der Waals grown atom by atom in such a way that regions of favorable
a connection has already been made between the intermo- energy are found and overlap with other molecules is
lecular potential and the phase behavior of the molecules. To avoided.9–12 This growing scheme introduces a bias that can
determine the phase diagram of a given model proved to be be removed exactly by adjusting the acceptance rules.9,12 A
an extremely difficult task. Exact analytical solutions have similar approach has been used by Cracknell et al.13 to per-
been obtained for only a few important but exceptional form Gibbs-ensemble simulations using a rotational-bias in-
cases. For exact data for a given model one therefore has to sertion of water molecules. The combination of the Gibbs-
rely on the results of computer simulations. The calculation ensemble technique with the configurational-bias Monte
of a phase diagram via computer simulations used to be an Carlo method has been applied successfully to determine the
elaborate task which required many simulations.1 An impor- vapor–liquid coexistence curve of chains of Lennard-Jones
tant step forward was the development of the Gibbs- beads6 and alkanes.8,14,15
ensemble technique by Panagiotopoulos.2– 4 By this method The description of alkanes has received considerable in-
data on phase coexistence can be obtained from one single terest. Many different models for these molecules have been
simulation. The Gibbs-ensemble technique has been applied proposed.16 –25,11 One of the first models for liquid butane
successfully to determine the vapor–liquid and liquid–liquid was developed by Ryckaert and Bellemans.17 This model
coexistence curves of various model fluids. Panagiotopoulos5 assumes that every CH3 or CH2 group can be described as
recently reviewed the applications of the Gibbs-ensemble one single interaction site. The dispersive interactions of
technique. these ‘‘united atoms’’ are described by a Lennard-Jones po-
tential. The bond lengths and bond angles are kept fixed.
a! Changes in the dihedral angle are described with a torsion
Author for correspondence.
b!
Present address: Department of Chemistry, University of Minnesota, potential which has been fitted to yield the experimental dis-
Minneapolis, Minnesota 55455. tribution of gauche/trans conformers. This type of model

2126 J. Chem. Phys. 102 (5), 1 February 1995 0021-9606/95/102(5)/2126/15/$6.00 © 1995 American Institute of Physics
 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria 2127

was extended to longer chain lengths by Jorgensen and N ~ b ,V 1 ,n 1 ; j N !

co-workers.19,20 At high densities the assumption that the n
CH3 and CH2 groups can be modeled as united atoms tends V 1 1 ~ V2V 1 ! N2n 1
to fail. Indeed simulations of alkane crystals26,27 and dense }
~ N2n 1 ! !n 1 !
monolayers of alkanes28,29 show that united atom models do
not describe the packing of alkane molecules correctly, and 3exp@ 2 b U 1 ~ n 1 !# exp@ 2 b U 2 ~ N2n 1 !# , ~1!
that it is important to take the hydrogen atoms explicitly into where n 1 denotes the number of particles in box 1, V 1 the
account. Because of the increased number of interaction volume of box 1, j N denotes the scaled ~with respect to the
sites, these ‘‘all atoms models’’16 are much more demanding box length! positions of the particles, and U(n i ) is the inter-
with regard to their use in a simulation. An interesting com- molecular potential.
promise between the united atom approach and the all atoms In a Gibbs-ensemble simulations the following Monte
models is the use of anisotropic potentials as proposed by Carlo moves are used: displacement of particles in the boxes,
Toxvaerd.21 The anisotropy is introduced to make the inter- changes in the volume, and exchange of particles between
action between the CH3 and CH2 groups dependent on the the two boxes. In the next section we use the configurational-
conformation, i.e., the interaction is different when the CH2 bias Monte Carlo method for the exchange of particles. To
group points with its ‘‘H-side’’ or ‘‘C-side’’ towards another introduce the notation, we consider the acceptance rules for
CH2 group. Toxvaerd has shown that such an anisotropic this step in some detail. The derivations of the acceptance
model gives a better description of the equation of state of rules for the other moves are given elsewhere.3,30
dense alkanes under high pressure than that of some of the Let us assume the system to be in a state o with n 1

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

isotropic models.21 particles in box 1 with volume V 1 and consider the move to
Most of these alkane models have been fitted to liquid the state n which has n 1 11 particles in box 1 with the same
properties such as heats of vaporization and liquid densities volume. The acceptance rule for this move is3

F
at standard conditions. In this work, we address the question
how accurately the phase behavior of n-alkanes can be mod- V 1 ~ N2n 1 !
acc~ o→n ! 5min 1,
eled over a large range of temperatures and chain lengths. ~ V2V 1 !~ n 1 11 !

G
The critical properties of the n-alkanes are of interest for
petrochemical applications. We use the Gibbs-ensemble 3exp ~ 2 b DU 1 ! exp ~ 2 b DU 2 ! , ~2!
simulations to estimate the chain length dependence of these
critical properties at conditions where experiments are not where DU 1 5U 1 (n)2U 1 (o) is defined as the energy differ-
possible. ence in box 1 between state n and state o. This acceptance
In Secs. II and III we describe the simulation techniques rule can be derived by imposing the condition of detailed
and in Sec. IV the models that have been studied in this balance
work, as well as the results of the simulations. In Sec. V, we
present the simulation results of the critical properties. Some K ~ o→n ! 5K ~ n→o ! , ~3!
preliminary results of this work have been described where K(o→n) is the flow of configurations from o to n.
earlier.14,15 This flow of configurations is equal to the product of the
probability of being in state o, the probability of generating
state n and the probability of acceptance
K ~ o→n ! 5N ~ n 1 ,V 1 ; j N ! 3 p ~ o→n ! 3acc~ o→n ! . ~4!
II. SIMULATION TECHNIQUES
For the reverse move, the removal of a particles from box 1,
the flow is given by
In this section a description is given of the simulation
techniques that are used in this work. A more extensive de- K ~ n→o ! 5N ~ n 1 11,V 1 ; j N ! 3 p ~ n→o ! 3acc~ n→o ! .
scription of the Gibbs-ensemble technique4,30 and the ~5!
configurational-bias Monte Carlo method11,12,31 can be found Since it is decided at random whether to remove or insert a
in the literature. particle, we have p(o→n)5 p(n→o). Substitution of Eqs.
A. The Gibbs ensemble technique ~4! and ~5! and distribution ~1! together with the acceptance
rule ~2! shows that indeed detailed balance ~3! is obeyed.
Simulations in the Gibbs ensemble are performed using
two boxes, each box having periodic boundary conditions.
B. Configurational-bias Monte Carlo
The boxes are kept at a constant temperature T, the total
volume of the two boxes is fixed at V, and a fixed number of In the conventional Gibbs-ensemble scheme particles are
N particles are distributed over the two boxes. The two boxes inserted at random positions. For a Lennard-Jones fluid the
are coupled via Monte Carlo rules that allow the exchange of probability that an attempt to insert a particle in the liquid
particles and changes in the volume in such a way that the phase does not result in an overlap with one of the other
two boxes remain in thermodynamic equilibrium with each liquid particles is of the order of 0.005.33 At similar condi-
other. tions, the probability that a chain of n atoms is successfully
The probability of finding a particular configuration in inserted will be of the order of 0.005 n . As a consequence the
the Gibbs ensemble is given by4,5,32 number of attempts to insert a particle increases enormously

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995

2128 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria

w l~ n ! 5 ( exp@ 2 b u ext
l ~ b j !# . ~8!
j51

Out of these k trial positions, we select one with prob-

ability
exp@ 2 b u ext
l ~ bi !#
l ~ bi ! 5
p ext . ~9!
w l~ n !
~3! Step 2 is repeated M 21 times until the entire molecule
is grown and the Rosenbluth factor of the molecule can
be calculated
M

W~ n !5 ) w l~ n ! . ~10!
l51

This algorithm biases the insertion of a molecule such

that regions with favorable energy are found and overlap
FIG. 1. Insertion of a chain molecule; the arrows indicate the k trial orien- with other atoms is avoided. The probability that a par-

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

tations to insert the fourth atom. ticular conformation is generated is given by
M

P~ n !5 ) p int ext
l ~n!pl ~n!
l52
for larger chain lengths. This limits the applicability of the
Gibbs-ensemble technique in its original form to very short M
exp$ 2 b @ u int
l ~ n ! 1u l ~ n !# %
ext
chain molecules.
The configurational-bias Monte Carlo technique has
5 ) Cw l ~ n !
l52
been developed to insert chain molecules in moderately
dense systems.10–12 Here we give a brief description of this exp@ 2 b U ~ n !#
5 , ~11!
method, a more extensive discussion is given else C M 21 W ~ n !
where.11,12,31
Let us divide the potential energy of an atom into two where the total energy of the inserted molecule is
contributions;12,34 ~1! the internal energy u int which includes M M
parts of the intramolecular interactions, and ~2! the external U5 ( u l5 ( u intl u ext
l . ~12!
energy u ext which contains the intermolecular interactions l51 l51
and those intramolecular interactions that are not part of the
internal energy. The division is to some extent arbitrary and To perform a move, we have to calculate the Rosenbluth
depends on the details of the model. Note that in some factor of the old configuration. This is done via the following
implementations of the configurational-bias Monte Carlo steps:
technique this division is not used.11 In Sec. IV B 2 we make ~1! A particle is selected at random;
a detailed comparison of the advantages of using this sepa- ~2! The energy of the first atom is determined u 1 (o) to-
ration. gether with
Instead of a random insertion of a molecule, we use the
following procedure to ‘‘grow’’ a molecule atom by atom; w 1 ~ o ! 5exp@ 2 b u 1 ~ o !# . ~13!
~1! The first atom is inserted at a random position, and the ~3! For the next atom, l, k21 trial orientations are generated
energy u 1 (n) is calculated together with with a probability given by Eq. ~7!. These trial orienta-
tions together with the actual position of the atom l form
w 1 ~ n ! 5exp@ 2 b u 1 ~ n !# . ~6! the set $b8%k ~see Fig. 2! for which we determine the
factor
~2! To insert the next atom l, k trial orientations are gener-
ated ~see Fig. 1!. The set of k trial orientations are de- k

noted by $b%k 5b1 , b2 ,...,bk . These orientations are not w l~ o ! 5 ( exp@ 2 b u l ~ b8j !# . ~14!
generated at random, but with a probability which is a j51
function of the internal energy
~4! Step 2 is repeated M 21 times until we have retraced the
exp@ 2 b u int entire chain and the Rosenbluth factor of the chain can
l ~ bi !#
l ~ bi ! 5
p int . ~7! be calculated
C
M
Of each of these trial orientations the external energy is W~ o !5 ) w l~ o ! . ~15!
calculated @ u ext
l (bi ) # together with the factor l51

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995

 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria 2129

For the reverse move, the insertion of a chain in box 2, we

have
exp@ 2 b U ~ o !#
P ~ n→o u $ b8 % k ! 5 . ~18!
C M 21 W ~ o u $ b8 % k !
Note that there are many different ways to generate a par-
ticular configuration n or o, namely all sets of trial orienta-
tions that include the selected orientation. Detailed balance
implies that we have to sum over all of these. We can, how-
ever, impose a much stronger condition, super detailed bal-
ance, which states that for all sets of trial conformations
individually detailed balance should be obeyed.12 If super-
detailed balance is obeyed, then detailed balance is certainly
obeyed. By definition, for super-detailed balance we have to
consider the same set of trial orientations for the moves
o→n and n→o, so
N ~ o ! 3 p ~ o→n u $ b% k , $ b8 % k ! 3acc~ o→n u $ b% k , $ b8 % k !
FIG. 2. Calculation of the Rosenbluth factor of the old conformation; the
arrows give the set of directions for which the Rosenbluth factor of the 5N ~ n ! 3 p ~ n→o u $ b8 % k , $ b% k !

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

second atom is calculated.
3acc~ n→o u $ b8 % k , $ b% k ! . ~19!
Substitution of Eqs. ~1!, ~17!, and ~18! gives as condition for
The above algorithms for the new and old configuration the acceptance rule
form the basis of the configurational-bias Monte Carlo tech-
nique. They need to be supplemented with acceptance rules acc~ o→n u $ b% k , $ b8 % k ! V 1 ~ N2n 1 ! W~ n !
5 . ~20!
that remove the bias from the insertion step. These accep- acc~ n→o u $ b8 % k , $ b% k ! ~ V2V 1 !~ n 1 11 ! W ~ o !
tance rules depend on the type of move and type of en- Since acceptance rule ~16! obeys this condition, we have
semble. For example, in Refs. 10–12 acceptance rules are demonstrated that the correct distribution is sampled.
derived for a move in which part of the molecule is regrown. We have outlined the general scheme of the Gibbs-
ensemble technique combined with configurational-bias
Monte Carlo. For a given model it is important to tune the
C. Configurational-bias Monte Carlo and the Gibbs technique optimally as will be discussed in the next sections.
ensemble
In the Gibbs ensemble, we use the configurational-bias
III. DESCRIPTION OF THE SIMULATIONS
Monte Carlo technique to make the exchange of chain mol-
ecules between the two boxes possible. The simulations have been performed in cycles. Each
Let us assume the system to be in state o with n 1 par- cycle consists of R randomly selected Monte Carlo moves
ticles in box 1 with volume V 1 and we try to generate state n ~we usually take R equal to the total number of molecules!.35
by moving a particle from box 2 into box 1. We use the The type of moves we perform are ~1! displacement of
algorithms of the previous section to grow a chain in box 1 a randomly selected particle; ~2! rotation of a randomly
and to calculate the Rosenbluth factor of the old conforma- selected particle around the middle atom; ~3! regrowing
tion of the chain which is removed from box 2. We then of parts of a randomly selected molecule; ~4! change of
accept this move with probability volume of the two boxes; and ~5! exchange of particles be-

F G
tween the two boxes. The relative probability that a particu-
V 1 ~ N2n 1 ! W~ n !
acc~ o→n ! 5min 1, . ~16! lar move is attempted is set to p 1 :p 2 : p 3 : p 4 : p 5
~ V2V 1 !~ n 1 11 ! W ~ o ! 50.222:0.222:0.222:0.006:0.328. At low temperatures
We now have to demonstrate that this acceptance rule and for long chain alkanes the relative probability of attempt-
indeed removes the bias from the insertion step and hence ing an exchange of particles was increased to ensure a suffi-
the method indeed samples the correct distribution of con- cient number of successful exchanges.
figurations. As in Sec. II A, we impose the condition of de- Note that moves ~1! and ~2! do not change the internal
tailed balance ~3!. The main difference is that in the structure of the molecule. In these moves the maximum dis-
configurational-bias Monte Carlo scheme the probability of placement and maximum rotation are adjusted in such a way
generating a particular conformation does depend on the par- that 50% of the moves are accepted.
ticular configuration of the molecules and the probability of For move ~3!, the partial regrowing of a molecule, we
generating the reverse move will be different. The probabil- select a molecule at random and choose the number of atoms
ity of generating conformation n is given by @see Eq. ~11!#, that are to be regrown. With equal probability we regrow the
atoms at the end or beginning of that part of the molecule
exp@ 2 b U ~ n !# that does not get regrown. We use the configurational-bias
P ~ o→n u $ b% k ! 5 . ~17!
C M 21 W ~ n u $ b% k ! Monte Carlo technique for this move with acceptance rules

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995

2130 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria

as given in Ref. 12. The number of trial orientations range TABLE I. Comparison of nonbonded interaction parameters used to model
from six for C5 to ten for C48 . The total number of molecules alkanes. The energy parameter can be converted to ~kJ/mol! by multiplica-
was 200 for the short chain alkanes and 100 for the long tion by 0.008 315.
chain alkanes. e CH3 e CH2 s CH3 s CH2
During a volume change ~4!, we rescale the coordinates ~K! ~K! Å Å Ref.
keeping the internal conformation of the molecule fixed. The
C5 –C15 90.5 49.3 3.94 3.94 8,25
maximum volume change is set such that 50% of the moves
C24 and C71 49.3 49.3 3.94 3.94 25
are accepted. For the exchange step ~5!, we have used the C5 –C8 104.0 49.7 3.923 3.923 42,43
algorithm as described in the previous section. The number C16 50.5 50.5 4.045 4.045 79
of trial orientations was equal to the number used for partial C5 –C8 102.0 51.3 3.983 3.863 42,43
regrowing of a molecule. decanoate 77.2 51.8 3.74 3.74 48
90.5 55.3 3.86 3.98 80
Most simulations were started with equal initial densities
C5 –C8 96.0 56.7 4.123 3.723 42,43
of the two boxes. The initial density was chosen such that if C5 –C8 116.0 56.8 3.70 3.70 42,43
at the given temperature the simulations would give coexist- C1000 57.0 57.0 4.28 4.28 81
ing liquid and gas densities equal to the experimental densi- C5 –C8 85.6 57.07 3.905 3.905 41
ties, the equilibrium volumes of the two boxes would be 88.1 59.4 3.905 3.905 19
C5 –C8 88.1 59.38 3.820 3.820 41
equal. Such an equilibrium configuration was subsequently
C50 60.1 60.1 3.80 3.80 82
used for some simulations at higher and lower temperatures. C5 –C8 92.0 65.5 4.323 3.523 42,43
We found that systems could get easily trapped in undesired, C4 72.0 72.0 3.923 3.923 18

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

far from equilibrium configurations when using initial densi- C4 84.0 84.0 3.923 3.923 17
ties which would yield widely differing gas and liquid vol- C4 –C48 114.0 47.0 3.93 3.93 this work
umes. Therefore, we used several systems with new initial
densities during the determination of the entire coexistence
curve. To generate the initial state, we placed the alkanes on
a lattice. For long chain alkanes it was important to ‘‘melt’’
r l 2 r g 5B ~ T2T c ! b , ~22!
this lattice using ordinary N,V,T simulations before the
Gibbs ensemble simulations were started. Immediately start- where b is the critical exponent. For small molecular fluids
ing with the Gibbs ensemble simulations made the system such as the Lennard-Jones fluid5,30,39 the data can be fitted
initially move far away from equilibrium and subsequently very well with an Ising-type critical exponent ~b50.32!. For
very long simulations were required to reach equilibrium. the short chain alkanes, C5 –C10 , we could fit the simulation
For short chains, however, we could start directly with the data and experimental data well with such an Ising exponent,
Gibbs-ensemble simulations. whereas a classical exponent ~b50.5! could not fit the data.
During the simulations the number of particles in the For the long chain alkanes, the data were not sufficiently
two boxes and the volumes of the boxes were stored. From accurate to distinguish between the two exponents. To be
these data we constructed a histogram of the densities and an consistent with the short chain lengths, we used the Ising
x – y plot.4,30 At sufficiently low temperatures, the two boxes exponent for all molecules. Note that this exponent is con-
of the Gibbs ensemble do not change identity. Once one of sistent with the experimental value for polymer solutions as
the boxes contains the liquid phase it will keep it during the determined by Dobashi et al.40
simulation. At these low temperatures, average densities in We tested our program by comparing our results with the
the two boxes are used as estimates of the coexistence den- Gibbs-ensemble simulations of an 8 bead Lennard-Jones
sities. The accuracy is estimated using the standard block polymer of Mooij et al.,6 with which they were in excellent
averaging techniques.1 Close to the critical point the boxes agreement. In addition, we compared our results with the
may switch identity. At those conditions the density histo- simulation data reported by Laso et al.8 The agreement was
grams are used and the coexistence densities are determined again very good. A more extensive discussion will be given
from the maxima of this density histogram. Estimates of the in the next section.
accuracy are made by dividing the simulations in blocks. The
x – y plots are used to judge the reliability of a simulation, for
each sample two points are plotted on the x – y plane IV. MODELS AND RESULTS
~x5n 1 /N, y5V 1 /V! and (12x,12y). From these plots one Over the last two decades various models to describe the
can observe whether a simulation was reliable.4 interaction between alkanes have been developed. Table I
The critical point was determined36 by fitting the coex- shows the energy and size parameters of the models that use
istence densities to the law of rectilinear diameters,37 a Lennard-Jones potential for the nonbonded interactions. A
comparison of the various models shows that for the size
r l1 r g parameter the variation can be as much as 20%, although
5 r c 1A ~ T2T c ! , ~21!
2 most models use a value of approximately 3.9 Å. For the
energy parameters there is little consensus on the preferred
where r l ( r g ) is the density of the liquid ~gas! phase, r c the values. Some models use the same value for e for the methyl
critical density, and T c the critical temperature. Furthermore, and methylene units while others use different ones. These
the results were also fitted to the scaling law for the density38 differences can amount to a factor of 2.

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995

 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria 2131

These large differences in parameters motivated this TABLE II. Results of the Gibbs-ensemble simulations for the n-alkanes as
study. The availability of techniques to determine the phase described with the OPLS model. T is the temperature, r g , r l are the densi-
ties of the gas and liquid phase, respectively, and ‘‘acc’’ is the probability of
behavior of the alkanes over a large temperature range and
a successful exchange between the two boxes. The subscripts give the ac-
for various chain lengths allows for an extensive comparison curacy of the last decimal~s!, i.e., 0.081530 is 0.081560.0030.
with experimental data.
In this study, we focus on united atom models. It is well T rg rl acc
known that united atom models fail to describe solid alkanes ~K! ~g cm23! ~g cm23! ~%!
correctly.26,27 Also for dense liquids Toxvaerd21 observed pentane C5
that it is impossible to describe the equation of state of the 350 0.00779 0.571 0.5
alkanes consistently with a united atom model. Recently, Pa- 402 0.0261 0.511 2.2
435 0.0595 0.471 4.1
dilla and Toxvaerd41 argued that it is even impossible to de- 450 0.072 0.434 8.1
scribe the second virial coefficient of the alkanes with such a 460 0.082 0.393 9.2
model using realistic parameters. This observation is surpris- 470 0.101 0.362 10.7
ing since the assumptions underlying the united atom model
octane C8
should hold very well for alkanes in a low density gas phase. 424 0.0033 0.6207 0.2
López Rodrı́guez et al.,42,43 however, have shown that it is 464 0.0083 0.5815 0.8
possible to describe the second virial coefficient accurately 488 0.0147 0.541 1.5
using a united atom model. It is therefore interesting to in- 512 0.0137 0.54510 1.8
536 0.02420 0.512 3.0
vestigate how well a united atom model can describe vapor–

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

560 0.0437 0.47510 5.1
liquid equilibria. 576 0.05510 0.422 5.8
584 0.082 0.432 6.5
592 0.092 0.402 8.9
A. The OPLS model 596 0.114 0.355 12.3

dodecane C12
550 0.0062 0.6152 0.4
1. The model
600 0.0102 0.5703 1.4
One of the most popular models used in simulations of 625 0.02010 0.551 1.4
650 0.0246 0.50510 3.1
alkanes and monolayers44 – 46 is based on the OPLS model of
665 0.04810 0.492 3.4
Jorgensen et al.19 This model has been further refined by 680 0.052 0.46510 4.0
Hautman and Klein44 to include bond bending. The OPLS 700 0.08520 0.44520 5.5
model uses a united atom description in which CH2 and CH3
hexadecane C16
groups are considered as one united atom. The nonbonded 625 0.0089 0.611 0.1
interactions between united atoms of different molecules and 675 0.0229 0.581 0.7
within a molecule ~if two atoms are more than four atoms 725 0.03510 0.521 1.8
apart! are described with a truncated Lennard-Jones potential 750 0.04910 0.471 3.4

FS D S D G
770 0.082 0.452 3.9
sij 12
sij 6
780 0.093 0.442 4.0
u LJ~ r i j ! 54 e i j 2 . ~23! 790 0.132 0.412 4.2
rij rij
The energy parameters of CH2 and CH3 groups are, respec-
tively, e CH2 5 59.4 ~K! and e CH3 5 88.1 ~K!. Throughout this
work, we use e i j 5 Ae i e j as the combining rule for the energy 5
parameters of the unlike interactions. The size parameters of u tors~ f ! 5 ( c k cosk ~ f ! , ~25!
the methylene and methyl groups are assumed to be equal k50
and have the value s53.905 Å. The potential is truncated at
where f is the dihedral angle. The parameters are c 0 51116
11.5 Å. No tail corrections have been applied. Note that tail
~K!, c 1 51462 ~K!, c 2 521578 ~K!, c 3 52368 ~K!,
corrections can have a significant effect on the phase dia-
c 4 53156 ~K!, and c 5 523788 ~K!.
gram, for example, for the Lennard-Jones fluid including tail
corrections this results in a critical temperature which is
;30% higher than without these corrections.39 The intramo- 2. Results and discussion
lecular interactions consist of bond bending and torsion.47 For the OPLS model, we used the bond bending and the
The distance between the atoms has been fixed to 1.53 Å. torsion potentials for the internal energy. In Appendix A, the
For the bond bending the van der Ploeg and Berendsen details on how the trial orientations are generated are de-
potential48 is used scribed.
1 The results of these simulations are shown in Table II. In
u bend~ u ! 5 k ~ u 2 u 0 !2, ~24! Fig. 3 the vapor–liquid curve as obtained from the simula-
2 u
tions using the OPLS model is compared with experimental
with k u 562 500 ~K rad22! and equilibrium angle u05112 data. This model has been fitted to the thermodynamic data
~deg!. For the torsion potential the original Ryckaert and of short chain alkanes at room temperature. Figure 3 shows
Bellemans17 potential is used that for pentane the agreement with experimental data is

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995

2132 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria

excellent. For larger chain lengths the deviation from experi-

mental data becomes significant. For hexadecane the coex-
istence curve is shifted approximately 100 ~K!. In Table III
the estimated critical points are compared with the experi-
mental values. For pentane the critical temperature is within
3% of the experimental value, for hexadecane the critical
point is overestimated by 10%. Since the difference between
the simulation results and the experimental data increases
with chain length, a simple rescaling of the energy param-
eters is not expected to give a better description.

B. The de Pablo model

1. The model
Laso et al.8 used the model introduced by de Pablo
et al.25 to calculate a coexistence point of several alkanes.
The model of de Pablo also uses a Lennard-Jones potential to
describe the nonbonded interactions between united atoms.

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

The energy parameters are e CH2 5 49.3 ~K! and e CH3
5 90.5 ~K! and the size parameters have the same value for
the methylene and methyl groups, namely s53.94 Å. For the
truncation of the potential the minimum image convention
was used. The bond length was fixed at 1.53 Å and the bond
angle was constrained to 112 deg. The Jorgensen torsion po-
tential was used,
FIG. 3. Vapor–liquid equilibria of various alkanes as obtained from the
Gibbs-ensemble simulations ~open triangles! using the OPLS model. The u tors~ f ! 5c 0 10.5c 1 ~ 11cos f ! 10.5c 2 ~ 12cos 2 f !
small dots are experimental data ~for C5 –C9 the data are taken from Ref. 84
for larger chain lengths data are estimated from an equation of state!. The 10.5c 3 ~ 11cos 3 f ! , ~26!
large dots are the experimental critical points ~Ref. 76!. The filled triangle is
the estimate of the critical point based on the simulation data. with c 1 5355 ~K!, c 2 5268.19 ~K!, and c 3 5791.3 ~K!.

TABLE III. The critical points of the various models as calculated from the Gibbs-ensemble simulations. The
experimental critical points are from Ref. 76 and the experimental critical pressures from Refs. 83,54. T c is the
critical temperature, r c the critical density, and P c the critical pressure. The subscripts give the accuracy of the
last decimal~s!.

T c ~sim! T c ~exp! r c ~sim! r c ~exp! P c ~sim! P c ~exp!

~K! ~K! ~g cm23! ~g cm23! ~MPa! ~MPa!

OPLS model
C5 4814 469.7 0.23410 0.230
C8 6165 568.6 0.2348 0.232
C12 7429 658.2 0.2348 0.226
C16 8086 723.0 0.25410 0.219

Toxvaerd model
C5 4393 469.7 0.2258 0.230
C8 5327 568.6 0.23213 0.232
C12 5923 658.2 0.20510 0.226

de Pablo model
C8 58411 568.6 0.22110 0.232
C24 82314 n.a. 0.21914 n.a.

New model
C5 4944 469.7 0.2235 0.230 3.95 3.369
C6 5234 507.0 0.2265 0.233 3.25 3.014
C7 5564 539.8 0.2325 0.233 3.15 2.734
C8 5773 568.6 0.2295 0.232 2.75 2.485
C10 6046 617.5 0.2294 0.228 2.35 2.099
C12 6597 658.2 0.2237 0.226 2.35 1.810
C16 7195 723.0 0.2186 0.219 1.95 1.401
C24 7968 n.a. 0.2059 n.a. 1.35 n.a.
C48 92411 n.a. 0.19514 n.a. 1.06 n.a.

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995

 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria 2133

TABLE IV. Results of the Gibbs-ensemble simulations for the de Pablo

model ~see also the caption to Table II!.

T rg rl acc
~K! ~g cm23! ~g cm23! ~%!

octane C8
473 0.0175 0.531 0.8
498 0.0215 0.5051 2.6
523 0.041 0.461 4.7
548 0.062 0.422 7.1
563 0.082 0.384 14.3

tetracosane C24
650 0.0042 0.53610 0.3
675 0.0094 0.521 1.0 FIG. 4. Schematic sketch of Toxvaerd’s model ~a!. The interaction sites of
700 0.0184 0.50710 1.1 the nonbonded interactions ~filled circles! are displaced ~to the position of
750 0.03510 0.44510 2.4 the valence electrons! with respected to the center of the carbons ~open
775 0.05410 0.402 3.1 circles!. ~b! gives the approximated united atom model in which the inter-
action sites are at the same position as the center of the carbons.

2. Results and discussion the scheme advocated by de Pablo et al., however, the tor-

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

sional potential appears in the Rosenbluth factor and hence
The vapor–liquid equilibrium densities as obtained from
in the acceptance rule. As a consequence, the probability of
the model of de Pablo et al. are given in Table IV. Our data
acceptance will be less than unity. Moreover, since the
are in excellent agreement with those reported by Laso et al.8
Rosenbluth factor is a product of the Boltzmann factors of
The estimated critical points are listed in Table III. The
the torsional potentials, this probability will decrease rapidly
model proposed by de Pablo et al. gives a better overall de-
with chain length. For a system with external interactions,
scription than the OPLS model ~see Sec. IV A 2!.
our scheme has the additional advantage that we only calcu-
Laso et al. used a similar method to calculate the vapor–
late the ~expensive! external interactions for trial orientations
liquid curve. This method has been developed by de Pablo
that already have an ‘‘optimal’’ torsion potential. In the
et al.22 and is referred to as continuum-configurational-bias
scheme of de Pablo et al.,22 most trial orientations have such
Monte Carlo. This method also combines the Gibbs-
a high torsional potential that they have a very low probabil-
ensemble technique with the Rosenbluth algorithm to insert
ity of being accepted, yet for all these orientations the non-
chain molecules. An important issue pointed out by Laso
bonded interactions have to be calculated. These two factors
et al.8 is that the continuum-configurational-bias Monte
make our scheme already an order of magnitude more effi-
Carlo becomes computationally expensive for systems of
cient for C15 and up to several orders of magnitude for the
pure alkanes of more than about twenty segments. Compari-
longer chains. In addition, the approach of de Pablo et al. is
son of the acceptance probability of octane with the corre-
very inefficient for potentials that are strongly peaked such
sponding one of tetracosane ~C24! shows that in our version
as bond bending and bond vibration.
this probability does not decrease significantly. To see the
reason for this, it is instructive to compare the two schemes
in some detail. C. The Toxvaerd model
The difference between our algorithm and the scheme 1. The model
proposed by de Pablo et al.22 is the method in which the trial
orientations are generated. In the model used by Laso et al.8 Toxvaerd21,49 introduced an anisotropic potential to
the bond angle and bond length are fixed and therefore one model the effects of hydrogen on the thermodynamic prop-
has to generate only the torsional angle f. De Pablo et al. erties without increasing the number of interaction sites. In
generate the first torsional angle ~f1! at random and the other this model, the interaction site of the nonbonded Lennard-
n21 angles are calculated from Jones potential is displaced with respect to the center of mass
of the carbon atoms ~see Fig. 4!,

FS D S D G
f i11 5 f i 12 p /n.
sij 12
sij 6
Laso et al. used n512 which gives twelve equally spaced u LJ~ R i j ! 54 e i j 2 , ~27!
Rij Rij
trial orientations. Note that in the Rosenbluth factor of the de
Pablo scheme the torsional potential has to be included. where R i j is the distance between the interaction sites. The
To compare the efficiency of the scheme of de Pablo relation between Ri and the centre of mass ri of atom i is
et al. and the scheme utilized in this work, we consider a given by
model which has only internal interactions ~i.e., only the tor- ri 20.5 ~ ri11 2ri !
sional potential!. For our algorithm this implies that there are Ri 5ri 1d . ~28!
u ri 20.5 ~ ri11 2ri ! u
no external interactions and hence the Rosenbluth factors of
all generated conformations are by definition one. Therefore, Padilla and Toxvaerd49 use d50.37 ~Å! for the CH2 groups
all conformations that are generated will be accepted with and d50.275 ~Å! for the CH3 groups. The parameters of the
probability one irrespective of the length of the molecule. In Lennard-Jones potential are e CH3 5 120 ~K!, e CH2 5 80 ~K!,

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995

2134 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria

TABLE V. Results of the Gibbs-ensemble simulations for the Toxvaerd M k

model ~see also the caption to Table II!.
W̄ ~ n ! 5exp@ 2 b ū 1 ~ n !# )( exp@ 2 b ū j ~ n i !# . ~30!
T rg rl acc j52 i51
~K! ~g cm23! ~g cm23! ~%!
In addition, we also calculate the difference in energy be-
pentane C5 tween the Toxvaerd potential and the approximate potential
300 0.0051 0.6085 0.1
of the molecule in the selected conformation
325 0.0082 0.5805 0.3
350 0.0123 0.541 0.9
d U ~ n ! 5U ~ n ! 2Ū ~ n ! . ~31!
375 0.0285 0.512 1.5
400 0.04710 0.452 3.5
For the old conformation, we determine the Rosenbluth fac-
410 0.06310 0.431 4.5
420 0.072 0.392 7.6 tor W̄(o) using the approximate potential and we calculate
430 0.092 0.342 10.5 the energy difference between the two potentials of the old
conformation
octane C8
425 0.0155 0.55512 0.4 d U ~ o ! 5U ~ o ! 2Ū ~ o ! . ~32!
450 0.02710 0.52510 0.8
475 0.04710 0.491 1.4 If the move is the regrowing of part of a molecule, it is
500 0.05010 0.42510 5.1
520 0.092 0.351 9.6
accepted with a probability

S W̄ ~ n !
D

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

dodecane C12
450 0.0043 0.611 0.03 acc~ o→n ! 5min 1, exp$ 2 b @ d U ~ n ! 2 d U ~ o !# % .
475 0.0115 0.5806 0.06 W̄ ~ o !
500 0.0105 0.5327 0.3 ~33!
525 0.0156 0.492 0.6
550 0.0298 0.452 1.5 If the move involves an exchange of a molecule between the
two boxes, the acceptance rule is

S
575 0.04512 0.382 6.0
590 0.167 0.283 7.0
V 1 ~ N2n 1 ! W̄ ~ n !
acc~ o→n ! 5min 1,
~ V2V 1 !~ n 1 11 ! W̄ ~ o !

and s CH35sCH2 5 3.527 ~Å!. The potential is truncated at

12.0 ~Å! and the usual tail corrections are applied.1
3exp$ 2 b @ d U ~ n ! 2 d U ~ o !# % . D ~34!
The intramolecular interactions include bond bending
and torsion. For the bond bending Eq. ~24! is used with
k u 562 500 ~K! and u05113.3 deg. For the torsion, Eq. ~25!
was used with the parameters proposed by Padilla and
Toxvaerd,49 namely c 0 51038 ~K!, c 1 52426 ~K!,
c 2 581.6 ~K!, c 3 523129 ~K!, c 4 52163 ~K!, and
c 5 52252 ~K!. The bond length was fixed to 1.539 ~Å!.50

2. Simulation details
The configurational-bias Monte Carlo scheme can not be
applied directly to the Toxvaerd model because we have to
know the position of atom l11 to determine the position of
the interaction site of atom l. We have used the following
algorithm to make configurational-bias Monte Carlo simula-
tions for this model possible.
Let us define an approximate potential, denoted by ū,
which is identical to the Toxvaerd model but the interaction
site of the nonbonded interaction is at the position of the
carbon atoms. Hence the approximate model is an ordinary
united-atom model ~see Fig. 4! for which we can use the
configurational-bias technique as described in Sec. II B. The
probability of generating conformation n is given by
exp@ 2 b Ū ~ n !#
P ~ o→n ! 5 , ~29!
W̄ ~ n !
where Ū(n)5 ( i51
M
ū i . The bar above the symbols indicate
that this property is calculated with the approximate poten- FIG. 5. Vapor–liquid equilibria of Toxvaerd’s model ~see also the caption to
tial. The Rosenbluth factor is given by Fig. 3!.

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995

 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria 2135

TABLE VI. Results of the Gibbs-ensemble simulations using the model TABLE VI. ~Continued.!
presented in this work ~see also the caption to Table II!.
T rg rl acc
T rg rl acc ~K! ~g cm23! ~g cm23! ~%!
~K! ~g cm23! ~g cm23! ~%!
dodecane C12
pentane C5 450 0.0055 0.6255 0.03
350 0.0075 0.5565 0.4 475 0.0054 0.6025 0.2
365 0.0095 0.5435 0.6 500 0.0125 0.5765 0.5
375 0.0145 0.5285 0.8 525 0.0155 0.5455 1.2
385 0.0225 0.525 1.1 550 0.0145 0.5185 2.4
400 0.02115 0.501 1.6 575 0.02410 0.4941 2.0
585 0.04410 0.46610 3.5
415 0.0265 0.4825 2.1
600 0.06010 0.43720 4.8
425 0.03010 0.47510 2.5
615 0.05510 0.451 3.6
435 0.04410 0.45810 3.5
440 0.04610 0.44610 5.4 hexadecane C16
445 0.04710 0.43010 5.1 550 0.005810 0.5815 0.2
455 0.05510 0.411 6.8 575 0.005610 0.551 0.9
465 0.072 0.402 6.7 600 0.0145 0.52610 1.0
475 0.091 0.3801 7.4 615 0.0225 0.51410 1.5
625 0.021 0.51810 1.2
hexane C6 640 0.0305 0.48310 2.4

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

350 0.0032 0.5965 0.1 650 0.0335 0.47710 2.2
375 0.0083 0.5777 0.2 660 0.02210 0.442 4.8
385 0.0093 0.5657 0.3 675 0.0475 0.4237 6.1
400 0.0105 0.54410 0.6 685 0.061 0.39610 6.5
415 0.0175 0.5275 1.0
425 0.0175 0.5155 1.3 tetracosane C24
625 0.0043 0.5595 0.2
450 0.0285 0.481 2.7
650 0.0105 0.54710 0.3
465 0.0425 0.4645 3.2
665 0.0135 0.53210 0.5
475 0.0505 0.4435 5.1
675 0.0105 0.51710 0.8
500 0.073 0.393 7.3
685 0.01110 0.50310 1.2
heptane C7 700 0.01610 0.481 1.1
375 0.004250 0.6055 0.1 715 0.01610 0.45415 1.6
385 0.004850 0.5955 0.2 725 0.02210 0.43515 2.3
735 0.03410 0.451 1.0
415 0.008730 0.5603 0.4
750 0.04010 0.401 3.7
425 0.0113 0.5515 0.6
435 0.0143 0.5434 0.8 octatetracontane C48
450 0.0213 0.53110 1.1 800 0.0053 0.46510 0.2
475 0.0333 0.4974 2.1 825 0.0115 0.45310 0.3
490 0.0417 0.4797 2.9 850 0.02010 0.41910 0.4
500 0.04610 0.4627 4.0 875 0.03610 0.39315 1.1
515 0.0605 0.4325 5.3 890 0.05520 0.353 1.4
900 0.063 0.323 1.9
octane C8
400 0.0041 0.6053 0.1
425 0.0075 0.5857 0.2
435 0.0085 0.5737 0.8 In Appendix B it is proven that this scheme indeed samples
450 0.0125 0.5527 1.0
the desired distribution of configurations.
460 0.0133 0.5454 1.0
475 0.0155 0.5267 2.2 This method is similar to what can be used for systems
485 0.0283 0.5214 1.4 with ‘‘expensive’’ potentials. In such a model one can grow
490 0.0275 0.5177 1.8 the molecules with an approximate potential which is very
503 0.0355 0.4917 2.9 ‘‘cheap.’’ The correct energy of the conformation, as given
510 0.0333 0.4685 3.5 by the expensive potential, is only calculated once, namely
523 0.0475 0.4577 4.7 for the selected conformation and not for every trial orienta-
535 0.0525 0.4398 5.8
tion.
550 0.081 0.401 6.8

decane C10
3. Results
425 0.0022 0.62010 0.1
450 0.0052 0.59510 0.2 We have calculated the vapor–liquid curves of pentane,
475 0.0102 0.57410 0.6 octane, and dodecane. The results are presented in Table V.
500 0.0145 0.53110 1.5 The simulation results are compared with experimental data
515 0.0185 0.51110 2.3 in Fig. 5. The critical points are given in Table III. For pen-
530 0.0305 0.49910 2.5
tane Toxvaerd’s model predicts the critical point at a much
550 0.0487 0.45415 4.7
575 0.071 0.42510 5.8 lower temperature than the experimental one. Note that al-
though for small chain lengths results obtained with the

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995

2136 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria

1. The model
The nonbonded interactions between the united atoms
are described with a Lennard-Jones potential where the en-
ergy parameters are e CH2 5 47.0 ~K! and e CH3 5 114 ~K!.
The size parameters have the same value for the methylene
and methyl groups, namely s53.93 Å. The potential was
truncated at 13.8 Å and the usual tail corrections were ap-
plied.
The bond-bending potential48 is of the form of Eq. ~24!
with k u 562 500 ~K rad22! and equilibrium angle u05114
deg. The torsion potential19 is of the form of Eq. ~26! with
parameters c 1 5355 ~K!, c 2 5268.19 ~K!, and c 3 5791.3
~K!. We also tested the Ryckaert and Bellemans torsion po-
tential ~25!, but no significant differences in the phase behav-
ior could be observed.

2. Results and discussion

In Table VI, the results of the simulations are summa-

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

rized. In Fig. 6, the results of the simulations are compared
with data for the n-alkanes for which either experimental
data are available or can at least be estimated with some
reliability ~C5 –C16!, the overall agreement with experimental
data is surprisingly good. In Table III the estimated critical
properties are listed. This table shows that for pentane the
critical temperature is slightly overestimated.
FIG. 6. Vapor–liquid equilibria of the model presented in this work ~see These results show that it is possible to model the phase
also the caption to Fig. 3!.
behavior of the n-alkanes over a large temperature range
with a united atom model. The density appears to be suffi-
ciently low so that it is not necessary to model the hydrogens
OPLS model agree better with experimental data than those explicitly. To obtain this agreement, a large difference be-
obtained with Toxvaerd’s model, the difference between ex- tween the energy parameters of the CH2 interactions and the
perimental data and the predictions of Toxvaerd’s model CH3 interaction was required. Similar conclusions have been
does not increase with chain length. This suggests that by obtained by López Rodrı́guez et al.,42,43 Almarza et al.,51 and
rescaling of the parameters a better quantitative agreement de Pablo and co-workers.8,25 Padilla and Toxvaerd,41 how-
may be obtained. ever, argued that a ratio of e CH3 /eCH2 that is much larger than
1.5 is unphysical. López Rodrı́guez et al.42,43 showed that a
high ratio of e CH3 /eCH2 is to some extent due to the assump-
tion s CH35sCH2 . If s CH3 is taken as larger than s CH2 , the
D. New model
ratio e CH3 /eCH2 would not be as large. It would be interesting
Comparison of the parameters of the nonbonded interac- to investigate this in further detail.
tions of the OPLS model with those used by de Pablo et al.25
shows that in their study the ratio e CH3 /eCH2 is much larger V. CRITICAL PROPERTIES OF THE ALKANES
than the corresponding ratio for the OPLS model. Further-
Alkanes are thermally unstable above approximately 650
more, the e CH2 has a much larger value in the OPLS than in
~K!, which makes experimental determination of the critical
the model of de Pablo et al. Since this parameter determines points of alkanes longer than decane ~C10! extremely diffi-
to a large extent the value of the critical temperature, it be- cult. Long alkanes, however, are present in mixtures of prac-
comes clear why the OPLS model predicts a much higher tical importance for the petrochemical industry. In these mix-
critical point than the model of de Pablo et al. A similar set tures, the number of components can be so large that it is not
of parameters has been obtained by López Rodrı́guez practical to determine all phase diagrams experimentally.
et al.42,43 from a study of the virial coefficients of alkanes. To One therefore has to rely on predictions made by equations
describe these virial coefficients accurately López Rodrı́guez of state. The parameters of these equations of state are di-
et al. had to introduce a low value of e CH2 and a large dif- rectly related to the critical properties of the pure compo-
ference between e CH2 and e CH3 . We used the observations of nents. Therefore, the critical properties of the long-chain al-
López Rodrı́guez et al. and de Pablo et al. as a starting point kanes are essential in the design of petrochemical processes,
in investigating whether a united-atom model can give a even if they are unstable close to the critical point. However,
good description of a large range of alkanes. The model de- as experimental data are scarce and contradictory, we had to
scribed below gave a good overall description of the phase rely on semiempirical methods to estimate the critical
behavior. properties.52–55

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995

 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria 2137

FIG. 7. Critical temperature T c as a function of carbon number N c . n, the FIG. 9. Critical pressure P c as a function of carbon number N c . n, the
simulation data and d, the experimental data from Refs. 76,55. simulation data and d, the experimental data from Refs. 83,54.

The experimental data available to Kurata and Isida

More than 40 years ago, the critical properties and the
showed that the critical density was independent of chain

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

equation of state of n-alkanes were already the topic of theo-
length, suggesting that both the theory of Prigogine and of
retical investigations. Most of these studies were aimed at
Flory–Huggins were not directly applicable to these systems.
establishing empirical relationships between the critical
To take into account this experimental fact, Kurata and Isida
properties and carbon number56 – 62 based on the compilation
made an ad hoc assumption on the scaling behavior of the
of the available experimental data by Egloff.63 A more fun-
critical properties such that r c }n 1/3 /(n 1/3 11). Kreglewski
damental approach was taken by Prigogine and
and Zwolinski,73 Nakanishi et al.,74 and more recently
co-workers,64 – 66 who extended the cell method to n-mers to
Tsonopoulos52 used Kurata and Isida’s empirical corrections
derive an equation of state for these components. Hijmans67
of the Flory scaling relations to correlate the various critical
used the results of Prigogine and co-workers to derive phe-
properties successfully.
nomenological relations for the chain length dependence of
Prigogine’s treatment of r-mers was later revisited by
the thermodynamic properties. Prigogine’s treatment of
Flory, Orwoll, and Vrij,75 this study showed that the lattice
r-mers predicts that the critical temperature scales as
treatment inherent in the cell theory ~which by fixing the
T c }n/(n 1/2 11) 2 and the critical density as r c }n 21/2 ,
nearest neighbors of a given molecule exactly at their mean
where n is the number of monomeric units in the chain.
separation, suppresses the randomness which is the foremost
Kurata and Isida68 assumed the vapor–liquid equilibria
characteristic of the liquid state!. Flory and co-workers75
to be identical to a solution of rodlike polymers in a solvent
then continued and derived a continuum theory. This theory
of small molecules. The chain length dependence of the criti-
was recently used by Tsonopoulos and Tan55 to describe the
cal properties of the n-alkanes is, with this assumption, iden-
more recent experimental data successfully.
tical to the dependence of a polymer solution. Interestingly,
Experimentally, the critical properties of n-alkanes up to
the Flory–Huggins theory69–72 for polymer solutions also
C18 have been studied by Anselme et al.76 ~see Figs. 7 and
predicts that the critical density decreases with chain length,
8!. The most often used extrapolations assume that the criti-
i.e., r c }n 21/2 .
cal density is a monotonically increasing function of the car-
bon number, approaching a limiting value for the very long
alkanes.52,55 In contrast to the expectations which are based
on these extrapolations, the experimental data of Anselme
et al.76 indicate that the critical density has a maximum for
C8 and then decreases monotonically. The experimental data
of Steele ~as reported in by Tsonopoulos and Tan55!, how-
ever, do not provide any evidence for such a maximum ~see
Fig. 8!.
Since we can use our simulation technique to study
phase behavior of the longer alkanes at conditions where
experiments are not ~yet! feasible, we are in a position to
make predictions of the critical properties of these mol-
ecules. Figure 7 shows that our calculations of the critical
temperatures are in very good agreement with both the data
of Steele and Anselme et al. The simulation results for the
critical densities ~Fig. 8! show the same trend as observed by
FIG. 8. Critical density r c as a function of carbon number N c . n, the
simulation data and d, the experimental data from Anselme et al. ~Ref. 76! Anselme et al. and therefore strongly support these experi-
and h, the data of Steele ~as published in Ref. 55!. ments. At this point it is interesting to note that Mooij et al.6

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995

2138 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria

and Sheng et al.77 used Monte Carlo simulations to study the ACKNOWLEDGMENTS
vapor–liquid curve of a polymeric bead-spring model for
We would like to thank M. J. Anselme, D. Frenkel,
various chain lengths. These studies also show a decrease of
W. G. Heitman, A. K. van Helden, H. P. C. E. Kuipers, J.-P.
the critical density as a function of chain length. This indi-
Ryckaert, S. Schreuder, and M. Sprik for their contributions
cates that the decrease of the critical density with chain
to this work.
length is a more general feature of chain molecules that does
not depend on the details of a particular model.
APPENDIX A: GENERATION OF TRIAL ORIENTATIONS
The results for the critical pressure are presented in
Table III. The critical pressure was calculated from fitting the In this Appendix, we demonstrate the way we generate
vapor pressure data of the simulations to the Clausius– the trial orientations in the configurational-bias Monte Carlo
Clapeyron equation. This equation was then used to extrapo- scheme.
late to the critical point. Comparison with the experimental Let us first consider the general case with flexible bond
data ~see Fig. 9! shows that, considering the accuracy of the length, bond bending, and torsion. The probability that we
data, the agreement between the simulations and experiments generate a trial configuration b is given by
is very good. exp@ 2 b u int~ b!# db
P ~ b! db5 , ~A1!
VI. CONCLUDING REMARKS C

In this work we have used the Gibbs-ensemble technique where C is a normalization constant which is defined by
in combination with the configurational-bias Monte Carlo
E

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

method to determine the vapor–liquid curves of various C[ db exp@ 2 b u int~ b!# . ~A2!
b
n-alkanes. Different alkane models have been compared and
a new model is introduced that can describe the vapor–liquid Note that in the configurational-bias Monte Carlo scheme we
curve over a large temperature range for a large number of do not have to calculate this constant.
alkanes. It is convenient to represent the position of a atom using
Whereas the conventional Gibbs-ensemble technique is the bond length r, bond angle u, and torsional angle f ~see
limited to butane or pentane, the combination with Fig. 10!. With these coordinates the volume element db is
configurational-bias Monte Carlo allows for the simulation given by
of chains as long as C48 . On an IBM/340 workstation such a db5r 2 cos~ u ! drd u d f . ~A3!
simulation takes approximately 1 week of cpu time, for oc-
tane ~C8! the corresponding cpu time is approximately 12 h. The internal energy is the sum of the bond vibration poten-
Note that the increase of cpu time for the long chain alkanes tial, the bond-bending potential, and the torsion potential,
is mostly due to the increase of the number of atoms ~for C48 u int~ r, u , f ! 5u vib~ r ! 1u bend~ u ! 1u tors~ f ! . ~A4!
we use ;5000 atoms and for C8 ;1000!. Since the probabil-
ity of a successful exchange between the liquid and vapor Substitution of Eqs. ~A4! and ~A3! into Eq. ~A1! gives
phase does not lessen significantly in our scheme, we expect P ~ b! db5 P ~ r, u , f ! r 2 drd u d f
that it is possible to determine the coexistence curve of even
longer chains. 5exp@ 2 b u bond-vib~ r !# r 2 dr
This work demonstrates that for modeling vapor–liquid 3exp@ 2 b u bend~ u !# cos~ u ! d u
coexistence a relatively simple united-atom model is suffi-
cient to obtain a very good agreement with experimental data 3exp@ 2 b tors~ f !# d f . ~A5!
and it is not necessary to take the hydrogen atoms explicitly In our simulations we have used an alkane model with fixed
into account. To get this agreement it was necessary to make bond length, therefore in our case the first term in Eq. ~A5! is
the energy parameters of the nonbonded potential of the a constant.
CH3 –CH3 interaction very different from the corresponding For the second carbon atom there are no internal inter-
value for the CH2 –CH2 interaction. This observation is in actions other than the constraints of the bond length. The
agreement with the conclusions of other simulation distribution of trial orientations, Eq. ~7! reduces to
studies.8,25,42,51
For petrochemical applications knowledge of the critical
properties of the n-alkanes is of interest even at temperatures
where these molecules are thermally unstable. Even qualita-
tive aspects, such as the chain-length dependence of the criti-
cal properties, are poorly understood for these systems. Our
calculations show that, in contrast to the traditional view, the
critical density of the long alkanes decreases rather than in-
creases with carbon number. The simulations presented in
this work show that it is possible to use simulations as an
‘‘engineering tool’’ to generate reliable data for the critical
properties of the n-alkanes at conditions where experiments FIG. 10. Definition of the bond length r, bond angle u, and torsional angle
are not ~yet! feasible. f.

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995

 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria 2139

P 2 ~ b! db}cos~ u ! d u d f . ~A6! APPENDIX B: PROOF OF ALGORITHM FOR

TOXVAERD’S POTENTIAL
Hence, the trial orientations are randomly distributed on the In this Appendix we prove that the algorithm of Sec.
surface of a sphere. The algorithm that we have used for IV C 2 for the Toxvaerd potential gives the desired distribu-
generating random vectors on the surface of a sphere is de- tion of configurations. The flow of configurations from state
scribed in Ref. 1. o to state n is given by
For the third atom, the internal energy contains the bond-
K ~ o→n ! 5N ~ o ! 3 p ~ o→n ! 3acc~ o→n ! . ~B1!
bending energy as well. This gives for the distribution of trial
orientations Imposing detailed balance and substitution of Eqs. ~29! and
~30! gives as condition for the acceptance rule
P 3 ~ b! db}exp@ 2 b u bend~ u !# cos~ u ! d u d f . ~A7! acc~ o→n ! exp@ 2 b U ~ n !#
5
acc~ n→o ! exp@ 2 b U ~ o !#
To generate k trial orientations that are distributed according
to Eq. ~7! we generate again a random vector on a unit exp@ 2 b Ū ~ n !# W̄ ~ o !
sphere and determine the angle u. This vector is accepted 3
W̄ ~ n ! exp@ 2 b Ū ~ o !#
with a probability exp@2 b u bend~u!#. If rejected, this proce-
dure is repeated until a value of u has been accepted. In Ref. W̄ ~ n !
5 exp$ 2 b @ d U ~ n ! 2 d U ~ o !# % . ~B2!
78 it is shown that this acceptance/rejection method indeed W̄ ~ o !
gives a distribution of trial orientations given by Eq. ~7!.

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

Acceptance rule ~33! obeys this condition which proves that
Note that the term cos u is taken into account by generating
the correct distribution is sampled.
a random vector on a sphere. In this way, k ~or k21 for the
case of the old conformation! trial orientations are generated. 1
M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids ~Claren-
An alternative scheme would be to generate angle u uni- don, Oxford, 1987!.
formly ~P@0,p#! and the bond-bending energy corresponding 2
A. Z. Panagiotopoulos, Mol. Phys. 61, 813 ~1987!.
to this angle is calculated. This angle u is accepted with a 3
A. Z. Panagiotopoulos, N. Quirke, M. Stapleton, and D. J. Tildesley, Mol.
Phys. 63, 527 ~1988!.
probability cos~u!exp@2 b u bend~u!#. If rejected, this proce- 4
B. Smit, Ph. de Smedt, and D. Frenkel, Mol. Phys. 68, 931 ~1989!.
dure is repeated until a value of u has been accepted. The 5
A. Z. Panagiotopoulos, Mol. Simul. 9, 1 ~1992!.
selected value of u is supplemented with a randomly selected 6
G. C. A. M. Mooij, D. Frenkel, and B. Smit, J. Phys. Condensed Matter 4,
angle f. These two angles determine a new trial orientation. 7
L255 ~1992!.
G. C. A. M. Mooij, Ph.D. thesis, Rijksuniversiteit Utrecht, The Nether-
For the fourth and higher carbon atoms, the internal en- lands, 1993.
ergy includes both bond-bending and torsion energy. This 8
M. Laso, J. J. de Pablo, and U. W. Suter, J. Chem. Phys. 97, 2817 ~1992!.
gives for Eq. ~7!, 9
J. Harris and S. A. Rice, J. Chem. Phys. 88, 1298 ~1988!.
10
J. I. Siepmann and D. Frenkel, Mol. Phys. 75, 59 ~1992!.
11
J. J. de Pablo, M. Bonnin, and J. M. Prausnitz, Fluid Phase Equil. 73, 187
l ~ b ! db}exp@ 2 b u bend~ u !#
p int ~1992!.
12
D. Frenkel, G. C. A. M. Mooij, and B. Smit, J. Phys. Condensed Matter 4,
3exp@ 2 b u tors~ f !# cos~ u ! d u d f . ~A8! 3053 ~1992!.
13
R. F. Cracknell, D. Nicholson, N. G. Parsonage, and H. Evans, Mol. Phys.
71, 931 ~1990!.
We again generate a random vector on a sphere and calculate 14
J. I. Siepmann, S. Karaborni, and B. Smit, J. Am. Chem. Soc. 115, 6454
the bond-bending angle u and torsion f. These angles are ~1993!.
accepted with a probability exp$2b@u bend~u!1u tors~f!#%. If 15
J. I. Siepmann, S. Karaborni, and B. Smit, Nature 365, 330 ~1993!.
these angles are rejected, a new vectors are generated until
16
D. E. Williams, J. Chem. Phys. 47, 4680 ~1967!.
17
J. P. Ryckaert and A. Bellemans, Chem. Phys. Lett. 30, 123 ~1975!.
one gets accepted. 18
J. P. Ryckaert and A. Bellemans, Faraday Discuss. Chem. Soc. 66, 95
Again the alternative scheme would be first to determine ~1978!.
a bond-bending angle u by generating u uniformly on @0,p# 19
W. L. Jorgensen, J. D. Madura, and C. J. Swenson, J. Am. Chem. Soc.
and calculating the bond-bending energy corresponding to 106, 6638 ~1984!.
20
W. L. Jorgensen and J. Tirado-Rives, J. Am. Chem. Soc. 110, 1657 ~1988!.
this angle. This angle u is then accepted with a probability 21
S. Toxvaerd, J. Chem. Phys. 93, 4290 ~1990!.
cos~u!exp@2 b u bend~u!#. This procedure is continued until we 22
J. J. de Pablo, M. Laso, and U. W. Suter, J. Chem. Phys. 96, 2395 ~1992!.
have accepted an angle. Next generate a torsion angle ran- 23
P. M. Rodger, Mol. Phys. 76, 1385 ~1992!.
24
domly on @0,2p# and accept this angle with a probability A. Berker, S. Chynoweth, U. C. Klomp, and Y. Michopoulos, J. Chem.
Soc. Faraday Trans. 88, 1719 ~1992!.
exp@2 b u tors~f!#, again repeating this until a value has been 25
J. J. de Pablo, M. Laso, J. I. Siepmann, and U. W. Suter, Mol. Phys. 80, 55
accepted. In this scheme the bond angle and torsion are gen- ~1993!.
erated independently which can be an advantage in cases 26
J. P. Ryckaert and M. L. Klein, J. Chem. Phys. 85, 1613 ~1986!.
27
J. P. Ryckaert, M. L. Klein, and I. R. McDonald, Phys. Rev. Lett. 58, 698
where the corresponding potentials are sharply peaked.
~1987!.
In the de Pablo model the bond angle is fixed. For the 28
M. A. Moller, D. J. Tildesley, K. S. Kim, and N. Quirke, J. Chem. Phys.
generation of the trial orientations this implies that Eqs. ~A7! 94, 8390 ~1991!.
and ~A8! must be replaced by an algorithm that generates
29
S. Karaborni and S. Toxvaerd, J. Chem. Phys. 96, 5505 ~1992!.
30
B. Smit, in Computer Simulation in Chemical Physics, NATO ASI, edited
orientations on the surface of a cone. The approximated po- by M. P. Allen and D. J. Tildesley ~Kluwer, Dordrecht, 1993!, pp. 173–
tentials for Toxvaerd’s model and the model introduced in 209.
31
this work are of the same form as the Jorgensen potential. J. I. Siepmann, in Computer Simulation of Biomolecular Systems: Theo-

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995

2140 Smit, Karaborni, and Siepmann: Vapor–liquid phase equilibria

54
retical and Experimental Applications, edited by W. F. van Gunsteren, P. A. S. Teja, R. J. Lee, D. Rosenthal, and M. Anselme, Fluid Phase Equi-
K. Weiner, and A. J. Wilkinson ~Escom Science, Leiden, 1993!, pp. 249– libria 56, 153 ~1990!.
264. 55
C. Tsonopoulos and Z. Tan, Fluid Phase Equilibria 83, 127 ~1993!.
32
D. Frenkel, in Computer Modelling of Fluids, Polymers and Solids, NATO 56
L. Grunberg and A. H. Nissan, Nature 161, 170 ~1948!.
ASI, edited by C. R. A. Catlow ~Kluwer, Dordrecht, 1990!. 57
L. Grunberg and A. H. Nissan, Trans. Faraday Soc. 44, 1013 ~1948!.
33
B. Smit, Ph.D. thesis, Rijksuniversiteit Utrecht, The Netherlands, 1990. 58
S. S. Mitra, J. Chem. Phys. 21, 2234 ~1953!.
34
D. Frenkel, in Computer Simulation in Chemical Physics, NATO ASI, 59
Y. P. Varshi, J. Chem. Phys. 21, 1400 ~1953!.
edited by M. P. Allen and D. J. Tildesley ~Kluwer, Dordrecht, 1993!, pp. 60
L. Grunberg, J. Chem. Phys. 22, 157 ~1954!.
93–152. 61
Y. P. Varshi, J. Chem. Phys. 22, 150 ~1953!.
35
B. Smit and D. Frenkel, J. Chem. Phys. 94, 5663 ~1991!. 62
Y. P. Varshi and S. N. Srivastava, J. Chem. Phys. 22, 349 ~1954!.
36
B. Smit and C. P. Williams, J. Phys. Condensed Matter 2, 4281 ~1990!. 63
G. Egloff, Physical Constants of Hydrocarbons ~Reinhold, New York,
37
J. S. Rowlinson and F. L. Swinton, Liquids and Liquid Mixtures, 3rd ed. 1939!.
~Butterworth, London, 1982!. 64
I. Prigogine, N. Trappeniers, and V. Mathot, J. Chem. Phys. 21, 559
38
J. S. Rowlinson and B. Widom, Molecular Theory of Capillarity ~Claren- ~1953!.
65
don, Oxford, 1982!. I. Prigogine, N. Trappeniers, and V. Mathot, Discuss. Faraday Soc. 15, 93
39
B. Smit, J. Chem. Phys. 96, 8639 ~1992!. ~1953!.
40
T. Dobashi, M. Nakata, and M. Kaneko, J. Chem. Phys. 72, 6685 ~1980!. 66
I. Prigogine, The Molecular Theory of Solutions ~Interscience, New York,
41
P. Padilla and S. Toxvaerd, Mol. Phys. 75, 1143 ~1992!. 1957!.
42
A. López Rodrı́guez, C. Vega, J. J. Freire, and S. Lago, Mol. Phys. 73, 691 67
J. Hijmans, Physica 27, 433 ~1960!.
~1991!. 68
M. Kurata and S. Isida, J. Chem. Phys. 23, 1126 ~1955!.
43
A. López Rodrı́guez, C. Vega, J. J. Freire, and S. Lago, Mol. Phys. 80, 69
P. J. Flory, J. Chem. Phys. 10, 51 ~1942!.
1565 ~1993!. 70
P. J. Flory, J. Chem. Phys. 12, 425 ~1944!.
44
J. Hautman and M. L. Klein, J. Chem. Phys. 91, 4994 ~1989!. 71
M. L. Huggins, Am. Chem. Soc. 64, 1712 ~1942!.
45
J. I. Siepmann and I. R. McDonald, Phys. Rev. Lett. 70, 453 ~1993!. 72
M. L. Huggins, J. Phys. Chem. 46, 151 ~1942!.

Downloaded from http://pubs.aip.org/aip/jcp/article-pdf/102/5/2126/8104268/2126_1_online.pdf

46
J. I. Siepmann and I. R. McDonald, Mol. Phys. 79, 457 ~1993!. 73
A. Kreglewski and B. J. Zwolinski, J. Phys. Chem. 65, 1050 ~1961!.
47 74
Note that the model which is studied in this work includes bond bending K. Nakanishi, M. Kurata, and M. Tamura, J. Chem. Eng. Data 5, 210
potentials and has a slightly different potential for the torsion than that in ~1960!.
the original work of Jorgensen et al. ~Ref. 19!. We found that the phase 75
P. J. Flory, R. A. Orwoll, and A. Vrij, J. Am. Chem. Soc. 86, 3507 ~1964!.
76
behavior is mainly determined by the nonbonded potential and is not very M. J. Anselme, M. Gude, and A. S. Teja, Fluid Phase Equilibria 57, 317
sensitive to the details of the intramolecular potentials. It can therefore be ~1990!.
77
expected that differences between the model used in this work and the Y.-J. Sheng, A. Z. Panagiotopoulos, S. K. Kumar, and I. Szleifer, Macro-
original alkane model of Jorgensen et al. are very small. molecules 27, 400 ~1994!.
48
P. Van der Ploeg and H. J. C. Berendsen, J. Chem. Phys. 76, 3271 ~1982!. 78
W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Nu-
49
P. Padilla and S. Toxvaerd, J. Chem. Phys. 94, 5650 ~1991!. merical Recipes: The Art of Scientific Computing ~Cambridge University,
50
Toxvaerd uses a slightly different bond-bending potential with equilibrium Cambridge, 1986!.
angles that depend on the molecule. Furthermore, in Toxvaerd’s original 79
S. Chynoweth and Y. Michopoulos, Mol. Phys. 81, 133 ~1994!.
model the bond length depends on chain length. These differences should 80
W. L. Jorgensen, J. Am. Chem. Soc. 103, 335 ~1981!.
have little effect on the phase behavior. 81
D. Brown and J. H. R. Clarke, Macromolecules 24, 2075 ~1991!.
51
N. G. Almarza, E. Enciso, and F. J. Bermejo, J. Chem. Phys. 96, 4625 82
D. Rigby and R.-J. Roe, J. Chem. Phys. 87, 7285 ~1987!.
~1992!. 83
D. J. Rosenthal and A. S. Teja, AIChE J. 35, 1829 ~1989!.
52
C. Tsonopoulos, AIChE J. 33, 2080 ~1987!. 84
B. D. Smith and R. Srivastava, Thermodynamics Data for Pure Com-
53
K. Magoulas and D. Tassios, Fluid Phase Equilibria 56, 119 ~1990!. pounds: Part A Hydrocarbons and Ketones ~Elsevier, Amsterdam, 1986!.

J. Chem. Phys., Vol. 102, No. 5, 1 February 1995