Third-order thermodynamic perturbation theory for effective potentials that model complex fluids

We have performed Monte Carlo simulations to obtain the thermodynamic properties of fluids with two kinds of hard-core plus attractive-tail or oscillatory potentials. One of them is the square-well potential with small well width. The other is a model potential with oscillatory and decaying tail. Both model potentials are suitable for modeling the effective potential arising in complex fluids and fluid mixtures with extremely-largesize asymmetry, as is the case of the solvent-induced depletion interactions in colloidal dispersions. For the former potential, the compressibility factor, the excess energy, the constant-volume excess heat capacity, and the chemical potential have been obtained. For the second model potential only the first two of these quantities have been obtained. The simulations cover the whole density range for the fluid phase and several temperatures. These simulation data have been used to test the performance of a third-order thermodynamic perturbation theory TPT recently developed by one of us S. Zhou, Phys. Rev. E 74, 031119 2006 as compared with the well-known second-order TPT based on the macroscopic compressibility approximation due to Barker and Henderson. It is found that the first of these theories provides much better accuracy than the second one for all thermodynamic properties analyzed for the two effective potential models.


I. INTRODUCTION
Complex fluids, which are present in a variety of fields such as soft-matter physics, biophysics, and colloid science, have received much attention in past decades.The basic feature ͓1͔ of the complex fluids is the extremely large asymmetry in size, charge, number of degrees of freedom, and shape between the constituent particles.The multicomponent character and the large asymmetry between the particles of the different species in these mixtures make it cumbersome to carry out the statistical thermodynamic treatment considering explicitly each species in the mixture.Instead, theoretical investigation about the phase behavior of the complex fluids often has to resort to a so-called effective singlecomponent macrofluid approximation.In the singlecomponent approximation, the large colloidal particles interact with each other through an effective potential consisting of a direct colloid-colloid interaction and a solvent-induced potential.After obtaining the effective potential between colloidal particles in one way or another, the complex fluids can be handled by means of statistical mechanics theory for a single-component "atomic fluid." One of the subjects of more remarkable interest at present, in the field of complex fluids, is the study of their phase diagram ͓2-9͔.To theoretically tackle this problem with success, one needs a theory capable of providing reliable results for the pressure and the chemical potential.
Two main features of the effective potential make it necessary to revise the traditional liquid-state theory to adapt it to deal with these effective potentials.One is that the effective potential is very short ranged, to an extent that depends on the size asymmetry ͓10͔ of the particles in the original complex fluid.The other is an oscillatory behavior arising from the packing effect ͓11͔ of the solvent particles.For long-ranged and monotonically decaying potential functions, the traditional thermodynamic perturbation theories ͑TPTs͒ are in general whole satisfactory.Among them, one of the most outstanding is the well-known second-order TPT based on the macroscopic compressibility approximation ͑MCA͒ due to Barker and Henderson ͓12͔.However, truncating the perturbative expansion at second order may not provide enough accuracy for near critical densities and/or low temperatures, especially for short-ranged and/or discontinuous tail potentials.Therefore, for the short-ranged square-well ͑SW͒ potentials considered here, as well as for the potentials with oscillatory tail in which we are also interested here, it is worth analyzing to what extent the traditional second-order MCA TPT is enough accurate.Recently, one of us ͓13͔ has developed a third-order version of the TPT that, even at second order, seems to provide better accuracy than the MCA does.It has been also shown ͓14͔ that for conventional atomic fluid potential the third-order TPT is more accurate than the second-order MCA TPT ͓12͔.
The aim of this paper is to analyze the relative performance of these two theories for predicting the thermodynamic properties of models of colloidal dispersions.To this end, first we compare both theories with simulation data for a number of thermodynamic properties of short-ranged square-well fluids.Many authors have reported simulation data for the equation of state and the excess energy of SW fluids for different values of the reduced well width , but those including data for short ranges-say, ഛ1.2-are quite scarce ͓15-18͔.Still more scarce are the simulation data for the constant-volume excess heat capacity ͓19͔, and there seems to be an absolute lack of simulation data for other thermodynamic properties such as the excess chemical potential for these ranges.Therefore, we have performed Monte Carlo simulations for short-ranged SW fluids with two different values of to determine the above-mentioned and other derived thermodynamic properties.They are presented in the next section and further on used to test the relative accuracies of the two perturbation theories we mentioned previously.We also have performed Monte Carlo simulations for the presently proposed model potential with oscillatory tail with three different sets of parameters which are also presented in the next section and used in this paper to test the reliability of the third-order TPT for obtaining the thermodynamic properties of models of colloidal dispersions.
The organization of the present paper is as follows.In Sec.II the short-ranged model potential and the oscillatory tail model potential are presented.The description of the Monte Carlo simulations is presented in the Sec.III in which we also summarize the third-order TPT and the traditional second-order MCA TPT.The simulation and theoretical results are presented in Sec.IV for the SW potential and in Sec.V for the proposed model potential.A detailed comparison between simulation data and theoretical predictions is also performed in these sections.Finally, we end the paper summarizing the main conclusions arising from this work in Sec.VI.

