Thermodynamic perturbation theory and equation of state developments

An alternative way of utilizing the thermodynamic perturbation theory of Wertheim for the development of equations of state for associating fluid models is presented and detailed for water. The approach makes use of general features of the parameter of non-saturation to avoid the necessary solution of an algebraic equation and unbinds the results from a tight dependence on details of the simple reference fluid model used in the perturbation theory.


Introduction
Because of its indispensable role, both for life and in industry, water has been in the focus of research for centuries. It exhibits a number of anomalies which are, more or less, well understood but for applications it is also necessary to express its properties by means of closed analytic expressions. Such expressions may be obtained by a parametrization of available experimental data as, e.g., the IAPWS-90 equation [1] but the applicability of such equations is limited only to the range of the data used in the parametrization. The other possibility is the use of statistical mechanics which provides a link between the atomistic view of matter and the observed macroscopic behavior. Regardless of details, equations of state (EoS), both strictly theoretical or semiempirical, possess the form of a perturbed equation, where is the compressibility factor, is the pressure, is the inverse temperature, = 1/ B , where B is the Boltzmann's constant and is the temperature, is the number density, = / , ref is the compressibility factor of a certain reference system, and Δ is a correction term. The above form of the EoS results from a decomposition of the total intermolecular interaction model into a reference part, ref , and a perturbation part, Δ , where (1,2) stands for a set of generalized coordinates. Within the spirit of perturbation theories of liquids, the reference system has to reproduce the structure of the original liquid, and its properties should be available in an analytic form. As regards the associating liquids, simple models exhibiting hydrogen bonding (H-bonding), usually referred to as primitive models (PM), are employed to meet this requirement. Such simple models emerged in the end of the 1980's [2][3][4][5] and, nearly simultaneously, theories of these models were developed, a theory of Dahl and Andersen [6] and the thermodynamic perturbation theory (TPT) of Wertheim [7,8]. These two ingredients together, models and theory, opened up the way to the development of molecular-based equations following either (i) truly molecular approach which starts with a realistic force field and then, by means of a sequence of well defined approximations, a simple (primitive) model amenable to a theoretical treatment is constructed or (ii) a van der Waals (vdW)-like approach which does not refer to any specific interaction model but only makes use of the previously obtained knowledge and constructs a reference intuitively/empirically. Both approaches use the TPT to obtain the properties of the reference model. Despite the difference, both approaches may end up, formally, with the same EoS.
Despite all the effort invested into the development of an accurate and reliable EoS for water, no satisfactory results have been achieved so far. A theoretical route seems to have culminated with an accurate description of a short-range reference, i.e., primitive model descending from the TIP4P force field [9,10]; the equation reproduces the main anomalies exhibited by water but the attempts to use it as a reference for the development of a truly molecular-based EoS for water were not fully satisfactory [11][12][13]. The vdW-like approach is represented by statistical association fluid theory (SAFT) equations. As Llovel and Vega stated in their recent review [14], there were 47 SAFT equations for water available, all of them based on a sort of a primitive model. Nonetheless, none of them meets the demands for widespread and reliable applications which only witnesses to their empirical nature.
Application of the TPT is not straightforward. The property which links the model with the theory lies in certain integrals , which should to be evaluated for the model used. They involve correlation functions of the reference model and the Mayer function of the H-bonding interaction (for details see the next section). To evaluate these integrals analytically, certain approximations should be used. Once the integrals are known, then it is necessary to evaluate the key function of the theory -a non-saturation parameter . This evaluation lies in solving an algebraic equation whose order depends on the order of the used TPT.
Attempting to improve the performance of the EoS developed by the intuitive vdW-like approach, there are no other ways than (i) to modify/change the reference model, or (ii) to improve the performance of the TPT, and/or (iii) to play games with the correction terms in the same way as in the development of hundreds of cubic equations. As for the two former ways, it is evident that all the computations should be always performed again from the scratch and there is a question whether this effort is worth doing because the imperfectness of the reference model should be also reflected in the final results anyway. The theoretical approach offers another possibility. Results based on the actual underlying physics should be qualitatively the same regardless of their parameter values. It means, the resulting functions should exhibit the same course/functional dependence. Consequently, focussing on general features of the derived properties as, e.g., the parameter of non-saturation, we may skip over the details of the reference model and thus release the direct dependence of the EoS on the reference model. Furthermore, this approach will also introduce additional adjustable parameters, in addition to those defining the reference model, which should increase the chance of the final equation to perform better.
The purpose of this paper is to examine the properties of the derived functions and to find their general course, if possible at all. At a feasibility level, the suggested route towards an EoS is examined by considering a model descending from the TIP4P force field, the model with probably the most appropriate geometry of the water molecule. In the next section we provide the basic definitions and then in the subsequent section the results. The final section contains a summary where a future potential development is also briefly outlined.

