Scientific paper
Phase Diagram of the Lennard-Jones System of Particles from the Cell Model and Thermodynamic Perturbation Theory
Alan Bizjak, Tomaž Urbic* and Vojko Vlachy
Faculty of Chemistry and Chemical Technology, University of Ljubljana, A{ker~eva 5, 1000 Ljubljana, Slovenia * Corresponding author: E-mail: tomaz.urbic@fkkt.uni-lj.si
Received: 12-11-2008 Dedicated to Professor Josef Barthel on the occasion of his 80'' birthday
Abstract
The thermodynamic perturbation approach and the cell theory are used to determine the complete phase diagram of a system of particles interacting via the Lennard-Jones potential. The Barker-Henderson perturbation theory (J. A. Barker, D. Henderson, J. Chem. Phys. 1967, 47, 4714-4721) is invoked to calculate the liquid-vapour line, while the cell model approach is utilized to evaluate the solid-vapour and solid-liquid equilibria lines. The resulting phase diagram along with the triple point and the critical point conditions are determined. Although the liquid-vapour line is predicted quite well, the triple point parameters, T t*and Pt*are in relatively poor agreement with the available computer simulation data. This finding, together with our other calculations, suggest that a simple cell theory may not be adequate for characterizing the solid-fluid equilibrium in systems involving strongly correlated molecules.
Keywords:
1. Introduction
Over the last few decades there has been a fair amount of research into understanding and characterizing the vapour-liquid equilibrium in different types of fluids. In contrast, much less is understood about the solid-liquid equilibrium, which has inhibited construction of complete phase diagrams. Among the various model fluids, the system of particles interacting via the Lennard-Jones 6-12 potential has been of interest historically and plays a special role.1 This simple two-parameter potential is important perse, since it describes the interactions in rare gas-es2-4 quite well. Due to its intuitive simplicity the Lennard-Jones potential is often used as a starting point for construction of more elaborate fluid models (see, for example, reference [5]) and/or to mimic hydrophobic interactions.
The celebrated Lennard-Jones 6-12 potential reads1
(1)
where e is the depth of the potential well and a is the size parameter. These two parameters can be determined from experiments (see, for example, Tables 1 and 11 of Reference6). The first term in the square brackets has the right functional form for the dispersion interaction for large r values. The r12 form for the core repulsion is chosen for convenience and does not represent an optimal choice. For r = a the attractive and repulsive part are exactly balanced. The Lennard-Jones fluid has been extensively studied in the literature using different approaches such as the thermodynamic perturbation theory,37 integral equation theories,8-10 numerical machine simulations,1112 density functional theories,13 and cell theories.14 From the weight of evidence it is clear that the Lennard-Jones system remains one of the most interesting systems to test new theories. 1t is interesting that most of the theoretical studies and simulations to date have focused on the liquid-vapour transition and the critical point determination with the most recent paper in this regard being that by Betacourt-Cardenas and coworkers.15 Much less theoretical effort has been directed toward the accurate determination of the solid-liquid coexistence line.4,6 A complete phase diagram of the Lennard-Jones system of particles,
obtained by machine simulations was presented by Lofti et al.,16 and Agarwal and Kofke.17
Water is, without any doubt, a universal and the most important fluid. Its properties are thought to arise from the ability of water to form tetrahedrally coordinated hydrogen bonds and are subject of many studies (for recent reviews see references18-20). One approach uses atomistic, detailed simulation models, which includes a number of variables to describe van der Waals and Coulomb interactions, and hydrogen bonding.21 Another approach is based on simpler models having less structural details and consequently requiring much less computation.22-28 Recently, a simple model of water, where each molecule is a Lennard-Jones sphere having four hydrogen-bonding »arms« oriented tetrahedrally, has been proposed.5 The Monte Carlo method was used to determine the basic thermodynamic properties of this water-like fluid. Despite the simplicity of this model, the computer simulation of the complete phase diagram would be very time consuming and consequently alternative routes need to be explored.
There are two main methodologies to calculate the phase diagram of a model fluid. One uses an appropriate computer simulation method (see, for example, reference 29) Alternative approaches are more analytical and hence generally less time consuming, however, they contain various statistical-mechanical approximations. Among the most popular theoretical methods are the density functional30,31 and the cell-model theories,6,32,33 combined with the appropriate equation of state. In order to calculate the fluid-solid part of the phase diagram it is often necessary to apply two different theories: one to describe the fluid phase, and another to characterize the equilibrium solid. In this contribution we will explore two theoretical approaches: i) the Barker-Henderson theory will be used2,34 to determine the properties of the fluid phase, while ii) the thermodynamic parameters of the solid phase will be calculated by the cell theory.6,32,33 Both of these techniques, especially the Barker-Henderson perturbation theory, have been used before to determine parts of the phase diagram of Lennard-Jones fluid but, but it is unclear how accurate they are when working in tandem to obtain the solid-fluid part of the phase diagram. A calculation along these lines was first proposed by Henderson and Barker4 who estimated the triple point properties of Argon, but did not present a full phase diagram of the system. In passing we also note an earlier study by Cottin and coworkers,6 who used a semi-empirical equation of state (instead of the perturbation approach) to describe the fluid phase.
The interest in this work is therefore to test the above mentioned combination of theories on the Lennard-Jones system of particles. If the theory is working we are planning to apply it to systems involving strongly correlated molecules5. Phase diagram of the Lennard-Jones fluid is known from simulations16,17 as also from experi-
ments for systems behaving like Lennard-Jones ones.35-38 The remainder of this paper is organized as follows. In Section 11 we describe how to use the cell theory to calculate the chemical potential of the solid phase, and the thermodynamic perturbation approach to obtain the same information for the liquid phase. 1n Section 111 we will present the phase diagram of Lennard-Jones particles, and finally, in Section 1V some conclusions are drawn from this study.
2. Methods of Calculation
Methods of choice were the cell theory6,32,33 for description of the solid phase, and thermodynamic perturbation theory to characterize the liquid and vapour2 phases. 1n order to calculate the phase equilibrium curves we need to evaluate the Helmholtz free energy of each phase as a function of pressure and temperature.
2. 1. Solid Phase: Cell Theory
We start from the perfect lattice; it is well-known that Lennard-Jones (LJ) particles crystallize39 in the form face centered cubic (FCC) crystal lattice structure. The shape of the Wigner-Seitz cell for the face centered cubic crystal is rhombic dodecahedron.40 We chose one atom as the central one. The interaction energy of this central atom with all the other atoms of the solid phase is given by6
[/(r) = C/o + A[/(r)
(2)
where U0 is the lattice energy of a central atom in equilibrium configuration and can be calculated as6

(3)
1ndex 0 denotes the central particle and j all the other particles in such crystal. Further, AU(r) is the change in interaction energy of a central atom when displaced for r from its equilibrium position; U(r) can be calculated as
= ^ [f7,,(|rjo - rl) - Uu{rjo)] ■
(4)
1n this approximation we assume that all the neighboring atoms are fixed at their lattice points. The Helmholtz free energy A(N; V; T) of the solid within the framework of the cell theory is given by the expression41
(5)
where ß = 1 / (kBT) and q is the configurational partition function evaluated as

(6)
where r defines the position of a central atom in the Wigner-Seitz cell and must be evaluated numerically. In this study we utilized the Monte Carlo integration procedure to evaluate the integral of Eq. 6.32,33
2. 3. Fluid Phase: Thermodynamic Perturbation Theory
The key quantity here is the Helmholtz free energy for the system of interest. We have obtained this quantity using the information about the hard-sphere system of particles; in other words, the free energy of hard spheres AHlS, was taken as the reference quantity. The Barker-Henderson perturbation theory2,3,34 was than utilized to calculate the excess free energy ALJ of the Lennard-Jones fluid.
The excess Helmholtz free energy per particle is accordingly written as the sum of the hard-sphere contribution AHXS and the first-order correction due to the Lennard-Jones interaction potential
(7)
Uu{r)9Hs{r,ri)dr
where a^j is the Lennard-Jones size parameter, with
n = ß^f^V
(8)
being the packing fraction of the reference hard-sphere system of particles with diameter d. To calculate the hard-sphere term of the Helmholtz free energy, we integrated the Carnahan-Starling42 expression for the equation of state
ßA^it, = NI
1 z - I
dt]
where z is
(9)
(10)
The diameter of the hard-sphere reference system, d, was determined through
(11)
Pair distribution function of the hard-sphere system, gHS(r,n), was approximated by the expression of Gonzalez et al.43
2. 3. Construction of the Phase Diagram
Once the Helmholtz free energy is known, we can calculate pressure and chemical potential using the standard thermodynamic equations

(12)
Derivatives were made numerically. For a one-component system with two phases (here denoted by 1 and 2), the conditions of thermodynamic equilibrium read44:
(13)
The three equations must be satisfied simultaneously. In the calculated chemical potential we have to include both ideal and excess parts of the quantity. The condition contained in Eq. 13 is equivalent to construction of common tangents to the solid and fluid free energy curves, given as a function of the reduced volume. Figure 1 shows reduced as a function of the reduced f at 1* = 0.7.
Figure 1: The reduced chemical potential as a function of the reduced pressure at f = 0.7. Continuous lines belong to the fluid phase (obtained by the TPT theory) and the dashed line to the solid phase (cell theory). The values of ß* and P" where the dashed and solid lines cross correspond to the fluid-solid equilibrium and where solid lines cross to liquid-liquid. For the fluid-solid transition ß* = -4.13 and P* = 0.0022; for the fluid-fluid transition ß* = -3.81 and P* = 0.0032.
Vapour-liquid transitions were calculated by the Maxwell construction.44 The isotherms under the critical temperature are shaped as shown in Figure 2. The dotted curve, from point A to minima corresponds to the overheated liquid, and from maximum to point B to the under-
cooled gas. The dotted curve between the minimum and maximum is unphysical since this cannot be seen in reality; such a fluid would be characterized by a negative compressibility. The end points of the two-phase region (A;B) can be predicted from the condition that the shaded area above the tie-line must be equal to that below it. This condition arises from the equality of the free energies of liquid and vapour phases.
3. Results
All the results reported here are given in reduced units in terms of the Lennard-Jones size parameter, aLJ, and the depth eLJ, of the Lennard-Jones potential well,
viz., r = kBT/leLJl, P* = Pa'u, P* = Pa'LJ/lfLJl, and M* = ß/\eLJ\. First in Figure 3 we show the T vs p* dependence. The dashed line represent the vapour-liquid equilibrium as calculated from
Figure 2: Maxwell's construction of the coexisting points (horizontal line from point A to point B) for the liquid-gas phase equilibrium at T* = 1.2.
P
Figure 3: Phase diagram of the Lennard-Jones fluid: Reduced temperature as the function of the reduced density. The dashed line corresponds to the liquid-vapour coexistence line and the solid ones to solid-liquid and solid-vapour coexistence.
the Barker-Henderson thermodynamic perturbation theory. The highest point on the dashed line corresponds to the critical temperature. We have estimated the critical temperature to be T* « 1.35. Computer simulations17 yield T* « 1.31, while experimental data for argon are a bit lower T* « 1.26. Since our theoretical approach is analogous to those in the literature, it is not surprising that a good agreement with old2 and recent evaluations15 of T* was obtained. Critical density, p** , as calculated by the thermodynamic perturbation theory, is comparable with that obtained from simulations and experiment (cf Table 1). Continuous lines of Figure 3 show the solid-fluid boundaries.
Figure 4: Phase diagram of the Lennard-Jones fluid: P* - T dependence. Solid lines denote our calculation and the dashed lines
the computer simulation results.16,17
The next figure, Figure 4, shows the pressure-temperature dependence. Here the agreement between the theoretical prediction (solid line) on one side, and the simulation and experiment (dashed lines) on the other is less favorable. The temperature at the point where all three phase co-exist (triple point) was determined to be T* « 0.80. For argon the experimentally determined temperature at triple point is T*« 0.70. The computer simulation17 yields the value T*« 0.68. Additional detailed comparative results can be seen in Table 1. Clearly the error in the triple point calculation is considerably bigger than the error in the critical point evaluation. One reason for the poor agreement seems to lie in the approximations inherent to
Tabele I: Comparison between the our theoretical results, simula-tions16,17 and experimental data.45 The Lennard-Jones parameters for noble gases were taken from Ref.46
	T*	Pc*	Pc	T*	P*
theory	1.35	0.16	0.30	0.80	0.0086
simulation	1.31	0.12	0.30	0.68	0.001
neon	1.27	0.12	0.31	0.70	0.0019
argon	1.26	0.12	0.32	0.70	0.0016
krypton	1.22	0.11	0.30	0.68	0.0015
xenon	1.31	0.13	0.35	0.73	0.0018
the cell theory, where the movements of surrounding atoms on crystal lattice were neglected. The central particle moves in the average field produced by all other atoms, fixed at their lattice sites. This approximation may prove inadequate for more strongly correlated system of particles than that occurring in a Lennard-Jones fluid.
4. Conclusions
The principal achievement of this paper has been a calculation of the phase diagram of the Lennard-Jones system of particles utilizing the thermodynamic perturbation theory in conjunction with the cell theory. Our motivation was to test the performance of the combination of these theories on the well known Lennard-Jones system, before applying them to the molecular systems with directional forces. The predicted phase diagram is qualitatively similar to the experimental and simulation results. More precisely the approach used here predicts the liquid-vapour part of the phase diagram quite well (this is of course known from before), but is considerably less successful in describing the solid-gas and solid liquid branches of the phase diagram. The origin of the disagreement can be traced into the approximations inherent to the cell model approach. This indicates that some modifications are needed before the cell theory can be successfully applied to more complicated models of fluid. In the cell theory of solids used in the present work the central particle is restrained to move in the neighborhood of its lattice site in an average potential field created by all the other particles. The latter particles are fixed at their lattice sites. Such a mean field description may not be suffcient to capture correctly the behavior of strongly correlated particles, such as water molecules, which interact via highly directional forces. In fact, our preliminary calculations of the phase diagram of a water-like fluid5 indicate that the critical point of such a fluid may fall within the solid region. One way of solving this dilemma is to try and improve the cell theory by introducing the possibility of cooperative motions of two or more particles in the system47. An alternative, albeit more time consuming, way would be to use the Monte Carlo method to calculate the solid-fluid equilibrium lines.
5. Acknowledgements
This work was supported by the Slovenian Research Agency (Program »Physical Chemistry«, P1 0103-0201). Authors wish to thank Professor Lutful B. Bhuiyan for critical reading of the manuscript.
6. References
1.	J. E. Lennard-Jones, Proc. Phys. Soc. 1931, 43, 461-482.
2.	J. A. Barker, D. Henderson, Rev. Mod. Phys. 1976, 48, 587671.
3.	J. A. Barker and D. Henderson, J. Chem. Phys. 1967, 47, 4714-4721.
4.	D. Henderson, J. A. Barker, Molec. Phys. 1968, 14, 587-589.
5.	A. Bizjak, T. Urbic, V. Vlachy, K. A. Dill, Acta Chim. Slov. 2007, 54, 532-537.
6.	X. Cottin, P. A. Monson, J. Chem. Phys. 1996, 105, 1002210029.
7.	J. D. Weeks, D. Chandler, H. C. Andersen, J. Chem. Phys. 1971, 54, 5237-5247.
8.	R. O. Watts, J. Chem. Phys. 1968, 48, 50-55.
9.	A. M. Berezhkovskii, N. M. Kuznetsov, I. V. Fryazinov, JAMT 1972, 13, 227-233.
10.	J. A. Anta, E. Lomba, C. Martin, M. Lombardero, F. Lado, Mol. Phys. 1995, 84, 743-745.
11.	L. Verlet, Phys. Rev. 1967, 159, 98-103.
12.	J. K. Johnson, J. A. Zollweg, K. E. Gubbins, Mol. Phys. 1993, 78, 591-618.
13.	Y. Tang, J. Wu, J. Chem. Phys. 2003, 119, 7388-7397.
14.	J. E. Magee, N. B. Wilding, Mol. Phys. 2002, 100, 16411644.
15.	F. F. Betacourt-Cardenas, L. A. Galicia-Luna, S. I. Sandler, Fluid Phase Equilibria 2008, 264, 174-183.
16.	A. Lotfi, J. Vrabec, J. Fischer, Mol. Phys. 1992, 76, 13191333.
17.	R. Agrawal, D. Kofke, Mol. Phys. 1995, 85, 43-59.
18.	L. R. Pratt, Annu. Rev. Phys. Chem., 2002, 53, 409-436.
19.	W. Blokzijl, J. B. F. N. Engberts, Angew. Chem. Int. Ed. Engl., 1993, 32, 1545-1579.
20.	A. Ben-Naim, Biophys. Chem. 2003, 105, 183-193.
21.	B. Guillot, J. Mol. Liq. 2002, 101, 219-260.
22.	I. Nezbeda, J. Kolafa, Yu. V. Kalyuzhnyi, Mol. Phys. 1989, 68, 143-160.
23.	I. Nezbeda, J. Mol. Liq. 1997, 73/74, 317-336.
24.	K. A. T. Silverstein, A. D. J. Haymet, K. A. Dill, J. Am. Chem. Soc. 1998, 83, 150-151.
25.	T. Urbic, V. Vlachy, Yu. V. Kalyuzhnyi, N. T. Southall, K. A. Dill, J. Chem. Phys. 2000, 112, 2843-2848.
26.	T. Urbic, V. Vlachy, Yu. V. Kalyuzhnyi, N. T. Southall, K. A. Dill, J. Chem. Phys. 2002, 116, 723-729.
27.	T. Urbic, V. Vlachy, Yu. V. Kalyuzhnyi, N. T. Southall, K. A. Dill, J. Chem. Phys. 2003, 118, 5516-5525.
28.	K. A. Dill, T. M. Truskett, V. Vlachy, B. Hribar-Lee, Annu. Rev. Biophys. Biomol. Struct. 2005, 34, 173-199.
29.	D. Frenkel, B. Smit, Understanding Molecular Simulation, Academic Press, San Diego, 2002.
30.	Y. Singh, Phys. Rep. 1991, 207, 351-444.
31.	S. W. Rick, A. D. J. Haymet, J. Phys Chem. 1990, 94, 52125220.
32.	E. P. A. Paras, C. Vega, P. A. Monson, Mol. Phys. 1992, 77, 803-821.
33.	C. Vega, P. A. Monson, Mol. Phys. 1995, 85, 413-421.
34.	T. Boublik, 1. Nezbeda, K. Hlavaty, Statistical Thermodynamics of Simple Liquids and Their Mixtures, Elsevier, Amsterdam 1980.
35.	A. Michels, H. Wijker, H. Wijker, Physica 1949, 15, 627633.
36.	A. Michels, J. M. Levelt, W. de Graff, Physica 1958, 24, 659- 671.
37.	A. van 1tterbeek, O. Verbeke, Physica 1960, 26, 931-938.
38.	A. van 1tterbeek, J. de Boelpaep, O. Verbeke, F. Theeuwes, K. Staes, Physica 1964, 30, 2119-2122.
39.	D. G. Henshaw, Phys. Rev. 1958, 111, 1470-1475.
40.	N. W. Ashcroft, N. D. Mermin, Solid State Physics, Thomson Learning, Australia 1976.
41.	C. Vega, P. A. Monson, J. Chem. Phys. 1998, 109, 99389949.
42.	N. F. Carnahan, K. E. Starling, J. Chem. Phys. 1969, 51, 635-636.
43.	D. J. Gonzalez, L. E. Gonzalez, M. Silbert, Mol. Phys. 1991, 74, 613-627.
44.	R. Kubo, Thermodynamics; an advanced course with problems and solutions, North-Holland, Amsterdam 1968.
45.	D. R. Lide, CRC Handbook of Chemistry and Physics, CRC Press, Boca Raton, New York, London, Tokyo 1999.
46.	T. L. Hill, An 1ntroduction to Statistical Thermodynamics, Addison-Wesley Publ., London, 1960.
47.	J. A. Barker, Lattice Theories of the Liquid State, Pergamon, New York, 1963.
Povzetek
1zracunali smo fazni diagram za sistem, ki ga sestavljajo Lennard-Jonesove kroglice. Lastnosti trdne faze smo opisali s celično teorijo, parno in tekočo fazo pa s pomočjo termodinamske perturbacijsko teorije. V slednjem primeru smo uporabili Barker-Hendersonov približek (J. A. Barker, D. Henderson, J. Chem. Phys. 1967, 47, 4714-4721). 1zračunana ravnotežna črta za prehod kapljevina-plin se dobro ujema s simulacijo in z ustreznimi poskusi, medtem ko je ujemanje za fazni prehod tekoče-trdno precej slabše. Slabo ujemanje med teorijo in simulacijami je najbolj opazno v okolici trojne točke (T*in Pt*). Smatramo, da so za to krive poenostavitve celične teorije. Na osnovi teh in tudi drugih, še neobjavljenih, rezultatov je moč sklepati, da celična teorija, v obliki kot je uporabljena tukaj, verjetno ni primerna za opis sistemov močno koreliranih molekul.