II. MODEL POTENTIALS
The SW potential model where is the diameter of the particles, − the potential depth, and the range or potential width in units of the diameter , is one of the most widely used in fluid theory.This is because of its ability to describe, at least in a qualitative way, molecular interactions in a wide variety of fluids, together with its great simplicity, which makes easier a theoretical treatment.Therefore, this model potential is often used as a testbed to assess the performance of new theories.
On the other hand, among the systems most widely studied in the last decade in the field of soft-matter physics are colloidal dispersions.These systems are characterized by shortrange interparticle interactions which are often modeled by the square-well potential ͓20-27͔.Another potential model widely used is the hard-core Yukawa potential ͓28-32͔.However, a more detailed analysis of the behavior of the colloidal suspensions requires more sophisticated potential models.
In a colloidal dispersion, due to excluded volume effects, depletion forces arise between the larger particles due to the presence of the smaller ones.In the effective singlecomponent macrofluid approximation, the effect of the solvent particles on the solute particles can be treated as an excess potential of mean force acting on the latter kind of particles ͓33͔.The depletion potential can be obtained from theory, from computer simulation, or from experiment.For a real colloidal dispersion in the same situation, in addition to a depletion potential, another potential due to the direct interaction between the colloidal particles would be present.For a simple model of colloidal dispersions consisting in a binary mixture of hard spheres with diameters 1 ͑the solute͒ and 2 ͑the solvent͒, respectively, with 1 ӷ 2 , in the lowdensity limit of the smaller particles the depletion potential is well described ͓34͔ by the Asakura-Oosawa ͑AO͒ potential model ͓35͔.This is a monotonically decreasing function for separation distances ␦ between the larger spheres in the range 0 ഛ ␦ ഛ 2 and zero for ␦ Ͼ 2 .This is confirmed by experiment.At higher densities of the smaller particles, the depletion potential departs from the AO model and is not longer zero for ␦ Ͼ 2 , but exhibits a damped oscillatory behavior.The asymptotic form of the depletion potential is then ͓34͔ Again, this behavior is confirmed by experiment ͓36͔.
In the present paper, we propose a model potential retaining the main features of these interactions-namely, hardcore plus short-ranged decaying attraction plus long-range damped oscillation-and is given by where where r c is the distance at which the potential is cut and shifted and b, ␣, k s , n, and b are parameters defining the precise form of the potential and are connected with the solvent bath properties by means of very complicated relationships whose explicit mathematical form is beyond the scope of the present paper.The part for b ജ r / ജ 1 reflects the depletion attraction ͓35͔ induced by the finite size of the solvent particles, and the part for b Ͻ r / reflects the packing effect of the solvent particles.By properly choosing the parameters, the potential ͑3͒ and ͑4͒ not only mimics the shape of the potential arising in depletion interactions, but also other kinds of interactions, such as the effective interionic pair potentials in some alkali metals ͓37͔ and alloys ͓38͔.
The shape of the potential ͑3͒ and ͑4͒ is shown in Fig. 1 for three combinations of the potential parameters.

