Thermische Verfahrenstechnik
Analysis and Extension of Wertheim's
Thermodynamic Perturbation
Theory
Wertheim’s Thermodynamic Perturbation Theory
Von der Fakult¨at Energie-, Verfahrens- und Biotechnik der Universit¨at Stuttgart zur Erlangung der W¨urde
eines Doktors der Ingenieurwissenschaft (Dr.-Ing.) genehmigte Abhandlung
Vorgelegt von Wasilios Zmpitas
aus Leonberg
Hauptberichter: Prof. Dr.-Ing Joachim Groß
Mitberichter: Prof. Dr. rer. nat. Amparo Galindo Tag der m¨undlichen Pr¨ufung: 2019/06/19
Institut f¨ur Technische Thermodynamik und Thermische Verfahrenstechnik der Universit¨at Stuttgart
Diese Studie bietet eine nachvollziehbare und detaillierte Analyse der zugrunde liegenden Konzepte der von Wertheim in den Jahren 1984-1986 entwickelten Thermodynamic Pertur-bation Theory (TPT). Sie behandelt die fundamentalen Schritte, die ihren Ursprung in der Statistischen Thermodynamik haben, und f¨uhrt den Leser durch die anspruchsvolle Anwen-dung der Graphentheorie, aus der die single-chain approximation und die St¨orungstheorien erster (TPT1) und zweiter Ordnung (TPT2) hervorgehen. Diese Arbeit liefert einen zug¨ angli-chen Einstieg f¨ur Interessierte, die die grundlegenden Konzepte der TPT studieren m¨ochten. Weitergehend werden zwei Erweiterungen innerhalb Wertheims TPT entwickelt. Dazu geh¨ort zum einen die Entwicklung einer St¨orungstheorie dritter Ordnung (TPT3), und zum an-deren eine Theorie, die die Dimer-Dimer Graphensummen ber¨ucksichtigt und damit ¨uber die single-chain approximation hinaus geht. Beide Erweiterungen werden graphentheoretisch hergeleitet und als analytische Ausdr¨ucke in Abh¨angigkeit von 3- und 4-Teilchenverteilungs-funktionen formuliert. Die genannten Verteilungs4-Teilchenverteilungs-funktionen werden mit Hilfe von Monte Carlo Simulationen bestimmt. Die entwickelten Theorien werden im Limit der vollst¨andigen Assoziation auf Hartkugelfluide mit einem und zwei Wechselwirkungszentren angewendet und mit Simulationsdaten von reinen Hartkettenfluiden mit Kettenl¨angen 4, 8, 16 und 64 verglichen. Ein weiterer wesentlicher Punkt dieser Studie ist die Bewertung des Einflusses weiterer St¨orungstheorieordnungen hinsichtlich der Vorhersage des zweiten Virialkoeffizien-ten f¨ur Hartkettenfluide im Niedrig-Dichte-Limit. Es ist bekannt, dass TPT1, die in fast jeder heutigen Statistical Associated-Fluid Theory (SAFT) Zustandsgleichung verwendet wird, in der Vorhersage des Niedrig-Dichte-Limits im Vergleich zu Hartkettenmodellsim-ulationen deutliche Abweichungen aufzeigt.
This study provides a comprehensible and detailed analysis of the underlying concepts of Thermodynamic Perturbation Theory (TPT) developed by Michael S. Wertheim in the years 1984-1986. It covers the fundamental steps which are rooted in Statistical Thermodynamics, and leads the reader through the difficult applications of graph theory, which give rise to the single-chain approximation and the first- and second-order theories TPT1 and TPT2. This work provides an available access for people who want to read up and work on the basic concepts of TPT. It further shows futher developement in the TPT framework by extend-ing Wertheim’s theories beyond the second-order (TPT2) by both, derivextend-ing a third-order expansion in the single-chain approximation, and going beyond the single-chain approxima-tion by including the infinite dimer-dimer graph sum in the graph-theoretical developement. Both extensions are derived based on graph theory and formulated as closed expressions with respect to reference 3- and 4-particle distribution functions, which are evaluated by Monte Carlo simulations. The resulting theories are applied on a hard-sphere reference fluid with one and two attractice sites in the limit of complete association, and compared with simulation data for pure hard-chain systems with chain lengths 4, 8, 16 and 64. Assessing the influence of higher orders in terms of the prediction of the second virial coefficient of hard-chains in the low-density limit is also a vital point of this work. It is well known that TPT1, which is used in nearly every today’s Statistical Associated-Fluid Theory (SAFT) equation of state, disagrees significantly in this limit when compared to hard-chain fluid simulations.
Introduction 1
Fluid theories for attractive interactions of spherically symmetric fluids . . . 2
Fluid theories for highly directional interactions . . . 3
Fluid theories for non-spherical fluids . . . 5
Structure of this thesis . . . 6
1 An analysis and detailed pedagogical review of Wertheim’s Thermody-namic Perturbation Theory 7 1.1 Introduction . . . 8
1.2 Classical Thermodynamics . . . 8
1.3 Basics of graph theory . . . 12
1.4 Molecular Model . . . 18
1.5 Wertheim’s TPT . . . 21
1.6 Assumptions and approximations . . . 33
1.6.1 One attractive site A . . . 34
1.6.2 One attractive site: a brief discussion . . . 37
1.6.3 Two attractive sites A and B . . . 39
1.6.4 Two attractive attractive sites: a brief discussion . . . 48
1.6.5 Multiple attractive sites: SAFT formalism of W. G. Chapman . . . . 49
1.6.6 Multiple attractive sites - SAFT formalism: A brief discussion . . . . 53
1.7 Conclusion . . . 55
Appendices 57 1.A s-Particle Distribution Functions g(s) and G(s) . . . . 57
1.A.1 2-Particle Distribution Function g(2) and G(2) . . . . 57
1.1.2 3-Particle Distribution Function g(3) and G(3) . . . 58
1.1.3 4-Particle Distribution Function . . . 59
2.2 Molecular model . . . 64
2.3 Homonuclear, flexible chains . . . 65
2.3.1 n-th-order TPT . . . 66
2.4 Third-order TPT . . . 68
2.5 Conclusion . . . 73
Appendices 77 2.A Molecular simulation . . . 77
2.A.1 Distribution functions at contact distance . . . 77
2.A.2 Compressibility factor . . . 80
3 Extension of Wertheim’s thermodynamic perturbation theory to include higher order graph integrals 85 3.1 Introduction . . . 86
3.2 Model . . . 86
3.3 Thermodynamic Perturbation Theory . . . 87
3.3.1 One Interaction Site, Γ = {A} . . . 88
3.3.2 Two Interaction Sites, Γ = {A, B} . . . 91
3.3.3 Limit of complete association . . . 95
3.4 Results . . . 95
3.4.1 Simple relations for correlation integrals λ2 and λ11 . . . 95
3.4.2 Hard-spheres with two interaction sites . . . 98
3.5 Conclusion . . . 98
3.6 Acknowledgment . . . 99
Appendices 101 3.A Correlation functions G(3)hs(123) and G(4)hs(1234) . . . 101
3.B Graphical derivation of G(4)hs(1234) . . . 105
lar model. A molecular model is a simplified description of intermolecular and intramolecular energies, thus defining the forces within and between molecules. Physical models are partic-ularly advantageous, when considering properties at conditions not covered by experiments. In contrast, empirical models (that often give expressions for the excess molar Gibbs energy gE) adjusted to a limited range of experimental data usually extrapolate rather poorly to
conditions outside of the training range. The advantages of physically-based models are pro-nounced when dealing with mixtures at fluid conditions. Predictions are often at least qual-itatively correct without requiring mixture-parameters to be adjusted to experimental data of the targeted mixture. The explicit description of intermolecular interactions is especially relevant when modeling mixtures containing substances with strongly directional interac-tions. Highly directional interactions between molecules, like hydrogen-bonding, need to be described by a suitable approach in order to preserve predictive power in a thermodynamic model. At the same time, it is rather difficult to derive a meaningful model from Statistical Mechanics for directional interactions, due to the complex intermolecular configurations that appear in such systems.
Fluid theories for attractive interactions of spherically
symmetric fluids
The pioneering work of J. D.van der Waals[1] provided a first perturbation theoretical ap-proach to describe substances in fluid conditions, including gaseous and liquid states, and J. Williard Gibbs[2] gave a clear Statistical Mechanical treatment of various facilitating the development of fluid theories (integral equation theories or perturbation theories) in their current form.
A perturbation theory in statistical thermodynamics is usually formulated in terms of Helmholtz energy A, and relates the (pair) potentials of a reference fluid and a target fluid by introducing a coupling parameter λ. Ranging from 0 to 1, the coupling parameter ensures a smooth transition from the reference fluid (with pair potential u0(r) at λ = 0) to the target
fluid (with uλ(r) ≡ u0(r) + λ · up(r) at λ = 1), with
A = A0+ Z 1 0 ∂A ∂λ dλ (1)
where for simplicity a spherically symmetric fluid with only pairwise additive interactions is assumed and, furthermore, defining the coupling parameter λ to act linearly on the potential. From a simulation perspective a smooth transition is expedient for effectively sampling the states in between the limiting cases λ = 0 and 1. Depending on the coupling parameter’s current value, the perturbation formulation describes a fluid with properties which are in between the properties of the reference and target fluid. For deriving an analytical thermo-dynamic model, usually a Tayler expansion for (∂A/∂λ) around the reference state λ = 0 is 2
reciprocal Temperature (kBT )−1), because the derivation gives also rise to a prefactor (ρβ)n
for each order term n. Other perturbation concepts, such as the f -expansion[3], relate the Mayer functions[4] by a coupling parameter. The target fluid’s potential is then not written as a linear function of the coupling parameter λ, but the general form of this simple pertur-bation scheme remains unchanged. Truncating the Taylor expansion results in a model that approximates the target fluid’s behaviour based on reference fluid properties. The accuracy of the model predictions improves as the order of the Tayler expansion increases. If the prop-erties of both, reference and target fluid differ strongly, it is necessary to include higher order terms. In general, it is possible to use an arbitrary reference fluid. However, to ensure a fast converging Taylor expansion, the choice of a well-suited reference fluid that shows similar characteristics as the target-fluid is important. For this purpose, the Helmholtz energy and the structure, as captured in two- and higher-body correlation functions, of the reference state must be well known. An often-used and well-known reference fluid for perturbational applications is the hard-sphere model[5–7]. It is simple in two aspects: first, the properties of the hard-sphere fluid can be written as a function of density only (instead of temperature and density) and second, the hard-sphere fluid shows only a liquid-solid transition, but no liquid-vapor transition.
Theories for simple, spherical fluids, such as the hard-sphere fluid or the Lennard-Jones fluid, were developed in the 1950s to 1970s, catalyzed by the advent of molecular simulations. The first description of attractive fluids based on a perturbation expansion about a hard-sphere fluid was due to Robert Zwanzig[8] in 1954. These developments are well documented in several reviews[3, 9, 10].
Fluid theories for highly directional interactions
The first theory for dipolar fluids, i.e. fluids with directional interactions, was developed by John A. Barker in 1951[11], although, at the time it was not possible to adequately assess the theoretical approach, because properties of a reference fluid were not available with sufficient accuracy and in particular expressions for pair correlation functions were not yet known. About twenty years later, in 1972 the field has matured to proceed on the development of theories for dipolar interactions, with Stell et al. first developing a perturbation expansion to third order in Helmholtz energy with a Lennard-Jones fluid as a reference fluid[12]. A strong effect of the dipole moment on thermodynamic properties was found for liquid states. Even a third order theory was seen to be insufficient, which means that the arrangement of dipolar molecules is substantially different from the structure of the reference fluid (i.e. the fluid without the dipolar moment), which motivated introducing a Pad´e approximant[13] to the Helmholtz energy terms of second and third order.
Hydrogen bonds are highly directional interactions that are, compared to dipolar inter-3
bation approach for hydrogen-bonding interactions in a Helmholtz energy expansion further deteriorates. It is common to refer to directional interactions, mimicking for example hy-drogen bonds, as association. The first non-lattice approaches to associating fluids were developed by Anderson using graph theoretical expansions in density[14–16]. A strong ap-proach was also developed by Pratt and Chandler with a graph theoretical development written in terms of fugacities[17, 18]. Høye and Olaussen[19] proposed a graph theory in terms of density and in terms of monomer density, a theory that is related in several aspects to the later work of Wertheim.
In the years 1984-1986, Michael S. Wertheim developed a powerful framework for model-ing associatmodel-ing fluids [20–23], known as Thermodynamic Perturbation Theory (TPT). For-mulating a simple perturbation theory, eq. (1) in terms of a reference fluid with simple, spherically symmetric pair interactions (such as a hard-sphere fluid) for describing highly directional interactions does not lead to accurate predictions using a nth-order Taylor
expan-sion (with reasonable n). The reason is obvious: the difference between the intermolecular structure of the reference and the target fluid is very pronounced. In this context, Wertheim’s TPT allows to describe a target fluid with highly directional interactions with respect to a simple reference fluid, with a spherically symmetric pair potential, like the mentioned hard-sphere model for example. To achieve this task, Wertheim turned away from the general single-density conception and introduced a multiple-density perturbational formalism, which is developed in a graph-theoretical framework. The multiple-density formalism is general, but it is especially expedient for directional short-ranged potentials, like association or chem-ical bonding. Steric hindrances, which influence short-ranged interactions, can be directly introduced in Wertheim’s theory in an early stage of graph-theoretical development. This allows for applying beneficial bonding constraints that introduce major simplifications and graph cancellation, by prohibiting certain particle configurations. Such a cancelation scheme was first considered by Andersen [14, 15] for particles with acentrical interaction sites.
Graph theory is a very general term for describing a graphical approach with graphs. It has a wide range of applications and is usually applied in mathematics and computer science. For this work’s purpose it provides a framework for relating complex mathematical equa-tions into an accessible, graphical level, without loosing information. Further, mathematical operations are possible for modifying graphs, i.e. for modifying rather involved multidimen-sional integrals. In this work only graph-theoretical tools are included that are necessary to reproduce the derivation of Wertheim’s theory as well as extensions thereof, as introduced in this work. For more background we refer to excellent reviews on graph theory[17, 24, 25]. Wertheim’s original work is captured in his first four publications[20–23]. In his first two studies[20, 21] he derives a TPT for particles with one sole attractive interaction site. The resulting theory is stated in terms of a two-density formalism and recast in a later study[26] into an equation of state for associating hard-spheres with a single site. By restricting each attractive site such that it takes part in only one interaction, his study shows that the equa-4
the theory shows excellent agreement and reveals its great potential in treating highly direc-tional interactions. In his third and fourth publications[22, 23], Wertheim extends his TPT to fluids with multiple sites, which gives rise to an equation that is capable of describing more complex molecules (depending on the applied bonding constraints and the number of interaction sites). For particles with two sites, Wertheim formulates a first-order (TPT1) and a second-order (TPT2) equation of state for linear, chain-like configurations[27]. Linear structures are ensured by, again, restricting each attractive site such that it takes part in only one interaction. He thereby introduces the single-chain approximation into his graph theoretical development which considers only graph sums that consist of one chain-associate with chain length ν surrounded by the reference hard-sphere fluid. Graph sums that contain more than one chain are thus neglected. This approximation allows him to formulate closed expressions of the graph sums in terms of particle distribution functions. Depending on the highest order n of considered (n + 1)-body effects (which is given by the longest chain length ν of a chain-associate in a graph sum i.e), the resulting equation of state is denoted as TPT of n-th order, or simplier TPTn.
Fluid theories for non-spherical fluids
Theories for non-spherical fluids are important for modeling real, molecular substances. Significant work has been done for purely repulsive (i.e. hard body) fluids. The Scaled Particle Theory [28] first proposed in 1959 and later refined for non-spherical bodies[29, 30] and for mixtures [6, 31, 32] is a successful model for hard body fluids. The Scaled Particle theory, however, has limitations for strongly non-spherical, (flexible) chain-like fluids.
The TPT proposed by Wertheim for associating fluids can also regard spherical entities in the limit of irreversible association, thus leading to a model for non-spherical fluids composed of chains of spherical segments. Driving TPT1 and TPT2 into the limit of irreversible association leads to equations of state which describe a mixture of chain fluids with an average number of spherical segments ν. The number of spherical segments of a chain is here referred to as chain ‘length’. We note that Wertheim’s equations of state (TPT1, TPT2) are not capable of describing systems with mixtures of chain fluids explicitly; rather the different chain lengths of a mixture are captured in terms of an average chain length. The underlying molecular model of the original TPT does not introduce any bonding constraints with respect to the chain length. It simply includes two acentric interaction sites A and B within each particle with pair potential uAA = uBB = 0 and uAB < 0, which allows
attractive bonding only between two sites of different types. Such a molecular model leads to a polydisperse system of chain-associates with different chain lengths.
Wertheim’s efficient way to approach highly directional forces was successfully applied in perturbation theories in the year 1990 by Chapman et al.[33]. Chapman [34] restated 5
sized particle species with a different number of interaction sites on each. This allowed to define monodisperse molecular configurations of irreversibly associated hard-spheres. the study of Chapman, Jackson and Gubbins then led to the Statistical Associating Fluid Theory (SAFT), as an equation-of-state framework. From that point on the SAFT framework was used to develop various SAFT models. A number of reviews provide a condense overview of developments and applications in this area of research[35–39].
Structure of this thesis
In Chapter 1 this work gives a detailed review of Wertheim’s TPT. It is obviously uncommon for a thesis to provide a rather elaborate introduction to a theory that is more than 30 years old. However, for Wertheim’s TPT I have come to realize that the theory is well documented in all aspects that target applications of the theory to various types of fluids in homogeneous or inhomogeneous systems, whereas the theory is even among experts in the field much less well understood in details of how Michael S. Wertheim derived the theory. Chapter 1 is therefore providing a pedagogical review that is meant to clarify Wertheim’s derivation of TPT and facilitate further development in fluid theory. First a brief introduction is given to the basis in Statistical Thermodynamics. Then the relevant concepts of graph theory are illustrated and the Statistical Mechanical equations are reformulated in terms of graph sums for a general molecular model. The graph-theoretical approach simplifies the treat-ment of otherwise rather complex equations. Next, Wertheim’s molecular model is applied and equations are developed according the multiple-density formalism. Concise approximate equations are derived by introducing Wertheim’s assumptions and simplifications, like the ‘single chain approximation’. Based on these results, the final theories are formulated, in-cluding TPT1, TPT2. This work goes beyond the original work of Wertheim in providing a general expression of a n-th order theory TPTn. In the second Chapter a third-order expansion is developed and evaluated. Three- and four-particle distribution functions of the reference (hard-sphere) fluid for chain-like hard-sphere configurations are determined by Monte Carlo simulations and are introduced into the theory. The third Chapter extends the Wertheim’s theory by going beyond the single-chain approximation. A closed expression of the dimer-dimer graph sum is included in the theoretical development and is graphically derived in terms of 2-, 3- and 4-particle distribution functions of the reference fluid. All results are compared to simulation data for hard-chains with chain lengths 4, 8, 16 and 64. The low-density limit is assessed by using the second virial coefficient of flexible hard-chains. The so derived theory improves TPT mainly in the low-density limit, where the single chain approximation of the original TPT is a rather coarse approximation.
An analysis and detailed pedagogical
review of Wertheim’s Thermodynamic
Perturbation Theory
The content of this chapter is a literal quote of the publication
W. Zmpitas, J. Gross. Detailed pedagogical review and analysis of Wertheim’s thermody-namic perturbation theory. Fluid Phase Equilibria, 428:121-152, 2016
Additions or deletions compared to the published work are marked with angular brackets. The Thermodynamic Perturbation Theory (TPT) developed by Michael S. Wertheim in the years 1984 to 1986 provides a powerful framework for modeling directional attractive interac-tions. TPT is successful, for example, in describing fluids with hydrogen-bonding interactions and it leads to accurate models of non-spherical chain fluids. The theory of Wertheim is el-egant in many ways, but his original work is not easy to read and does not effortlessly reveal all of its beautiful features. Our study aims at providing a review of the graph-theoretical development that led Wertheim to the TPT. Our motive is to (1) develop appreciation for the ingenious development of Wertheim, (2) make the assumptions behind TPT transparent, and (3) thereby possibly facilitate further development in the theory of fluids. This review is mainly pedagogical and we intend to alleviate the lack of scientific novelty to a small extent by formulating the TPT of arbitrary order.
1.1
Introduction
[...] This review first formulates essential properties of thermodynamics and statistical me-chanics in the usual single-density formalism (section 1.2). In section 1.3, we introduce the basics of graph theory. We give graphical representations of the thermodynamic properties described in section 1.2. We only include those aspects of graph theory that are needed for the purpose of reproducing the development of Wertheim’s TPT. One essential step is the reduction of graph series, referred to as topological reduction, which we illustrate in this section. The molecular model is shown in section 1.4. Chapter 1.5 introduces the multiple-density formalism for the thermodynamic properties of section 1.2. The multiple-multiple-density formalism is exact at this stage. The subsequent sections present varying levels of approxi-mation leading to the TPT of different order. We thereby consider fluids with a single site and subsequently fluids with two sites (section 1.6). Finally, subsection 1.6.5 reviews the generalization to an arbitrary number of sites according to the SAFT of Chapman et al.
1.2
Classical Thermodynamics
Perturbation theories are usually formulated in terms of Helmholtz energy A, which char-acterizes a system with fixed number N of particles, volume V and temperature T . By defining a well known reference fluid with Helmholtz energy AR, the influence of any further
interaction type between reference particles, referred to as perturbation, can be expressed as ∆A = A − AR. In this review, properties of the reference system are denoted by subscript
R. By neglecting cross interactions between multiple present perturbations, the Helmholtz energy difference is obtained as
A − AR= ∆A1+ ∆A2 + ... (1.1)
In the more general case of inhomogeneous systems, the Helmholtz energy A[ρ] is a functional of the entire density profile ρ(r). In this notation, we don’t make the dependence of the Helmholtz energy on temperature explicit. We define r as a generalized positional coordinate. From density functional theory [40], one gets for a pure fluid
A[ρ] = Z
ρ(r) (µ − Φ(r)) dr + Ω[ρ] (1.2) relating the intrinsic part of the Helmholtz energy A[ρ] and the grand potential Ω[ρ]. Fur-ther, µ is the fluid’s chemical potential (including external forces), and Φ(r) describes the interaction of an external field with particles located at r. We emphasize that the term ’in-trinsic’ indicates that the intrinsic Helmholtz energy A[ρ] only contains the potential energy caused by inter-particle interactions (as well as the intra-molecular interactions, for the case of more complex molecular fluids), whereas the potential energy of the external field is not covered.
Strictly speaking, it is justified to speak of A[ρ] as the intrinsic Helmholtz energy and Ω[ρ] as the grand potential only for the system’s equilibrium density profile ρ(r) = ρeq(r).
In the equilibrium state, functional Ω[ρ] is minimized regarding ρ(r) and its minimum value equals the grand potential Ω characterizing a system with volume V , temperature T and the fluid’s chemical potential µ. The density profile is thus no longer a degree of freedom and is determined by the functional derivative of eq. (1.2)
δA[ρ] δρ(r) ρeq(r) = µ − Φ(r) (1.3) The calculus of functionals is inevitable for inhomogeneous systems and we refer to literature[41].
For an equilibrated system, the grand canonical (V, T, µ) ensemble provides a relation between the grand potential and microscopic properties of a fluid with
Ω = −kT ln Ξ (1.4) where k is Boltzmann’s constant. The grand canonical partition function
Ξ = X N ≥0 " exp(βN µ) N !Λ3N Z exp − βUN({r}) exp −β N X m=1 Φ(rm) ! N Y m=1 drm # (1.5) represents the sum of all possible microstates for a pure fluid in a system with fixed (V, T, µ) and external potential Φ(r). In this equation and henceforth, it is convenient to introduce the reciprocal temperature β = 1/(kT ). The outer sum in eq. (1.5) runs over the particle number N from 0 to ∞. The de Broglie wavelength Λ captures the total translational kinetic energy of particles. In our notation, Λ is defined to also contain the rotational and other intra-molecular degrees of freedom.
The system’s potential energy UN({r}) arises from interactions between particles, that
are assumed only pairwise additive, reading UN({r}) =
N
X
m>n
u(rmn) (1.6)
where u(rmn) is the pair potential between the particles m and n, and rmn = rn− rm defines
the arrangement of the particles (or molecules) towards one another. For spherically sym-metric pair potentials one simply has u(rmn) = u(|rmn|) with |rmn| as the core-core distance.
For a concise notation, we use {r} = {r1, r2, ...} to express the positional coordinates of all
particles.
In what follows it is useful to restate eq. (1.5) by factoring particle-specific quantities into the fugacity
We note that the so-defined fugacity differs from the fugacity f (r)/p = z(r)/ρ(r) used commonly in engineering applications. We introduce the notation Ξ[z] to explicitly indicate the functional dependence of the partition function regarding the fugacity profile z(r). By inserting eq. (1.6) and (1.7) into eq. (1.5) we obtain
Ξ[z] = X N ≥0 " Z 1 N ! N Y m=1 z(rm) ! exp −β N X m>n u(rmn) ! N Y m=1 drm # (1.8) For a system with N indistinguishable particles, we obtain the averaged probability density ρ(r) of finding particles at location r according to
ρ(r) = * N X m=1 δ(|r − rm|) + Ξ[z] (1.9) where the bracket h...iΞ[z] denotes the grand canonical ensemble average. The Dirac delta
function δ(|r − rm|) is defined as usual; when multiplied onto an integrand, the Dirac delta
function extracts the value of the integrand for rm = r. The probability density of finding s
particles with configuration (r0, r00, ..., rs) is analogously ρ(s)(r0, r00, ..., rs) = * N X m=1 N X n6=m ... N X k6=m,n,... δ(|r0− rm|)δ(|r00− rn|) · ... · δ(|rs− rk|) + Ξ[z] (1.10) From eq. (1.8) and (1.10), a general expression of the s-particle probability function is deduced as ρ(s)(r0, r00, ..., rs) = 1 Ξ[z] δsΞ[z] Qs i=1δln z(ri) (1.11) where ρ(s)(r0, r00, ..., rs) is proportional to the s-fold functional derivative of the grand canon-ical partition function Ξ[z] with respect to logarithmic fugacity.
To derive an expression for the logarithmic grand canonical partition function ln Ξ[z], as needed in eq. (1.4), we start with its variation
δln Ξ[z] = Z
δln Ξ[z]
δln z(r) δln z(r)dr (1.12) and replace the functional derivative on the right-hand side with ρ(r), as given in eq. (1.11) for the case of s = 1, to get
ln Ξ[ρ] = − Z ρ(r) lnρ(r) z(r)− 1 dr + c(0)[ρ] (1.13) with the definition
δc(0)[ρ] = Z
lnρ(r)
We note that ln Ξ[ρ] in eq. (1.13) is now expressed in terms of the density profile ρ(r). The transition from fugacity z(r) to density ρ(r) as a variable of the grand canonical partition function will be important in the following section, where for a graph theoretical formalism the transition is performed by a topological reduction. Functional c(0)[ρ] is defined by its variation in eq. (1.14). Its functional derivative with respect to density is thus identified to
δc(0)[ρ] δρ(r) = ln
ρ(r)
z(r) ≡ c(r) (1.15) with c(r) being the single-particle direct correlation function.
By inserting eq. (1.4), (1.7) and (1.13) into eq. (1.2), the intrinsic Helmholtz energy becomes
βA[ρ] = Z
ρ(r) ln Λ3ρ(r) − 1 dr − c(0)[ρ] (1.16) where the first term represents the ideal gas contribution and the functional (−c(0)[ρ]) is equal to the intrinsic residual Helmholtz energy. Eq. (1.16) clearly reveals the meaning of the single-particle direct correlation function c(r), as defined in eq. (1.15), as the intrinsic residual chemical potential.
The Helmholtz energy difference between a fluid and the reference fluid is β(A[ρ] − AR[ρ]) = −
c(0)[ρ] − c(0)R [ρ] (1.17) We emphasize, however, that in this perturbational formalism, the target fluid and reference fluid have the same ideal gas contribution.
A wonderful feature of the thermodynamic formalism of this section is that all equations also hold for more complex molecules. We so far limited all arguments and explanations to r being a spacial coordinate. We can alternatively define r as a generalized positional coordinate. We let r define a position in space, an orientational vector and a configurational vector, specifying the coordinates of all interaction sites of a molecule. The density in this formalism ˆρ(r) represents a probability density of finding molecules in condition r, i.e. in a certain position, orientation, and conformation. Let us decompose r into r = (˜r, ω), where ˜
r simply represents the positional vector, pointing, say, to the center of the molecules’ mass and ω denotes the orientational and configurational vector. Then ρ(˜r) = R ρ(r)dω denotesˆ the number density of molecules. With this generalization in mind, we find that all equations of this section remain valid, with only one limitation. The de Broglie wavelength is for a generalized framework defined to capture rotational and conformational degrees of freedom. Doing so is expedient only if the probability distribution of intramolecular degrees of freedom is (assumed) density independent. Otherwise, one hides a density dependent function in the de Broglie wavelength, which is counteracting clarity and intuition.
1.3
Basics of graph theory
In this section, we summarize basic concepts of graph theory. For more depth, we refer to good reviews available in literature[17, 24, 25]. We derive graphical representations of thermodynamic properties introduced in section 1.2. For brevity, we introduce the shorthand notation
e(mn) = exp (−βu(rmn)) (1.18)
for the Boltzmann factor. We start with eq. (1.8), giving the grand canonical partition function Ξ[z] in terms of fugacity z(r). The first four terms of the infinite sum read
Ξ[z] = 1 + Z z(r)dr + 1 2! Z z(r1)z(r2)e(12)dr1dr2 + 1 3! Z z(r1)z(r2)z(r3)e(12)e(13)e(23)dr1dr2dr3+ ... (1.19) An explicit formulation beyond the first few terms is an elaborate task and it is very helpful to introduce a simple graphical representation of each term of the sum. The fugacity z(r) of a particle, which is located at r, is represented as a point . Points with a white background pattern are referred to as root points, because they have a fixed location. A colored background pattern, i.e. , indicates an integration over the particle location. Points with a colored background pattern are called field points. In fact, the graphical illustration is not only restricted to the fugacity z(r), but can also be applied to other properties, such as density ρ(r). To distinguish between properties, we will use different background colors of the points and refer to them for example as z-root points and z-field points in case of fugacities z(r). The Boltzmann factor e(mn), which connects two interacting particles m and n, is illustrated with a ’zigzag’ line and will be referred to as e-bond. As a consequence, eq. (1.19) can be illustrated as a sum of graphs, in which each z-field point is connected by e-bonds to any other z-field point within the particular graph. The first six terms of the sum then turn into
Ξ[z] = 1 + + 1 2! + 1 3! + 1 4! + 1 5! + ... (1.20) The prefactor of a graph corrects for the permutability of particles and prevents any config-uration of, say, two particles (rm, rn) to be revisited as the (identical) configuration (rn, rm).
In the language of graph theory, by labelling each particle the denominator of a prefactor represents the number of topologically identical graphs that are obtained by swapping field point labels. We will call this denominator the symmetry number σk of graph k.
For the further development, eq. (1.20) is recast into a more convenient form by intro-ducing the Mayer f -function [4]
In graphical representation, each e-bond between two z-field points is replaced with e(nm) = 1 + f (nm)
= + (1.22)
Here, the straight line connecting two z-field points is referred to as f -bond. The graphical substitution of each e-bond in eq. (1.20) with (1.22) leads to unconnected graphs, consisting of two or more subgraphs that are not connected by f -bonds. A missing f -bond between different subgraphs means, that their integrations can be independently performed. Thus, unconnected graphs can be illustrated as product of their subgraphs. As an example, we show that 1 2 = × 1 2 1 2 Z z(r1)z(r2)z(r3)f (23)dr1dr2dr3 = Z z(r1)dr1× 1 2 Z z(r2)z(r3)f (23)dr2dr3 (1.23) As a consequence, the infinite graph sum, eq. (1.20), can be written as a product of infinite sums, where each sum contains only one type of subgraph k with its symmetry number σk,
as Ξ[z] = " 1 0!+ 1 1! + 1 2! 2 + 1 3! 3 + ... # × " 1 0!+ 1 1! 1 2 + 1 2! 1 2 2 + 1 3! 1 2 3 + ... # × " 1 0!+ 1 1! 1 2 ! + 1 2! 1 2 !2 + 1 3! 1 2 !3 + ... # × ... (1.24) Each sum in eq. (1.24) corresponds to an infinite series for which
∞ X i=0 1 i! Ik σk i = exp Ik σk (1.25) with Ik representing subgraph k. Each subgraph is a connected graph, i.e. each z-field point
is bound to at least one other z-field point with an f -bond. After substituting each sum in eq. (1.24) according to eq. (1.25), the grand canonical partition function reads
Ξ[z] = exp ∞ X k=0 Ik σk ! (1.26)
By applying the natural logarithm to eq. (1.26), we obtain the infinite sum of connected graphs ln Ξ[z] = +1 2 + 1 2 + 1 6 + 1 2 + 1 6 + ... (1.27) The number density ρ(r) can now be determined using eq. (1.11) with s = 1. We obtain
ρ(r) = z(r)δ ln Ξ[z] δz(r) = + +1 2 + + 1 2 + + + ... (1.28) The functional differentiation of eq. (1.27) (with respect to fugacity z(r)) and multiplication by z(r) yields a sum of connected graphs with one z-root point in each graph. Except for the monomer graph , any other graph in eq. (1.27) contains more than one z-field point. To differentiate a product of fugacities, the product rule must be applied. In order to better illustrate the mathematical procedure that corresponds to the differentiation in graph theory, we take the second graph of eq. (1.27) as an example, for which
z(r) δ δz(r) 1 2 = z(r) δ δz(r) 1 2 Z z(r1)z(r2)f (12)dr1dr2 = z(r) Z z(r0)f (r0, r)dr0 (1.29) After differentiating and summing up identical graphs, the prefactors in eq. (1.28) still ac-count for the symmetry of the graphs. By factoring out z(r) (i.e. the z-root point), eq. (1.28) is restated as ρ(r) = z(r) × + +1 2 + + 1 2 + + ... ! (1.30) The dashed root point in each graph represents only the fixed location r of the factored out fugacity and equals the value of 1.
The expression of ρ(r) allows us to graphically illustrate the functional c(0)[ρ], as
de-fined in eq. (1.14). We start by determining the intrinsic residual chemical potential c(r), eq. (1.15). By dividing the number density ρ(r), eq. (1.30), by the fugacity z(r) (or func-tionally differentiating eq. (1.27) with respect to z(r), as shown in eq. (1.11) for s=1), we obtain ρ(r) z(r) = + + 1 2 + + 1 2 + + + ... (1.31)
Some graphs in eq. (1.31) (in particular the third and seventh graph) are unconnected graphs and can be decomposed into a product of their subgraphs. This is, from a graphic point of view, not immediately obvious. The missing volume integration over a root point’s location leads to independent subparts, that are not directly f -bond -connected to each other. To identify these types of graphs, one introduces the term ’articulation point ’. Generally, points, that upon deletion cause a graph to become an unconnected graph, i.e. a graph consisting of unconnected subgraphs, are referred to as articulation points. We differ between articulation field points and (dashed) articulation root points - both types are relevant in graph theory. Graphs with a dashed articulation root point in eq. (1.31) can be expressed as a product of their unconnected subgraphs, which can independently be integrated over. The principle is shown for the third graph on the right-hand side, with
1 2 = 1 2 2 1 2 Z z(r2)z(r3)f (12)f (13)dr2dr3 = 1 2 Z z(r0)f (r0, r)dr0 2 (1.32) As a consequence and similar to the prescription leading to eq. (1.27), we can rewrite eq. (1.31) as a product of sums, where each sum contains only one type of subgraph. We obtain ρ(r) z(r) = " 1 0!+ 1 1! ( ) + 1 2! ( ) 2 + 1 3! ( ) 3 + ... # × " 1 0!+ 1 1! ! + 1 2! !2 + 1 3! !3 + ... # (1.33) × " 1 0!+ 1 1! 1 2 ! + 1 2! 1 2 !2 + 1 3! 1 2 !3 + ... # × ... which is again a product of infinite sums, whose limit value is given with eq. (1.25). We see that separating graphs bound to the dashed articulation root point according to the scheme of eq. (1.32), grouping them as in eq. (1.33) and taking the natural logarithm leads to the intrinsic residual chemical potential c(r), eq. (1.15), as
c(r) = + + 1 2 + + 1 2 + 1 2 + ... (1.34) with solely connected graphs without dashed articulation root point.
The intrinsic residual chemical potential c(r) is needed in order to determine the residual Helmholtz energy c(0)[ρ] by integration of eq. (1.14) using eq. (1.15). We note, however, that
we strive for c(0)[ρ] formulated as a functional of the number density ρ(r). The function c(r), as derived in eq. (1.34), is still expressed in terms of fugacities z(r), i.e. each field point in a graph is a z-field point. To accomplish the transition from fugacity z(r) to number density ρ(r), a topological reduction is applied by substituting a particular, infinite sum of graphs each consisting of one z-root point and z-field points with ρ(r), according to the relation given in eq. (1.30).
The graph sum in eq. (1.34) contains articulation z-field points, whereas we already elim-inated dashed articulation root points in c(r). We first identify graphs without articulation field points which are referred to as irreducible graphs. Examples for irreducible graphs are the first, third and sixth graph in eq. (1.34). Graphs with at least one articulation field point are on the other hand called reducible graphs. We find reducible graphs in c(r) (the second, forth and fifth) that have an articulation z-field point that upon deletion, leads to unconnected subgraphs of which one is equal to the irreducible graph . By summing up all graphs of this type (including the corresponding irreducible graph itself), we obtain a reduction of the form
+ + + 1
2 + 1
2 + ... = (1.35) with the ρ-field point
= + + + 1 2 + 1 2 + ... Z ρ(r)dr = Z z(r) × + + +1 2 + 1 2 + ... ! dr (1.36) For clarity, the articulation z-field point, on which the topological reduction is applied, is highlighted by a green background pattern, i.e. , in eq. (1.35) and (1.36). The graph sum in eq. (1.36) equals the number density ρ(r), as derived in eq. (1.30), but integrated over the position of the z-root point. This is equal to replacing the z-root point by a z-field point. The topological reduction, illustrated above, reduces an infinite sum of reducible graphs into one single irreducible graph by eliminating articulation z-field points and introducing ρ-field points. We emphasize that it is also possible to apply a topological reduction on other articulation z-field points as on the one chosen in eq. (1.35) leading to reducible graphs consisting of both z- and ρ-field points. However, these graphs are always further reducible to irreducible graphs by eliminating the remaining articulation z-field points. By applying this procedure illustrated in eq. (1.35) and (1.36) based on each irreducible graph of eq. (1.34), we obtain c(r) = +1 2 + 1 2 + 1 2 + 1 2 + 1 6 + ... (1.37)
As a consequence of the applied topological reduction, c(r) is now (graphically and analyti-cally) expressed in terms of the number density ρ(r).
The intrinsic chemical potential c(r) equals the functional derivative of the residual Helmholtz energy c(0)[ρ] with respect to the number density ρ(r), as shown in eq. (1.15). Since we have already obtained an expression for c(r) with eq. (1.37), the functional deriva-tive can graphically be reversed to obtain an expression for c(0)[ρ] as
c(0)[ρ] = 1 2 + 1 6 + 1 8 + 1 4 + 1 24 + ... (1.38) It is easy to show, that functionally differentiating eq. (1.38) with respect to number density ρ(r) yields eq. (1.37). A graphical illustration of the logarithmic grand canonical partition function ln Ξ[ρ] in terms of ρ-field points is derived by inserting eq. (1.37) and (1.38) into eq. (1.13). The first few graphs are
ln Ξ[ρ] = −1 2 − 1 3 − 1 8 " 3 − 6 − # − ... (1.39) and can directly be identified as the graphical representations of second, third and fourth virial coefficient of the virial expansion [4, 42].
1.4
Molecular Model
All derivations of the previous sections 1.2 and 1.3 hold for arbitrary pair potentials u(rmn)
between two particles m and n. For molecular fluids, rmn carries a distance measure but
also defines the orientation and conformation of both molecules. The system’s overall inter-molecular potential energy reads
UN({r}) = N
X
m>n
u(rmn) (1.40)
Perturbational approaches are based on the division of a target-fluid’s intermolecular pair potential u(rmn) into a reference part (R) and a perturbation contribution (P ). In many
cases, the reference fluid covers the repulsive part of the potential and the perturbation contribution captures the attractive portion of the interaction potential. Such a division is useful (1) if the properties of the repulsive fluid are known to sufficient accuracy and (2) if the structure of a fluid (i.e. g(2)(r
1, r2)) is dominated by the repulsive interactions. The
decomposition of the pair potential is thus done according to
u(rmn) = uR(rmn) + uP(rmn) (1.41)
The corresponding Mayer f -function f (mn) of eq. (1.41) can be written as
f (mn) = fR(mn) + eR(mn)fP(mn) (1.42)
with fR(mn) = eR(mn) − 1 and fP(mn) = eP(mn) − 1, using the shorthand notation of
eq. (1.18) and (1.21).
Repulsive interactions are in the simplest case of spherically symmetric type, where the repulsive potential uR(rmn) solely depends on the distance rmn = |rn− rm| between the
spheres’ centers, i.e. uR(rmn) = uR(rmn). The repulsion can be described by a soft-core
potential, i.e. allowing partial core overlap, or by a hard-core potential, i.e. prohibiting any core overlap.
To account for attractive interactions, Wertheim uses a set Γ = {A, B, C, ...} of acen-trically located sites within each particle. The attractive potential uAB(rmAnB) between
site A of particle m and site B of particle n can generally have any functional form, with |rmAnB| = |ˆrn+ dB(Ωn) − ˆrm − dA(Ωm)| as the distance between the sites A and B, and
dA(Ωm) and dB(Ωn) as their position vectors originating in the centers of the particles m and
n with orientations Ωm and Ωn. We use ˆrn for the core position of a particle. However, we
will see in a later section that Wertheim, in his development, took advantage of a sufficiently short-ranged, spherical potential around the acentrically positioned sites, which allowed him to introduce steric effects efficiently.
The overall pair potential in eq. (1.41) between two identical particles m and n with the set Γ of sites on each particle thus reads
u(rmn) = uR(ˆrmn) + X A∈Γ X B∈Γ uAB(rmAnB) (1.43)
One gets for the Mayer f -function, cf. eq. (1.42), as f (mn) = fR(mn) + eR(mn) " Y A∈Γ Y B∈Γ (fAB(mn) + 1) − 1 # (1.44) where (mn) is used for compact notation, i.e. we don’t make explicit that coordinates rmn
contain both, the vector pointing from one to the other particle core, ˆrmn, as well as the
orientations Ωm and Ωn, defining rmAnB.
The repulsive prefactor eR(mn) of the bracket in eq. (1.44) suppresses attractive fAB
-bond contributions between two sites A and B of particles m and n whenever cores of the particles overlap. In the case of a hard-core potential, core overlap means that eR(mn) goes
to zero and the value of the bracket is irrelevant. With a graph theoretical mindset, the second term allows the introduction of highly directional attractive interactions by the right choice of the range of the site-site potential, uAB(rmAnB), and the geometric location of the
sites within the particles.
In his first two publications [20, 21], Wertheim considers the simple case of particles with one single site A. Then, eq. (1.44) reads
f (mn) = fR(mn) + eR(mn)fAA(mn) (1.45)
which we illustrate with
= A A + A A (1.46)
between two particles. To distinguish between the different interaction functions fR(mn),
eR(mn) and fAA(mn), we choose to add labelled, small black nodes within each field point
and root point illustrating the acentric sites. In the remainder of this review, we adopt Wertheim’s notation and refer to the field point’s and root point’s shell as hyperpoint and to the attractive sites simply as attractive sites. We continue representing f -bonds by straight lines and e-bonds by ’zigzag’ lines , as introduced in section 1.3. Attractive fAB-bonds
connect two attractive sites A and B of distinct particles, whereas repulsive fR-bonds and
eR-bonds only connect the particles’ hyperpoints. As eq. (1.44) demands, hyperpoints of two
attractively connected particles are always connected with an eR-bond to each other.
In Wertheim’s third paper [22], an exact TPT for particles with multiple attractive sites is derived. The explicit case of polymerisation for two attractive sites is shown in his fourth paper [23]. For particles with two attractive sites A and B, eq. (1.44) takes the more
elaborate, graphical form = A B B A + A B B A + A B B A + A B B A + A B B A + A B B A + A B B A + A B B A + A B B A + A B B A + A B B A (1.47) + A B B A + A B B A + A B B A + A B B A + A B B A
illustrating each possible bonding state between two particles. In this sum, we grouped terms with the same number of attractive bonds in the same column. Bonding constraints, however, can prohibit some of these bonding states. E.g. assuming that attractive sites of the same type are not interacting with each other, i.e uAA(rmAnA) = uBB(rmBnB) = 0, and
defining the attractive potential range and the attractive sites’ geometry such that double bonds between particles are prohibited, eq. (1.47) simplifies to
f (mn) = fR(mn) + eR(mn)fAB(mn) + eR(mn)fBA(mn) = A B B A + A B B A + A B B A (1.48) Introducing steric limitations at an early stage eliminates a multitude of graph types and simplifies the further treatment.
For particles with more than two attractive sites, the graphical illustration is derived analogously, but results in a more elaborate treatment. For reasons of clarity and con-ciseness, we show in the following sections of this review the graphical derivations of the thermodynamic properties in terms of particles with two attractive sites and, moreover, we adopt the assumptions leading to eq. (1.48), i.e no interactions between attractive sites of the same type and no double bonded particles. The more general case, depicted in eq. (1.47) leads to an increased complexity without providing much more insight.
1.5
Wertheim’s TPT
Wertheim pursues the idea of introducing geometry of interactions at an early stage for fluids with highly directional forces, resulting in a multiple-density formalism and allowing an early consideration of steric effects. For this purpose, each f -bond in the logarithmic grand canonical function ln Ξ[z], eq. (1.27), is substituted by eq. (1.48). We obtain
ln Ξ[z] = A B + A B B A + 1 2 A B A B B A +1 2 A B A B B A +1 2 A B A B B A + ... + A B A B B A + A B A B B A +1 3 A B A B B A + B A A B B A A B + B A A B B A A B + ... +1 2 A B B A + 1 2 A B A B B A +1 6 A B A B B A +1 8 B A A B B A A B +1 4 B A A B B A A B + ... + A B A B B A + A B A B B A + A B A B B A +1 2 B A A B B A A B + B A A B B A A B + ... (1.49) Eq. (1.49) can be grouped in three subtotals. The first subtotal {S} corresponds to the first two rows of the sum. Each particle in {S} is connected to any other particle within the particular graph by a path of fAB-bonds (and consequently eR-bonds). These graphs will
be referred to as pure s-mer graphs, with s indicating the number of attractively bonded particles. The monomer (s=1) graph is also included in this first subtotal. The second subtotal {M } is illustrated in the third row and contains graphs with solely repulsively interacting monomers, i.e. the monomers’ hyperpoints are connected with fR-bonds. Each
graph in the third subtotal {R}, shown in the last row, consists of repulsively interacting s-mers with at least one of them an (s ≥ 2)-mer. E.g. the fourth and fifth graph of the last row show two dimers (s=2) being connected with fR-bonds.
The pure s-mer graph sum {S} can further be divided into the bare s-mer graph subtotal {B} lacking of any fR-bonds and the hindered s-mer graph subtotal {H} containing at least
one fR-bond, as {B} = A B + A B B A + 1 2 A B A B B A + 1 2 A B A B B A + A B A B B A + 1 3 A B A B B A + B A A B B A A B + ... (1.50) {H} = 1 2 A B A B B A +1 2 A B A B B A + A B A B B A + B A A B B A A B + ... (1.51)
Now let us take only the bare s-mer graph sum {B}, eq. (1.50). If the removal of a z-field hyperpoint (the shell of a z-z-field point, i.e. the attractive sites remain) along with its attached eR-bonds leads to an unconnected graph, the so-obtained subgraphs are said to be
constraint- connected at the attractive sites of the removed z-field hyperpoint. Constraint-connected subgraphs are only present for particles with two or more attractive sites. The last graph in eq. (1.50) is an example for such a case: the top right particle and the two bottom particles are constraint- connected to the attractive sites of the top left particle. Further, by removing all z-field hyperpoints and eR-bonds, only attractive sites and attractive fAB-bonds
are left. The chosen graph then turns to
B A A B B A A B −→ B A A B B A A B (1.52) Attractive sites, which are uninterruptedly connected by a path of fAB-bonds, are then said
to be located in the same bond-connected network. In the case of particles with one single attractive site, each graph in {S} contains exactly one bond-connected network, i.e. each attractive site within a graph is located in the same bond-connected network. To account for steric effects, Wertheim inserts eR-bonds between each two z-field hyperpoint of not
attractively connected particles within a bond-connected network (the z-field hyperpoint of two attractively connected particles are already connected with an eR-bond, as shown in the
scheme (1.52)). Introducing eR-bonds is possible, because for each graph in {B}, we find
similar graphs in {H} additionally including fR-bonds between not attractively interacting
particles. This is originally caused by the substitution of each e-bond according to eq. (1.21) and (1.44). These types of graphs can be combined according to the following scheme
B A A B B A A B + B A A B B A A B = B A A B B A A B (1.53) illustrating the combination of the last graph of eq. (1.50) and (1.51). We emphasize that there are also graphs in {H}, which are not combined with graphs of {B}. Such graphs in-clude fR-bonds connecting hyperpoints, which belong to different bond-connected networks.
The third graph in eq. (1.51) is an example for this case.
Introducing eR-bonds into bond-connected networks is an important step, because
site-specific bonding constraints can then be easily applied. For example, if an attractive site can’t be bonded to multiple attractive sites of distinct particles at the same time, as being the case for attractive site A of the bottom left particle in eq. (1.53), because of a sufficiently short interaction range of the site-site potential, any bond-connected network with more than two particles is prohibited for a repulsive hard-core potential. At least one eR-bond
will then be zero, due to overlap of the particles’ cores. More specific, the integral over all geometrical arrangements of such a bonding state will be zero, which in turn means the
respective graph is zero. The graph on the right-hand side of eq. (1.53) is an example for a graph vanishing under these conditions, whereas both original graphs on the left-hand side are non-vanishing. Generally, the insertion of eR-bonds can also be introduced between
constraint-connected subgraphs by additionally considering the remaining graphs in {H}. However, doing so is not helpful, because one can not easily make use of steric constraints that would eliminate graphs.
The bare s-mer graph subtotal {B} was only chosen to comprehensively illustrate this procedure. The insertion of eR-bonds is applied to any other graph with bond-connected
networks in eq. (1.49). Before proceeding we note that in the remainder of this section, we will not actually set any graph, such as the graph on the right hand side of eq. (1.53), to zero, because in this section we do not wish to limit the theory to a certain class of molecular models. Restricting the type of molecular model and introducing approximations is reserved for the following section.
We now turn to the number density, which in the formalism of Wertheim is decomposed into several contributions. The number density ρ(r) is determined by eq. (1.11). The dif-ference to the graph treatment of section 1.3 is that particles are now represented with two attractive sites. By introducing eR-bonds into bond-connected networks in eq. (1.49)
and functionally differentiating the resulting graphs with respect to logarithmic fugacity, we obtain ρ(r) = A B + A B B A + 1 2 A B A B B A + A B A B B A +1 2 A B A B B A + A B A B B A + ... + A B B A + A B A B B A + A B A B B A + A B A B B A +1 2 A B A B B A + ... + A B B A + A B A B B A + A B A B B A + A B A B B A + A B A B B A + ... + A B A B B A + A B A B B A + A B A B B A +... (1.54) Wertheim realized that it is possible to decompose ρ(r) into, in this case, four subtotals of graphs, where the z-root point is bonded only repulsively, or only attractively at its attractive site A, or only attractively at its attractive site B, or simultaneously at both attractive sites A and B. We arranged these four subtotals as lines in eq. (1.54): The z-root point in the first row is only repulsively connected to any z-field point within the graphs, whereby the z-field points may very well be additionally connected attractively to each other (see for example the last graph in the first row). All graphs of the second row have a z-root point with attractive site A bonded to one or more other attractive sites of distinct z-field point,
while additionally connected with fR-bonds to other z-field points. The same holds for the
third row, but now for bonded attractive site B. The last row contains z-root points bonded at both, attractive site A and site B. Depending on the bonding state, Wertheim introduces the multiple- density formalism
ρ(r) =X
γ⊂Γ
ργ(r) (1.55)
= ρ0(r) + ρA(r) + ρB(r) + ρAB(r) (1.56)
where the index γ = {0, A, B, AB} in density ργ(r) indicates the set of the z-root point’s
bonded attractive sites. For γ = 0, the particle has no bonded attractive sites, corresponding to the first line of eq. (1.54). The formalism can analogously be applied to particles with more attractive sites, leading to a general division of the number density ρ(r) into
ρ(r) = ρ0(r) + X A∈Γ ρA(r) + 1 2 X A∈Γ X B∈Γ\{A} ρAB(r) + 1 6 X A∈Γ X B∈Γ\{A} X C∈Γ\{A,B} ρABC(r) + ... (1.57)
The number density ρ(r) represents the probability of finding a particle at the location r regardless of its bonding state, whereas the single density ργ(r) represents the probability
of finding a particle with set γ = {A, B, ...} of bonded attractive sites at r. We emphasize that the analytical division of the number density, as shown in eq. (1.57), is universally valid for any applied bonding constraint. Even the case, where an attractive site A is not permitted to bond at all, is included, because each density ργ(r) with A ∈ γ simply turns to
zero. Graphically spoken, each graph with a bonded attractive site A turns to zero, because fAB(nm) = 0 for each B ∈ Γ\{A}.
Let’s proceed with eq. (1.11) and (1.12) describing the variation of the logarithmic grand canonical partition function. The variational equation can be extended to a multiple-density formalism, according to δln Ξ[z] = Z ρ(r)δ ln z(r)dr = δ Z ρ(r)dr − Z ρ0(r)δ ρ(r) ρ0(r) dr − Z ρ(r)δ lnρ0(r) z(r)dr (1.58) In the last term of eq. (1.58), Wertheim defines the function c0(r) as
c0(r) ≡ ln
ρ0(r)
z(r) (1.59) resembling a residual chemical potential (or single-particle direct correlation function) similar to eq. (1.15). This quantity should not be confused with the residual Helmholtz energy c(0)[ρ],
limited to graphs contained in the density contribution ρ0(r). Its graphical representation reads c0(r) = AB BA + A B A B B A + A B A B B A + A B A B B A + B A A B B A A B +1 2 B A A B B A A B + ... (1.60) The index 0 of function c0(r) demands that all dashed root points are without bonded
attractive sites. The z-field points, however, include attractive bonds.
Replacing ρ(r) according to eq. (1.57), the second term on the right-hand side of eq. (1.58) becomes Z ρ0(r)δ ρ(r) ρ0(r) dr =X A∈Γ Z ρ0(r)δ ρA(r) ρ0(r) dr+1 2 X A∈Γ X B∈Γ\{A} Z ρ0(r)δ ρAB(r) ρ0(r) dr+ ... (1.61) To determine the variations of the density ratios δ(ργ(r)/ρ0(r)) in eq. (1.61), we use the
graphical representations of the particular densities ργ(r) derived in eq. (1.54).
For z-root points with one bonded attractive site, we consider the ratio ρA(r)/ρ0(r) as
an example. The density ρ0(r) corresponds to the first row and ρA(r) to the second row of
graphs in eq. (1.54) as ρ0(r) = AB + A B B A + 1 2 A B A B B A + A B A B B A +1 2 A B A B B A + A B A B B A + ... (1.62) ρA(r) = AB BA + A B A B B A + A B A B B A + A B A B B A +1 2 A B A B B A + 1 2 B A A B B A A B + ... (1.63) In the following, the density ratio ρA(r)/ρ0(r) will be referred to as function cA(r) ≡
ρA(r)/ρ0(r) for an attractive site A. The graph sum of ρA(r), eq. (1.63), can be
decom-posed into a product of cA(r) × ρ0(r). A topological reduction is applied through the
following pattern of graphs: some graphs (in our example the second and sixth graph of eq. (1.63)) contain an articulation z-root hyperpoint; the deletion of the z-root hyperpoint leaves unconnected subgraphs behind, of which a part was only repulsively connected to the z-root point. By factoring out those repulsively connected subgraphs (including a monomer graph A
B ), the sum in eq. (1.63) can be decomposed into a product of two graph sums,
includes a solely repulsively connected z-root point, according to ρA(r) = BA BA + A B A B B A + A B A B B A + ... | {z } cA(r) × BA + A B B A + 1 2 A B A B B A + ... | {z } ρ0(r) (1.64) By identifying the right graph sum in eq. (1.64) as ρ0(r) and reversing the factorisation,
each root point in ρA(r) is now a ρ0-root point. Such a topological reduction (here, from
z-root points to ρ0-root points) was similarly performed in eq. (1.35) for field points in the
single-density formalism. From factoring ρA(r) into cA(r) and ρ0(r), we further identify
cA(r) = BA BA + A B A B B A + A B A B B A +1 2 A B A B B A + A B A B B A + B A A B B A A B + ... (1.65) where all graphs are free of dashed articulation root points. For any other attractive site B, function cB(r) is analogously derived as eq. (1.65). In our case, i.e. particles with two
attractive sites, cB(r) reads
cB(r) = AB BA + A B A B B A + A B A B B A + A B A B B A +1 2 A B A B B A + B A A B B A A B + ... (1.66) For z-root points with two bonded attractive sites A and B, the same topological reduction is applied to obtain ρ0-root points. The division by ρ0(r) leads to
ρAB(r) ρ0(r) = A B A B B A + A B A B B A + A B A B B A + B A A B B A A B + B A A B B A A B + ... (1.67) Because subparts of graphs that are constraint-connected to the attractive sites of the dashed root point were earlier not connected with eR-bonds to each other, ργ(r)/ρ0(r) (here
ρAB(r)/ρ0(r)) can further be partitioned into a product of cψ(r) functions, with ψ ⊆ γ\{0}.
For example, the first graph and the last two graphs of eq. (1.67) arise from multiplying cA(r), eq. (1.65), by cB(r), eq. (1.66). We have already applied such treatment on similar
graphs with dashed root points in eq. (1.32). The remaining graphs, such as the second and third graph in the example above, are grouped together to define function cAB(r) consisting
of not further decomposable graphs. We emphasize that ’not further decomposable’ is not equal to ’irreducible’. Some of the remaining graphs may not be further decomposable, but reducible to irreducible graphs by a topological reduction. Eq. (1.67) is then expressed as
ρAB(r)
ρ0(r)
The decomposition for dashed root points with more than two bonded attractive sites is performed equivalently. More generally, one obtains the relation
ργ(r) ρ0(r) = X {ψ}=P (γ) Y {ψ} cψ(r) (1.69)
where P (γ) denotes the possible partitioning of the set γ of attractive sites, e.g. P (AB) = {A, B}, {AB}. Eq. (1.69) does not apply to the case γ = 0, as defined in eq. (1.59).
For the further analytical treatment, Wertheim introduces the density function σα(r)
defined as
σα(r) =
X
γ⊆α
ργ(r) (1.70)
with the order γ ⊆ α ⊆ Γ of the subsets γ and α. We notice that the unbonded state 0 resembles the empty set and thus satisfies the relation γ ⊆ α ⊆ Γ for γ and α. In the two limiting cases α = 0 and α = Γ, we obtain, according to eq. (1.57),
σ0(r) = ρ0(r), σΓ(r) = ρ(r) (1.71)
One can further derive a direct relation between density σα(r) and cψ(r)-functions by
com-bining eq. (1.69) and (1.70) into σα(r) ρ0(r) = 1 + X γ⊆α\{0} X {ψ}=P (γ) Y {ψ} cψ(r) (1.72)
To ensure the validity of eq. (1.72) for the case α = 0, the term of the sum representing the contribution γ = 0 was set to unity.
We use eq. (1.72), now for the case α = Γ, to reformulate eq. (1.61) as Z ρ0(r)δ ρ(r) ρ0(r) dr =X A∈Γ Z σΓ−A(r)δcA(r)dr + 1 2 X A∈Γ X B∈Γ\{A} Z σΓ−AB(r)δcAB(r)dr + ... = X α⊆Γ\{0} Z σΓ−α(r)δcα(r)dr (1.73)
The subscript Γ − α denotes a subset of Γ without the attractive sites, that are contained in set α. For example, let Γ = {ABCD} and α = {AD}, then we get for the subset Γ − α = {BC}.
The insertion of eq. (1.73) and (1.59) into the second and third term of eq. (1.58) leads to δln Ξ[z] = δ Z ρ(r)dr − Z X α⊆Γ σΓ−α(r)δcα(r)dr (1.74)
After reversing the variation on both sides, the logarithmic grand canonical partition function reads ln Ξ[{ρ}] = − Z ρ(r) lnρ0(r) z(r) − 1 + X α⊆Γ\{0} σΓ−α(r)cα(r) dr + c(0)[{ρ}] (1.75) in terms of the number density {ρ}, with the definitions
δc(0)[{ρ}] = Z X α⊆Γ δσΓ−α(r)cα(r)dr (1.76) and δc(0)[{ρ}] δσΓ−α(r) ≡ cα(r) (1.77)
The notation {ρ} = {ρ0, ρA, ...} in eq. (1.75) reflects the multiple-density formalism of the
expression. As announced in section 1.3, we strive for c(0)[{ρ}] to be formulated as a
func-tional of number density {ρ}. The functions cα(r) are at this point still expressed in terms
of z-field points. The transition from fugacity z(r) to number density {ρ} is performed by a topological reduction, as earlier shown in eq. (1.35). In the present case, however, with different number density contributions ργ(r) present, the topological reduction, is performed
slightly differently.
We illustrate the topological reduction considering function cA(r), eq. (1.65). The
sub-script ψ of function cψ(r) denotes the dashed root point’s bonded attractive sites. In cA(r),
attractive site A of the dashed root point is bonded. Further, the dashed root point in any function cψ(r) is not a dashed articulation root point (due to the division by ρ0(r) and the
decomposition according to eq. (1.69)). The focus is therefore on the z-field points. As al-ready explained in the derivation of eq. (1.35), we first identify irreducible graphs. Then we sum up graphs which have an articulation z-field point, that upon deletion leads to uncon-nected subgraphs of which one equals the chosen irreducible graph. In the case of Wertheim’s multiple-density formalism, however, it is more beneficial to solely delete the hyperpoints of the articulation z-field points, because it is important to know, which attractive sites of the deleted particle are bonded within the irreducible part of a graph. The topological reduc-tion based on the irreducible graph A B
B A , which is the first graph of function cA(r) in
eq. (1.65), reads A B B A + A B A B B A + A B A B B A + B A A B B A A B + B A A B B A A B + ... = A B B A (1.78)
with the σA-field point A B = A B + A B B A + 1 2 A B A B B A + 1 2 A B A B B A + A B A B B A + ... + A B B A + A B A B B A + A B A B B A + A B A B B A + ... (1.79)
The z-field point, on which the topological reduction is applied, is highlighted with a green background pattern, i.e. A
B . It is connected with its attractive site B to the dashed root
point’s attractive site A in eq. (1.78). As a result of combining graphs for introducing eR
-bonds as illustrated in eq. (1.53), bond-connected networks always represent irreducible parts within a graph. Thus, if an attractive site of an articulation z-field point is bonded within an irreducible graph, the remaining subgraphs attached to the articulation z-field point can only be bonded repulsively to its hyperpoint and/or constraint-connected to its remaining unbonded attractive sites. In our example in eq. (1.78), the subgraphs are only bonded repulsively to the hyperpoint (first row of eq. (1.79), corresponding to ρ0) and/or
constraint-connected to attractive site A of the articulation z-field point (second row, corresponding to ρA), because its attractive site B is already bonded within the irreducible graph. If we sum
these subgraphs, as shown in eq. (1.79), we receive the graph sum of the density function σA(r) = ρ0(r) + ρA(r) integrated over the particle location r with the first row in eq. (1.79)
as ρ0(r) and the second row as ρA(r). The red colored point BA is thus referred to as σA
-field point. Generally, the subscript α of σα-field points denotes the set of attractive sites,
that are not bonded within the irreducible graph. Applying the topological reduction on each articulation z-field point of a function cψ(r), we obtain irreducible graphs with σα-field
points.
As a consequence, the functional c(0)[{ρ}] can be graphically determined with the
defini-tions shown in eq. (1.76) and (1.77). It should be noted, that, again, one applies the product rule for the functional differentiation of c(0)[{ρ}] to obtain the functions cψ(r). For particles
with two attractive sites A and B, c(0)[{ρ}] can be derived as c0(r) = AB BA + 1 2 A B A B B A + A B A B B A + ... cA(r) = AB BA + A B A B B A + A B A B B A + ... cB(r) = AB BA + A B A B B A + A B A B B A + ... cAB(r) = A B A B B A + ... c(0)[{ρ}] = 1 2 A B B A + 1 6 A B A B B A + A B B A + A B A B B A + A B A B B A + ... (1.80) We obtain an expression for c(0)[{ρ}] consisting of σ
α-field points. The prefactors in this
graph sum are best understood and confirmed by performing a functional differentiation of c(0)[{ρ}] with respect to σΓ−α(r), according to eq. (1.77), and allocating the resulting graphs
to the appropriate cα(r)-function. Graphs in c(0)[{ρ}], that contain at least two particles
which differ in their bonding states, are thereby allocated to several functions cα(r).
Considering all derivations up to here, the intrinsic Helmholtz energy A[{ρ}] in eq. (1.2), together with eq. (1.4), (1.7) and (1.75), turns into
βA[{ρ}] = Z ρ(r) ln Λ3ρ0(r) − 1 + X α⊆Γ\{0} σΓ−α(r)cα(r) dr − c(0)[{ρ}] (1.81) The intrinsic Helmholtz energy of the considered class of model fluid is herewith expressed as graph sums. In order to express the intrinsic Helmholtz energy as a functional of (multiple) densities, it is convenient to express the cψ(r)-functions in terms of σα(r)-functions, rather
than vice versa. For one (A) and two attractive sites (AB), the functions cψ(r) are easily
derived, from eq. (1.72), as cA(r) = σA(r) ρ0(r) − 1, cAB(r) = σAB(r) ρ0(r) −σA(r)σB(r) ρ2 0(r) (1.82)