Basic methodological details
The thermodynamic perturbation theory is not a general theory since it requires a corresponding form of intermolecular interactions, where 0 is a reference system interaction potential, usually the hard sphere model, and the second term is a site-site interaction giving rise to H-bonding represented by a short-range attractive site-site interactions HB between sites and from set Γ. For water, typically |Γ| = 4 if all sites (i.e., both those bearing a 33501-2 charge and uncharged ones) are counted but this value is not strictly defined and may also be both 2 and 3 in other models [15]. The basic quantity of the theory are fundamental integrals involving the Mayer functions of the H-bonding interaction between sites and , and correlation functions, , of the reference system. Denoting by (+) the hydrogen-like sites and by (−) the oxygen-like sites, then the 1st order theory requires only integral +− defined as where is the pair correlation function. The actual form of in the second order of the TPT depends on additional approximations. In the 2nd order TPT as developed by Slovak et al. [16], integral contains the triplet correlation function, : The basic quantity of the theory is a parameter of non-saturation, . It is obtained by solving an algebraic equation with coefficients given by integrals and whose order depends on the level of approximations. For single bonding models it is a quadratic equation; for double bonding models within the second order it is, in general, an equation of the sixth order but it can be reduced, after introducing an additional approximation, to a cubic equation; for details see, e.g., [16] or a pedagogical introduction to the TPT by Zmpitas and Gross [17]. An additional quantity appearing in the final expressions, but which is a function of , is a parameter 0 (in notation of Slovak et al. [16]). In a general case of molecules containing hydrogen-like (+)-sites satisfying the condition of single bonding, and one complementary (−)-oxygen-like site which may form up to two bonds, the Helmholtz free energy, , of the reference fluid assumes, within the TPT, the form Δ PM / = Δ 0 / + (1 − ) + ln 0 . (2.5) Integrals +− and +−+ were evaluated by simulations and consequently parametrized so that all information necessary for the analytic evaluation of parameters and 0 , and hence the Helmholtz free energy and all other thermodynamic properties of the primitive model, is available.

Results and discussion
To examine the above outlined idea on the application of the TPT, we consider the primitive model descending from the TIP4P force field as developed by Vlcek and Nezbeda [18] along with analytic parametrizations of integrals . We have evaluated parameters and 0 along 9 isotherms within the range (300-700) K for packing fractions up to a typical water value = 0.35.
Parameter obtained from numerical computations along with its parametrization is shown for three selected temperatures in figure 1 in dependence on density for three temperatures.
As it is seen, the logarithmic fit, provides an excellent fit to the numerical data and the same also applies to the exponential parametrization of parameter 0 , shown in figure 2. The value of 0 for = 300 is around zero and would not be thus distinguishable from a zero straight line within the scale of the graph. In fact, the exponential dependence on should not be 33501-3  so surprising because the underlying integrals themselves can be factorized into -and -dependent terms.
The above two functions, and 0 , contain two temperature dependent parameters which determine their curvature. The answer to the question, whether these parameters also exhibit a simple dependence on is provided in figures 3 and 4.
The result is appealing: all these coefficients exhibit a simple quadratic dependence on . In the end, we thus have at most 4×3 constants available for fitting the experimental data.

Conclusions
In common applications of the SAFT methodology, the reference (primitive) model is defined by three parameters, the size of the hard sphere, and the range and strength of the H-bonding attraction. To gain a flexibility and improve the performance of the results, when fitting the experimental data, a temperature dependence of the HS diameter is introduced and the restrictions imposed by the TPT are lifted. These adjustments are rather arbitrary without physical considerations. As an attempt to incorporate/offer at least some physical insight into this amendment, we suggest that a tight link to the reference simple model should be broken by focussing on a general behavior of the derived properties, the parameter of non-saturation, and 0 , instead of fitting directly the parameters of the reference model. To summarize, a rather complicated computation of parameters and 0 for the given reference model may be replaced by general functions of the known temperature and density functional forms ensured by the physical nature.