A. Monte Carlo simulations
We have performed Monte Carlo simulations in the canonical ͑NVT͒ ensemble for SW fluids with = 1.05 and reduced temperatures T * ϵ kT / = 0.50 and 0.80 and for = 1.15 and reduced temperatures T * = 0.68 and 0.90 The range of reduced densities considered was * = 0.1-1.0, with step 0.1, where * = 3 and = N / V is the number density.
To obtain the compressibility factor Z = pV / NkT, the excess energy U ex / N, and the constant-volume excess heat capacity C V ex / Nk, the system considered consisted of N = 500 particles, initially placed in a low-density fcc configuration.Next, the particles were allowed to move while increasing their diameters until the desired reduced density * was achieved.Then the system was equilibrated at the desired reduced temperature T * for N e =2ϫ 10 4 cycles, each cycle consisting in an attempted move per particle.The acceptance ratio was fixed at around 50%.The thermodynamic properties were determined from averages performed over the next N c =10 6 cycles, with 100 partial averages from which the statistical uncertainty was estimated as the standard deviation.Such a considerable number of cycles were needed in order to obtain satisfactory accuracy for the constant-volume excess heat capacity, which was determined from the fluctuations of the energy in the canonical ensemble, through the exact relationship C V = ͑͗U 2 ͘ − ͗U͘ 2 ͒ / kT 2 .If only the excess energy and the equation of state were needed, a much lower number of cycles would have been sufficient.Using the simulation data for g͑r͒, the compressibility factor Z was obtained from the virial equation and the excess energy U ex / N was determined from the energy equation.The reduced pressure was obtained as p * ϵ p 3 / = * T * Z.These quantities are listed in Tables I and  II.
To determine the excess Helmholtz free energy F ex / NkT, we have used two procedures.The first of them, labeled ͑1͒ in Tables I and II, is thermodynamic integration along isotherms, according to the expression To this end, the simulation data for the compressibility factor Z were fitted to a suitable polynomial in terms of * .The second procedure, labeled ͑2͒ in the same tables, uses the exact thermodynamic relationship For the excess chemical potential ex we used the simulation data obtained in the form described below.Two different procedures were used to determine the excess chemical potential.The first of them, labeled ͑W͒ in Tables I and II, was the Widom test particle insertion method ͓39͔.We took N = 1372, N e =2ϫ 10 4 , and N c =5ϫ 10 4 .Statistical uncertainty was determined as the standard deviation from 100 partial calculations performed every 500 cycles, each with N t =10 7 trial insertions.
The second procedure, labeled ͑TI͒ in the tables, used the exact thermodynamic relationship ͓40͔ where ␤ =1/ kT, the subscript 0 refers to the hard-sphere fluid with the same reduced density, and u ex = U ex / N. In order to apply the preceding equation, we performed MC NVT simulations to obtain the excess energy U ex / N for different values of ␤, from 0 to 1.6 for = 1.15 and from 0 to 2.0 for = 1.05, along each of the isochors considered.The simulations were carried out with N = 500, N e =2ϫ 10 4 , and N c =5 ϫ 10 4 .The simulation data for ␤U ex / N were fitted to a suitable polynomial in ␤ for each of the densities considered, from which the integration involved in Eq. ͑7͒ was carried out.For Z 0 we used the equation of state derived by Kolafa et al. ͓41͔ by fitting their simulation data.This equation and the simulation data in which it is based are considered the most accurate available at present.The use of the Kolafa equation for performing the thermodynamic integration of the simulation data is justified to avoid any additional source of error ͑although very small͒ in the "exact" data.However, the most simple, but still very accurate, Carnahan-Starling ͑CS͒ equation ͓42͔ is completely satisfactory for theoretical calculations.The chemical potential 0 ex for the hard-sphere fluid was determined from thermodynamic integration of its equation of state, Z 0 , using Eqs.͑5͒ and ͑6͒.
For fluids with the potential form ͑3͒ and ͑4͒, we have performed MC NVT simulations for three parameter sets: ͑a͒ b = 1.15, a = 1.0, k s = 12, n = 3, and b = 0.5; ͑b͒ b = 1.1, a = 0.55, k s = 12, n = 3, and b = 1.0; and ͑c͒ b = 1.1, a = 0.8, k s = 16, n = 4, and b = 0.4.The shape of the potential for these three sets of parameters is shown in Fig. 1.In all simulations the cutoff distance was fixed at r c = 4.The temperatures considered were T * = 0.5, 0.8, and 1.1 for the first two sets and T * = 0.4, 0.7, and 1.0 for the third set.In all cases the density range covered was * = 0.1-0.9 with step 0.1.To obtain the compressibility factor Z and the excess energy U ex / N, from the virial and the energy equations, respectively, we took N = 500, N e =2ϫ 10 4 , and N c =5ϫ 10 4 .Statistical uncertainty was determined as the standard deviation from 100 partial calculations.The reduced pressure p * was calculated from the simulation data for Z in the form indicated before.The excess Helmholtz free energy F ex / NkT was determined by means of thermodynamic integration ͑TI͒ from Eq. ͑5͒.The results are shown in Tables III-V.

