Generalized Coupling Parameter Expansion: Application to Square well and LennardJones fluids
Abstract
The coupling parameter expansion in thermodynamic perturbation theory of simple fluids is generalized to include the derivatives of bridge function with respect to coupling parameter. We applied seventh order version of the theory to SquareWell (SW) and LennardJones(LJ)fluids using Sarkisov Bridge function. In both cases, the theory reproduced the radial distribution functions obtained from integral equation theory (IET) and simulations with good accuracy. Also, the method worked inside the liquidvapor coexistence region where the IETs are known to fail. In the case of SW fluids, the use of CarnahanStarling expression for Helmholtz free energy density of HardSphere reference system has improved the liquidvapor phase diagram(LVPD) over that obtained from IET with the same bridge function. The derivatives of the bridge function are seen to have significant effect on the liquid part of the LVPD. For extremely narrow SW fliuds, we found that the third order theory is more accurate than the higher order versions. However, considering the convergence of the perturbation series, we concluded that the accuracy of the third order version is a spurious result. We also obtained the surface tension for SW fluids of various ranges. Results of present theory and simulations are in good agreement. In the case of LJ fluids, the equation of state obtained from the present method matched with that obtained from IET with negligible deviation. We also obtained LVPD of LJ fluid from virial and energy routes and found that there is slight inconsistency between the two routes. The applications lead to the following conclusions. In cases where reference system properties are known accurately, the present method gives results which are very much improved over those obtained from the IET with the same bridge function. In cases where reference system data is not available, the method serves as an alternative way of solving the Ornstein Zernike equation with a given closure relation with the advantage that solution can be obtained throughout the phase diagram with a proper choice of the reference system.
pacs:
61.20.Gy, 64.60.l, 64.70.pI Introduction
Describing accurately the thermodynamics and structure of fluids with short ranged interparticle potentials has been a long standing problem in liquid state theory. Various methods have been developed to attack the problem. They can be divided into two categories. One category is based on integral equation theory(IET) and the other is based on thermodynamic perturbation theory(TPT) (1). Methods based on IET take into account correlations between particles and aim at accurate description of structural properties of liquids i.e., the radial distribution function (RDF), direct correlation function (DCF) etc. The IETs involve solving the OrnzteinZernike equation(OZE) coupled with a closure relation. Thermodynamic properties like internal energy and pressure are obtained as integrals of RDF. But inside the two phase region, most of the closures do not have a solution(2). Apart from this, the IET method suffers from lack of thermodynamic consistency between the energy, virial and compressibility routes. Various closure relations were proposed to minimize the inconsistency(3). Also closures with adjustable parameters have been proposed to enforce thermodynamic consistency by adjusting them. The requirement of consistency makes the computations complicated. Despite the improvements, accurate description of thermodynamic and structural properties of simple liquids with narrow potentials is still lacking.
In earlier formulations of TPT, the structure and thermodynamic properties of reference system are assumed to be known and the thermodynamic properties of actual system are obtained as a small perturbation over the reference system. The structure of the reference system and the actual system are assumed to be the same. Within this approach, it is not possible to go beyond second order term in the perturbation series as they require higher order correlation functions(4). Zhou(5) developed a method in which the Helmholtz free energy of the system is written as a series over coupling parameter . This method requires derivatives of the radial distribution function (RDF) at . Zhou calculated the derivatives using finitedifference method and could obtain terms in the perturbation series up to sixth order(6). However, as one goes to higher and higher orders, errors due to numerical differentiation creep in and the data sets of derivatives of RDF require additional smoothing to remove fluctuations. So it is difficult to get derivatives of RDF to arbitrary order. On the other hand, for very shortranged interparticle potentials, the convergence of the free energy series becomes very slow and requires higher order terms to be evaluated.
Recently, we addressed the problem in a more fundamental way combining the ideas of IET and TPT. We assumed that the RDF and the DCF of a system can be written as Taylor series in coupling parameter at . We showed that the terms in Taylor series of both RDF and DCF (i.e., their derivatives) can be obtained by solving a simple set of equations obtained by differentiating the OZE and the generalized closure . The formula for free energy given by Zhou can be derived within our formalism apart from the RDF and DCF of the actual system. Within this approach we obtained DCF, RDF and LVPDs for Square Well (SW) fluids of various ranges up to seventh order. By comparing the third, fifth and seventh order versions of our theory, we concluded that the perturbation series for Helmholtz free energy has practically converged by seventh order. However, the series for RDF and DCF didn’t converge for low densities and low temperatures for narrow square well fluids by seventh order. Even though the Helmholtz free energy series has converged, we found that there was significant deviation of the obtained LVPD using seventh order TPT from simulation results for narrow SW fluids. The deviation and slow convergence could be because of two reasons: Firstly, we neglected the derivatives of bridge function in our calculations and secondly the error in the bridge function we have chosen for calculation.
Also, in our previous work, we were able to obtain the RDF and DCF inside the spinodal and close to the critical region without any problems of convergence of the iteration scheme even though our numerical scheme is very simple. This raises questions about the existence of the solution to OZE inside the spinodal region. According to Sarkisov(7) , solutions to OZE exist in the spinodal region also. However, the IET methods available are not convergent in that region. The RDF and DCF in the twophase region find applications in the density functional theory of fluids (8). One example is the calculation of surface tension. This calculation requires the DCF in the two phase region. As the IETs doesn’t have solution in some part of the LVPD, earlier calculations of surface tension were done using an interpolation formula given by Ebner (9) to obtain inside the twophase region.
In present paper, we address these issues using a generalized version of our method with application to SW fluids. We use the bridge function proposed by Sarkisov(7) with slight modification for SW fluids(11) and include its derivatives in the calculation. We obtain RDFs, isotherms and LVPDs for SW fluids of various ranges and compare with those obtained from IET and simulations. We also calculate the surface tension for SW fluids of ranges , and using the expression obtained from the square gradient functional for Helmholtz free energy and compare with simulation results. As an application of our theory to nonhard sphere reference systems, we apply it to LennardJones (LJ) fluid. We use the Sarkisov bridge function for both reference system and perturbation part. The RDFs, isotherms and LVPD for LJ fluid are obtained and compared with those obtained from IET and simulations wherever available. The paper is organized as follows. In Section II we discuss the method briefly. In Section III we apply the theory to SW fluids to LJ fluids and the results are analyzed. The paper is concluded in Section IV.
Ii Theory
Basic theory is the same as explained in Reference[12] except for a generalization to include derivatives of bridge function in the calculation. Consider a fictitious system at temperature and density with interaction potential given by
(1) 
is a coupling parameter. When , the fictitious system becomes reference system with interaction potential . A nonzero will add a perturbation to as shown in the above equation. corresponds to the potential of the system under consideration. We postulated that the RDF and DCF of the system with potential can be written as a McLaurin series in , that is,
(2) 
(3) 
Hereafter we shall denote partial derivative of any function as .
The radial distribution function can be written as
(4) 
where is the indirect correlation function defined as and is the total correlation function of the fictitious fluid. is called the potential of mean force. The expression for in Eq.(4) is exact(1). called the bridge function (3). Since , and hence , as well as are expanded in a series in , the correlation function is also a series in . The order coefficient in its series is given by . Generally, is chosen as a function of . Several approximations to in terms of are available (3) in literature.
The general expression for order derivative of from Eq.(4) is
(5) 
where is the binomial coefficient. The derivatives are given by
(6) 
where is the Kronecker delta. The derivatives can be computed using Eq.(5) in a recursive manner, using either from initial guess or previous iteration.
Since
(7) 
its order derivative is
(8) 
To get another set of relations between and , we consider the Ornstein Zernike Equation (OZE) in Fourier space:
(9) 
where and are the Fourier transforms of and , respectively.
Also, structure factor is defined as
(10) 
Therefore, using above equation in OZE i.e., Eq.(9) we get
(11) 
Now, differentiating Eq.(10) times, we get
(12) 
Using Eq.(11) the derivatives of are expressed as
(13) 
Thus in order to obtain up to order term in Taylor series expansion of both and , coupled linear equations in realspace and coupled linear equations in Fourier space have to be solved simultaneously with .
For example, to obtain Taylor series expansion up to second term, the set of four equations given by
(14)  
(15)  
(16)  
(17) 
have to be solved (where is ).
The CPE for Helmholtz free energy density of a homogeneous fluid at a given temperature is given by (1)
(18) 
where is the free energy density of the reference system. For notational simplicity we do not show temperature dependence of explicitly hereafter. Substituting Eq.(3) in Eq.(18) and integrating over , we get
(19) 
Here we have used the shortened notation for the derivatives which is readily obtained as . Thus the method provides the DCF, RDF as well as the free energy density. The Taylor series of RDF truncated up to second order gives third order CPE for . The numerical procedure we used is as explained in our previous paper(12) and is not repeated here. Eq.(19) has been first derived by Zhou(5) in a different way.
Iii Applications
iii.1 Square Well Fluids
Seventh order version of above mentioned theory is applied to SW fluids. We used the closure proposed by Sarkisov with slight modification by Mendoub(11).
(20) 
where
(21)  
and is the Hard sphere diameter. We obtained the RDF and DCF of the reference system using the same bridge function by solving the OZE using a similar procedure as explained above to maintain consistency.
Reduced units (, where is the well depth and is hard sphere diameter) are used throughout the paper. For convenience, we denote , by , respectively and , by , respectively. In Fig.(1) we compare obtained using order version of our TPT and that obtained through IET(11) for SW fluid of range for densities and and temperature . Clearly, there is negligible difference between results obtained using present method and those obtained from IET. Also, the agreement with simulation results is good. We observed that except for very low temperatures and low densities, results of fifth order and seventh order TPTs have negligible deviation showing that the Taylor series has converged. Convergence of the iteration scheme has been good in the whole phase diagram. obtained by our method for SW fluid of range in the spinodal region at and is shown in Fig.(1). Above observations show that the present method may be viewed as an alternative way of solving the OZE.
In Fig.(2) we give plots of LVPD for SW fluids of ranges and respectively obtained using seventh order version of TPT with and without including the derivatives of bridge function in the calculations. We used the CarnahanStarling(CS) expression for Helmholtz free energy density of the reference system. Our results are compared with simulations, those obtained from IET(11) and our previous results using MalijevskyLabik brdige function neglecting the derivatives of bridge function. Results of the present calculations are very much improved in the liquid part of the LVPD over those of Mendoub (11). The main reason for this improvement is the use of the CS expression for free energy density. A comparison of the LVPDs obtained by including derivatives of the bridge function in the calculation and those obtained by neglecting them shows that the correction due to the derivatives becomes important close to the critical region and along the liquid side of the phase diagram. Neglecting them leads to shifting of the critical point to high density region.
In Fig.(3), pressure isotherms obtained using , and order TPT for SW fluids of extremely short ranges and are plotted for temperatures and respectively. The figures depict that the order results are much closer to the simulation results. However, a comparison of the three versions shows that the series is indeed convergent at each density point. But it is converging to a point away from the simulation results. This is because of the inaccuracy of the bridge function for very short ranged fluids. Thus the apparent closeness of the third order TPT results to the simulation results is spurious.
Possibility of obtaining in two phase region opens up applications in density functional theory(DFT) of liquids. One simple application is calculation of surface tension of liquids. We use the obtained from seventh order version of TPT to obtain surface tension for SW fluids. Formula for surface tension is obtained by Yang et. al. (13) from squaregradient functional for Helmholtz free energy of inhomogeneous fluids. A brief derivation is as follows: The squaregradient approximation for Helmholtz free energy functional of an inhomogeneous liquid occupying volume at temperature is given by
(22) 
where is Helmholtz free energy density of homogeneous liquid. The second term is the effect of inhomogeneity. is the number density in an infinitesimal volume around ., and are all functions of even though the dependence is not shown explicitly.
The coefficient of gradient term also referred as influence parameter is
(23) 
We assume that liquidvapor interface is flat and that axis is normal to the interface pointing out into the vapor from the liquid. In such a case, Eq.(22) becomes
(24) 
where is the surface area. Grand free energy of the system is given as
(25) 
where is the chemical potential of the system and is the total number of particles
Minimizing we get
(26) 
where is the grand free energy density.
Integrating above Eq.(26) with boundary conditions
(27) 
gives
(28) 
where and are coexisting liquid and vapor densities for temperature under consideration. Surface tension can be calculated using the formula
(29) 
Using Eq.(28) in above equation gives
(30) 
In the above equations may be replaced by also as both have same value.
Above explained formalism is applied to SW fluids of ranges , and with obtained from seventh order TPT as explained above. Fig.(4) shows as a function of for different temperatures for SW fluid of range . From Fig.(4) it can be seen that becomes negative at high densities which is unphysical. This might be an artifact of the approximate bridge function used. However, within the binodal, where is required for the calculation of surface tension(S), it is positive. Surface tensions obtained from our calculation are plotted as a function of in Fig.(5). Our results compare well with simulations except close to the critical temperatures.
iii.2 Lennard Jones Fluid
The theoretical formalism explained in section II is general and can be applied to nonhard sphere reference systems also. As an example, we applied the theory to Lennard Jones(LJ) fluid. The Lennard Jones potential was split into reference and perturbation parts according to the Weeks Chandler and Anderson (WCA)(14) method. In this case, we used for , the expression proposed by Sarkisov (7) which is same as Eq.(20) with
(31)  
where is the minimum energy point of the LennardJones potential. The and of the reference system are obtained by solving the OZE with the same bridge function. Thus, in effect, we solved the OZE for LJ fluid with the above bridge function using the perturbation method. In Fig.(6), we compare the of LJ fluid obtained from seventh order TPT with simulation results. There is excellent agreement between the theory and simulation results for the cases shown except for a slight deviation in the case of lowest temperature. In Fig.(7), Equation of State(EOS) of LJ fluid for various temperatures obtained using seventh order TPT is compared with those obtained by solving OZE by Sarkisov(7). Pressure(P) is obtained using the virial formula given by
(32) 
Values obtained by our method matched with negligible deviation from those given by Sarkisov. In this case also, using our method we could obtain the pressure for all density points in the phase diagram, as at any point in the phase diagram could be calculated. Again, even close to the critical region, the same numerical procedure was sufficient whereas IETs face numerical convergence problems in this region. In Fig.(8), LVPD of LJ fluid is shown and compared with simulation results. We plotted the LVPD obtained in two ways. One is obtained by Maxwell construction of pressure isotherm obtained using Eq.(32) which is the so called virial route. Alternatively, LVPD is also obtained using the energy route. It is obtained as follows: Pressure isotherm of the reference fluid is obtained using Eq.(32) with . From this, Helmholtz free energy density of the reference fluid is obtained from the formula below
(33) 
Once the reference free energy is obtained, method described in section II can be applied to get the free energy density of the required system. The pressure isotherm is obtained by differentiating the Helmholtz free energy volume of the system. Maxwell construction is done to get the coexistence points. From the figure, it can be seen that the LVPD obtained from both the routes differ slightly along the liquid part of coexistence curve. Also, there is some deviation of both LVPDs from the simulation results. This is due to the bridge function used. Even though Sarkisov bridge function is supposed to be quiet accurate, there is some inconsistency between various thermodynamic routes. Imposition of the thermodynamic consistency between various routes as a constraint would solve the problem and may improve the accuracy of the results. This can be done in a straightforward way within the theory presented in the paper.
Iv Conclusion
We have generalized the TPT of fluids based on coupling parameter expansion to include the derivatives of the bridge function. The applications presented in the paper let us conclude the following points.(i) Inclusion of derivatives of the bridge function in the TPT has a significant effect on structural and thermodynamic properties of liquids. (ii) If accurate reference system data is available, it can be used in the TPT to improve the accuracy of the results over those of IET. (iii) The present method may be viewed as an alternative way of solving the OZE with the advantage that solution can be obtained in the whole phase diagram.
V Acknowledgements
I am indebted to Dr. S.V. G. Menon for introducing me to the subject when he was in Theoretical Physics Division, B.A.R.C. I thank Dr. D. M. Gaitonde for useful comments and discussions.
References
 J. P. Hansen, and I. R. McDonald, Theory of Simple Liquids, Academic Press, London, 2006.
 P. M. W. Cornelisse, C. J. Peters and J. de Swaan Arons, J. Chem. Phys. 106, 9820 (1997) and references therein.