B. Thermodynamic perturbation theory
In the three-order TPT ͓13͔, the excess Helmholtz free energy F ex of the system is given by where N is the particle number, = N / V is the number density, and V is the volume occupied by the system.F ex-ref is an excess Helmholtz free energy of a reference hard-sphere fluid with a potential u ref , and u per ͑r͒ is the perturbation part of the whole potential u͑r͒ given by In the present form of the third-order TPT ͓13͔, u ref ͑r͒ is the hard-sphere potential given by g͑r , ␣ , ͒ is the radial distribution function ͑RDF͒ of the bulk fluid with pair potential u͑r ; ␣͒ given by u͑r;␣͒ = u ref ͑r͒ + ␣u per ͑r͒.͑12͒ ͉ ‫ץ‬ ͑n−1͒ g͑r,␣,͒ ‫␣ץ‬ ͑n−1͒ ͉ ␣=0 is the ͑n −1͒th derivative evaluated at ␣ =0 of g͑r , ␣ , ͒ with respect to ␣, and ͉ ‫ץ‬ 0 g͑r,␣,͒ ‫␣ץ‬ 0 ͉ ␣=0 = g͑r ,0,͒ is the RDF of the hard-sphere fluid with density and diameter .
On the other hand, the excess Helmholtz free energy in second-order MCA TPT ͓12͔ is given by The RDF g͑r , ␣ , ͒ was obtained by solving the Ornstein-Zernike ͑OZ͒ integral equation theory ͑IET͒ approximating the bridge function by an accurate expression for the hardsphere fluid developed by Malijevský and Labík ͓43͔.This approximation is based on the fact that we are interested in the region ␣ close to 0, where it is expected that the bridge function will not be very different from that corresponding to the reference hard-sphere fluid.The derivatives ͉ ‫ץ‬ ͑n−1͒ g͑r,␣,͒ ‫␣ץ‬ ͑n−1͒ ͉ ␣=0 for n =1,2,3 were calculated numerically by finite differences.F ex-ref and the compressibility 1 ␤ ͑ ‫ץ‬ ‫ץ‬P ͒ ref of the reference hard-sphere fluid are calculated by means of the CS equation of state ͓42͔.The Verlet-Weis expression for g͑r ,0,͒ ͓44͔ was used in Eq. ͑13͒.The reader can consult Ref. ͓13͔ for further details.
After determining F ex / NkT, the other thermodynamic quantities are obtained by simple differentiation manipula-tion.Thus, the reduced excess chemical potential ␤ ex is given by from which the reduced pressure P * = P 3 / immediately follows: where ␤ = ␤ ex +ln * and ␤f = F ex / NkT +ln * −1.
The reduced excess energy U ex / N is obtained in the form which in turn allows us to obtain the reduced constantvolume excess heat capacity C V ex / Nk:

IV. RESULTS FOR SHORT-RANGED SQUARE-WELL FLUIDS
The results from the second-order MCA TPT and the third-order TPT for the thermodynamic properties of short-  I and II.Two facts emerge clearly from these figures: ͑a͒ the third-order TPT provides excellent agreement with simulation data, even at low temperatures, for all the thermodynamic properties considered, but in some cases for the constant volume excess heat capacity C V ex / Nk, and ͑b͒ the third-order TPT constitutes a strong improvement over the second-order TPT-MCA The fact that the second-order MCA TPT provides poorer agreement with the simulation data for the excess energy as the potential width becomes shorter, as seen in Fig. 6, is due to the fact that the relative importance of the higherorder perturbative terms increases as the potential width is reduced.Moreover, the MCA strongly underestimates the magnitude of the second-order perturbative term of the Helmholtz free energy, which is directly related to the corresponding term of the excess energy, for short-ranged SW potentials ͓45͔.
Another drawback of the second-order MCA TPT, also seen in Fig. 6, is that the predicted temperature dependence of the excess energy is too small.This is reflected in the poor prediction of the constant-volume excess heat capacity shown in Fig. 3.
In contrast, the fact that the third-order TPT provides very good agreement with the simulation data for the excess energy, even for short-ranged SW potentials, not only is due to include the third-order perturbative term, but also implies that the prediction of the second-order term by this theory is good.Now, the temperature dependence of U ex / N is more accurately reproduced.As a consequence, this theory predicts values of C V ex / Nk in much closer agreement with simulation data than those from the second-order MCA TPT, as seen in Fig. 3.However, to achieve complete accuracy with simulation data for C V ex / Nk, it would be necessary the incorporation of higher-order perturbative terms.
One might argue that we are comparing a second-order TPT with a third-order one and that, if the Barker-Henderson perturbation theory were extended to third or higher order, the differences between the two theories would be reduced.Such an extension is possible, because the procedure developed by Barker and Henderson to obtain the second-order perturbative contribution to the free energy in the macroscopic compressibility approximation was generalized ͓46,47͔ to higher-order terms and the resulting infinite perturbative series was summed up to obtain a closed expression for the excess Helmholtz free energy.However, the difference between the results obtained in this way and those resulting from the perturbative series truncated at the level of the second-order term are nearly negligible.On the other hand, instead of the macroscopic compressibility approximation we could have used the local compressibility approximation ͑LCA͒ ͓12͔ for the second-order term in the Barker-Henderson TPT.However, at least for SW fluids, the difference between the two approximations is small ͓19,45͔.In Figs.2-6 it seems that the relative accuracy of the third-order TPT increases as the well width decreases, and this becomes more apparent in Fig. 3.However, it is to be noted that, from recent simulations ͓48͔, the reduced critical temperatures for = 1.05 and for = 1.15 are T c * = 0.366 and T c * Ϸ 0.57, respectively.Therefore, the lowest reduced temperature considered for = 1.05 is more supercritical than that for = 1.15, and the same is true for the highest temperatures considered in the two cases.This explains what at first sight may appear quite surprising.It is worth mentioning here that the liquid-vapor coexistence becomes mestastable with respect to solid-gas coexistence for of the order of 1.25 or lower ͓48͔.
The fact that the third-order TPT accurately predicts the pressure and the chemical potential for short-ranged SW fluids suggests that it might be useful for predicting the liquidvapor coexistence.For the reasons just pointed out, it is difficult to obtain reliable simulation data for the coexistence curve below = 1.25.Therefore, we have considered two intermediate ranges-namely, = 1.25 and = 1.375.
In the third-order TPT, the coexistence curve is obtained by equating chemical potentials and pressures in both phases at a given temperature.These quantities are easily determined from the total free energy.The solution of the equi-  librium conditions corresponds to a double-tangent construction on the curve of the free energy per unit volume versus the density.To this end, the free energy curve is obtained numerically at different temperatures.As the numerical implementation of the third-order TPT involves the calculation of g͑r , ␣ , ͒ corresponding to the pair potential u͑r ; ␣͒ given by Eq. ͑12͒, then one might think that a problem would be the possible breakdown of the code when the considered temperature is below the critical temperature.Such a problem actually does not arises because, to obtain numerically the derivatives ‫ץ͉‬ ͑n−1͒ g͑r , ␣ , ͒ / ‫␣ץ‬ ͑n−1͒ ͉ ␣=0 with n =1,2,3, involved in Eq. ͑9͒, one only needs to calculate the values of g͑r , ␣ , ͒ for ␣ =0, Ϯ⌬ , Ϯ 2⌬, where ⌬ is a small quantity-for example, 0.005.For such small values of ␣, u͑r ; ␣͒ is very close to the hard-sphere potential or, equivalently, very close to the infinite-temperature limit.Therefore, the calculations of the perturbative terms are performed at a temperature that is always much higher than the critical temperature of the true potential u͑r͒.This is exactly why we can use the hard-sphere bridge function for the potential u͑r ; ␣͒ when ␣ is very small, as we mentioned before.From another viewpoint, solving the OZ integral equation for u͑r ; ␣͒ at temperature T is actually equivalent to solving the OZ integral equation for u͑r͒ at temperature T / ␣.When ␣ is sufficiently small, the equivalent temperature T / ␣ for the potential u͑r͒ is sufficiently high and always higher than the critical temperature of the potential u͑r͒.
The results for the liquid-vapor coexistence obtained from the third-order TPT and those from the Barker-Henderson second-order TPT in the MCA are compared in Fig. 7 with the simulation data ͓49-51͔.We have included for comparison the results from the self-consistent OZ approximation ͑SCOZA͒, recently extended to SW fluids ͓52͔, although it belongs to a different class of theories.The figure shows that the third-order TPT compares favorably with the other two theories.

V. RESULTS FOR THE OSCILLATORY POTENTIAL MODEL
In Figs.8-10, the results from the two perturbation theories for the thermodynamic properties of fluids with interparticle interactions obeying to the oscillatory potential ͑3͒ and ͑4͒ are compared with the simulation data from Tables III-V.The same conclusions can be drawn about the relative performance of both perturbation theories.Now, the superior quality of the third-order TPT over the second-order MCA TPT is particularly noteworthy for the excess energy, as  I and II.Squares and solid circles, nearly indistinguishable from each other at the scale of the figure, correspond to procedures ͑1͒ and ͑2͒, respectively.The curves have the same meaning as in Fig. 2.
shown in Fig. 10.In this figure stands out a quality also present, but less clearly, in Fig. 6 for the excess energy of the SW fluid-namely, the fact that the third-order TPT clearly differentiates the energy of isotherms with temperatures close to each other.In contrast, the second-order MCA TPT does not.This is because, as pointed before, in the latter theory the second-order term in the energy, from which depends the variation with temperature of the excess energy, U ex / N, is usually very small.
In the previous section we have pointed out that the thirdorder TPT seems accurate enough so as to provide reliable estimates of the liquid-vapor coexistence.Thus, we have used it to estimate the critical temperatures for the three sets of parameters considered for the oscillatory potential ͑3͒ and ͑4͒.They are T c * = 0.3365, 0.3692, and 0.2989 for sets I, II, and III, respectively.Therefore, the lowest temperature considered in the simulations and the theoretical calculations for each of the three sets is higher than, but close to, the critical temperature.
However, it is to be noted that, among the three sets of parameters considered, the third-order TPT, in combination with a first-order TPT ͓53͔ for the solid phase, only predicts a stable liquid phase for set II.For the two other sets of parameters, the liquid-vapor coexistence is metastable with respect to the solid-gas coexistence.It has been shown ͓54-56͔ that in models of colloidal fluids with short-range attractive and weak long-range repulsive competing interactions, a microphase separation can take place, with large stabilized clusters, at temperatures below the critical temperature, preventing gas-liquid separation.This situation is not present in the states we have studied for the oscillatory potential, because all the temperatures considered are supercritical.In any case, with the sets of parameters considered it is not expected that a stable microphase separation will occur even at subcritical temperatures, because the repulsive interactions, apart from the hard core, are not enough weak nor long ranged.Moreover, the long-range oscillatory tail of the potential hinders the formation of such microphases.