J. M. Bomont, Recent Advances in the Field of Integral Equation Theories,
Adv. in Chem. Phys. vol. 139, John Wiley and Sons  J. A. Barker and D. henderson, J. Chem. Phys. 47, 2856(1967); Rev. Mod. Phys., 48, 587, 1976.
 S. Zhou, Phys. Rev. E 74, 031119(2006)
 S. Zhou, Phys. Rev. E 77, 041110(2008),J. Chem. Phys. 130, 054103 (2009)
 G. Sarkisov J. Chem. Phys. 114, 9496 (2001)
 R. Evans , Adv. in Phys. 28, 143(1979)
 C. Ebner, W. F. Saam, and D. Stroud, Phys. Rev. A 14, 2264 (1976).
 N. F. Carnahan and K. E. Starling, J. Chem. Phys. 51, 635 (1969)
 E. B. El Mendoub, J. F. Wax, I. Charpentier, N. Jakse, Mol. Phys. 106, 2667(2008)
 A. Sai Venkata Ramana and S. V. G. Menon, Phys. Rev. E 87, 022101 (2013)
 Arthur J. M. Yang, Paul D. Fleming and Julian H. Gibbs, J. Chem. Phys. 63, 3732(1976)
 John D. Weeks, David Chandler and Hans C. Andersen, J. Chem. Phys. 54, 5237(1971)
 J. Largo, J. R. Solana, S. B. Yuste, and A. Santos, J. Chem. Phys. 122, 084510 (2005)
 L. Vega, E. de Miguel, L. F. Rull, G. Jackson and I. A. McLure, J. Chem. Phys. 96, 2296 (1992)
 Fernando D. Rio, E. Avalos, R. Espindola, L. F. Rull, G. Jackson and S. Lago, Mol. Phys. 100, 2531 (2002)
 S. Zhou, J.R. Solana Phys. Chem. Chem. Phys. 11, 11528(2009)
 Jayant K. Singh, David A. Kofke and Jeffrey R. Errington, J. Chem. Phys. 119, 3405(2003)
 R. LopezRendon, Y. Reyes, and P. Orea, J. Chem. Phys. 125, 084508 (2006)
 L. Verlet, Phys. Rev. 165, 201(1968)
 A. Lofti, J. Vrabec, and J. Fisher, Mol. Phys. 76, 1319 (1992)