VI. CONCLUSIONS
In the preceding sections we have provided compelling evidence of the great accuracy of the third-order thermody- namic perturbation theory ͓13͔, summarized in Sec.III, for most of the thermodynamic properties considered.In fact, this theory is much more accurate than the second-order Barker-Henderson perturbation theory in the MCA.The reasons why the performance of the third-order TPT is superior to that of the Barker-Henderson ͑BH͒ secondorder TPT are quite obvious: on the one hand, the BH perturbation theory predicts values for the second-order perturbative term in the energy that are too small in magnitude, especially for short-ranged potentials, as shown in Ref. ͓45͔, whereas the second-order term in the third-order TPT is much more accurate ͓13͔.This is due to the fact that, in the latter, the second-order term is associated with the first-order derivative of the radial distribution function g͑r , ␣ , ͒ evaluated at ␣ = 0, as required rigorously by the TPT, whereas in the BH second-order TPT the same term is evaluated from the RDF g͑r ,0,͒ of the reference fluid, which is a more crude approximation.On the other hand, the third-order TPT includes an additional term, and it is known that higher-order terms become increasingly important with decreasing potential range, temperature, and density.In contrast, adding higher-order terms to the BH TPT in the MCA does not result in any significant improvement, as they are nearly negligible.
From the preceding considerations, we can conclude that the performance of the third-order perturbation theory developed by one of us ͓13,14͔ in predicting the thermodynamic properties of fluids with short-ranged and oscillatory potentials is very satisfactory.This makes it particularly suitable to deal with complex fluids and colloidal suspensions.A remarkable fact is that the theory is particularly accurate for the pressure and the excess chemical potential, two quantities that play an essential role in studying phase transitions.
For the constant-volume excess heat capacity C V ex / Nk instead, the third-order TPT is still insufficient to provide enough accuracy at very low temperatures and higher-order terms would be needed.
FIG.1.Shape of the potential model defined by Eqs.͑3͒ and ͑4͒.
FIG.2.Excess chemical potential for the square-well fluids considered.Points: simulation data from TablesI and II.Squares and solid circles, nearly indistinguishable from each other at the scale of the figure, correspond to TI and W procedures, respectively.Dashed curves: Barker-Henderson second-order MCA TPT.Solid curves: Zhou third-order TPT.
FIG. 3. Constant-volume excess heat capacity for square-well fluids.Points: simulation data from TablesI and II.The curves have the same meaning as in Fig.2.
FIG. 5. Reduced pressure for the square-well fluids considered.Points: data from TablesI and II.The curves have the same meaning as in Fig.2.
FIG. 7. Liquid-vapor coexistence for square-well fluids.Points are the simulation data from Refs.͓49-51͔, respectively, for squares, triangles, and diamonds.Solid lines are the results from the third-order TPT, dashed lines are the results from the SCOZA ͓52͔, and the short-dashed lines are the results from the second-order MCA TPT.
FIG. 8. Excess Helmholtz free energy for fluids with the oscillatory potential ͑3͒ and ͑4͒.Points: data from Tables III-V.Dashed curves: Barker-Henderson second-order perturbation theory in the MCA.Solid curves: third-order TPT of this work.
FIG. 10.As in Fig.8for the reduced excess energy.

TABLE I .
Simulation results for the SW fluid with = 1.05 ͑see explanations in the text͒.The numbers between parentheses are the statistical uncertainties in the last decimal place.

TABLE II .
Simulation results for the SW fluid with = 1.15 ͑see explanations in the text͒.The numbers between parentheses are the statistical uncertainties in the last decimal place.

TABLE III .
Simulation results for a fluid with the model potential ͑3͒ and ͑4͒ with the first parameter set.

TABLE IV .
Simulation results for a fluid with the model potential ͑3͒ and ͑4͒ with the second parameter set.

TABLE V .
Simulation results for a fluid with the model potential ͑3͒ and ͑4͒ with the third parameter set.