Main Advances in Chemical Physics, Volume 160

Advances in Chemical Physics, Volume 160

,
4.0 / 0
How much do you like this book?
What’s the quality of the file?
Download the book for quality assessment
What’s the quality of the downloaded files?
The Advances in Chemical Physics series provides the chemical physics field with a forum for critical, authoritative evaluations of advances in every area of the discipline. This volume explores the following topics:
  • Thermodynamic Perturbation Theory for Associating Molecules
  • Path Integrals and Effective Potentials in the Study of Monatomic Fluids at Equilibrium
  • Sponteneous Symmetry Breaking in Matter Induced by Degeneracies and Pseudogeneracies
  • Mean-Field Electrostatics Beyond the Point-Charge Description
  • First Passage Processes in Cellular Biology
  • Theoretical Modeling of Vibrational Spectra and Proton Tunneling in Hydroen-Bonded Systems
Year:
2016
Edition:
1
Publisher:
Wiley
Language:
english
Pages:
360 / 371
ISBN 10:
1119165148
ISBN 13:
9781119165149
Series:
Advances in Chemical Physics
File:
PDF, 9.99 MB
Download (pdf, 9.99 MB)
1

Modeling the Effect of Damage in Composite Structures: Simplified Approaches

Year:
2015
Language:
english
File:
PDF, 11.91 MB
0 / 0
2

U-Boats at War in World War I and II: Rare Photographs from Wartime Archives

Year:
2009
Language:
english
File:
EPUB, 83.23 MB
0 / 0
Advances in Chemical Physics
Volume 160

Advances in Chemical
Physics
Volume 160

Series Editors
Stuart A. Rice
Department of Chemistry
and
The James Franck Institute,
The University of Chicago,
Chicago, IL, USA

Aaron R. Dinner
Department of Chemistry
and
The James Franck Institute,
The University of Chicago,
Chicago, IL, USA

Copyright © 2016 by John Wiley & Sons, Inc. All rights reserved
Published by John Wiley & Sons, Inc., Hoboken, New Jersey
Published simultaneously in Canada
No part of this publication may be reproduced, stored in a retrieval system, or transmitted
in any form or by any means, electronic, mechanical, photocopying, recording, scanning, or
otherwise, except as permitted under Section 107 or 108 of the 1976 United States Copyright Act,
without either the prior written permission of the Publisher, or authorization through payment
of the appropriate per‐copy fee to the Copyright Clearance Center, Inc., 222 Rosewood Drive,
Danvers, MA 01923, (978) 750‐8400, fax (978) 750‐4470, or on the web at www.copyright.com.
Requests to the Publisher for permission should be addressed to the Permissions Department, John
Wiley & Sons, Inc., 111 River Street, Hoboken, NJ 07030, (201) 748‐6011, fax (201) 748‐6008, or
online at http://www.wiley.com/go/permissions.
Limit of Liability/Disclaimer of Warranty: While the publisher and author have used their best efforts
in preparing this book, they make no representations or warranties with respect to the accuracy
or completeness of the contents of this book and specifically disclaim any implied warranties of
merchantability or fitness for a particular purpose. No warranty may be created or extended by sales
representatives or written sales materials. The advice and strategies contained herein may not be
suitable for your situation. You should consult with a professional where appropriate. Neither the
publisher nor author shall be liable for any loss of profit or any other commercial damages, including
but not limited to special, incidental, co; nsequential, or other damages.
For general information on our other products and services or for technical support, please contact
our Customer Care Department within the United States at (800) 762‐2974, outside the United States
at (317) 572‐3993 or fax (317) 572‐4002.
Wiley also publishes its books in a variety of electronic formats. Some content that appears in print
may not be available in electronic formats. For more information about Wiley products, visit our
web site at www.wiley.com.
Library of Congress Catalog Number: 58-9935
ISBN: 9781119165149
Set in 10/12pt Times by SPi Global, Pondicherry, India
Printed in the United States of America
10

9

8

7

6

5

4

3

2

1

Editorial Board
Kurt Binder, Condensed Matter Theory Group, Institut Für Physik, Johannes Gutenberg‐
Universität, Mainz, Germany
William T. Coffey, Department of Electronic and Electrical Engineering, Printing
House, Trinity College, Dublin, Ireland
Karl F. Freed, Department of Chemistry, James Franck Institute, University of Chicago,
Chicago, IL, USA
Daan Frenkel, Department of Chemistry, Trinity College, University of Cambridge,
Cambridge, UK
Pierre Gaspard, Center for Nonlinear Phenomena and Complex Systems, Université
Libre de Bruxelles, Brussels, Belgium
Martin Gruebele, Departments of Physics and Chemistry, Center for Biophysics and
Computational Biology, University of Illinois at Urbana‐Champaign, Urbana, IL, USA
Gerhard Hummer, Theoretical Biophysics Section, NIDDK‐National Institutes of
Health, Bethesda, MD, USA
Ronnie Kosloff, Department of Physical Chemistry, Institute of Chemistry and
Fritz Haber Center for Molecular Dynamics, The Hebrew University of Jerusalem,
Jerusalem, Israel
Ka Yee Lee, Department of Chemistry, James Franck Institute, University of Chicago,
Chicago, IL, USA
Todd J. Martinez, Department of Chemistry, Photon Science, Stanford University,
Stanford, CA, USA
Shaul Mukamel, Department of Chemistry, School of Physical Sciences, University of
California, Irvine, CA, USA
Jose N. Onuchic, Department of Physics, Center for Theoretical Biological Physics, Rice
University, Houston, TX, USA
Stephen Quake, Department of Bioengineering, Stanford University, Palo Alto, CA, USA
Mark Ratner, Department of Chemistry, Northwestern University, Evanston, IL, USA
David Reichman, Department of Chemistry, Columbia University, New York City, NY,
USA
George Schatz, Department of Chemistry, Northwestern University, Evanston, IL, USA
Steven J. Sibener, Department of Chemistry, James Franck Institute, University of
Chicago, Chicago, IL, USA

v

vi	Editorial Board	
Andrei Tokmakoff, Department of Chemistry, James Franck Institute, University
of Chicago, Chicago, IL, USA
Donald G. Truhlar, Department of Chemistry, University of Minnesota, Minneapolis,
MN, USA
John C. Tully, Department of Chemistry, Yale University, New Haven, CT, USA

Contents
Contributors List
Preface to the Series

ix
xi

Thermodynamic Perturbation Theory
for Associating Molecules

1

Bennett D. Marshall and Walter G. Chapman
Path Integrals and Effective Potentials in the Study
of Monatomic Fluids at Equilibrium

49

Luis M. Sesé
Spontaneous Symmetry Breaking in Matter Induced
by Degeneracies and Pseudodegeneracies

159

Isaac B. Bersuker
Mean Field Electrostatics Beyond the Point Charge
Description

209

Derek Frydel
First‐Passage Processes in Cellular Biology

261

Srividya Iyer‐Biswas and Anton Zilman
Theoretical Modeling of Vibrational Spectra and
Proton Tunneling in Hydrogen‐Bonded Systems

307

Marek Janusz Wójcik
Index

343

vii

Contributors List
Bennett D. Marshall, ExxonMobil Research and Engineering, Spring, TX, USA
Walter G. Chapman, Department of Chemical and Biomolecular Engineering, Rice
University, Houston, TX, USA
Luis M. Sesé, Departamento de Ciencias y Técnicas Fisicoquímicas, Universidad
Nacional de Educación a Distancia, Madrid, Spain
Isaac B. Bersuker, Institute for Theoretical Chemistry, Department of Chemistry,
University of Texas at Austin, Austin, TX, USA
Derek Frydel, Institute for Advanced Study, Shenzhen University, Shenzhen, China;
School of Chemistry and Chemical Engineering, Shanghai Jiao Tong University,
Shanghai, China; Laboratoire de Physico‐Chime Thèorique, ESPCI, CNRS Gulliver,
Paris, France
Srividya Iyer‐Biswas, Department of Physics and Astronomy, Purdue University,
West Lafayette, IN 47907, USA
Anton Zilman, Department of Physics and Institute for Biomaterials and Biomedical
Engineering, University of Toronto, Toronto, Ontario, Canada
Marek Janusz Wójcik, Faculty of Chemistry, Jagiellonian University, Kraków, Poland

ix

Preface to the Series
Advances in science often involve initial development of individual specialized
fields of study within traditional disciplines followed by broadening and overlap,
or even merging, of those specialized fields, leading to a blurring of the lines
between traditional disciplines. The pace of that blurring has accelerated in the
past few decades, and much of the important and exciting research carried out
today seeks to synthesize elements from different fields of knowledge. Examples
of such research areas include biophysics and studies of nanostructured materials.
As the study of the forces that govern the structure and dynamics of molecular
systems, chemical physics encompasses these and many other emerging research
directions. Unfortunately, the flood of scientific literature has been accompanied
by losses in the shared vocabulary and approaches of the traditional disciplines,
and there is much pressure from scientific journals to be ever more concise in the
descriptions of studies, to the point that much valuable experience, if recorded at
all, is hidden in supplements and dissipated with time. These trends in science
and publishing make this series, Advances in Chemical Physics, a much needed
resource.
The Advances in Chemical Physics is devoted to helping the reader obtain
general information about a wide variety of topics in chemical physics: a field that
we interpret very broadly. Our intent is to have experts present comprehensive
analyses of subjects of interest and to encourage the expression of individual
points of view. We hope that this approach to the presentation of an overview of a
subject will both stimulate new research and serve as a personalized learning text
for beginners in a field.
Stuart A. Rice
Aaron R. Dinner

xi

Thermodynamic Perturbation Theory
for Associating Molecules
BENNETT D. MARSHALL1 and WALTER G. CHAPMAN2
ExxonMobil Research and Engineering, Spring, TX, USA
Department of Chemical and Biomolecular Engineering, Rice University,
Houston, TX, USA
1

2

Contents
I.
II.
III.
IV.

Introduction
A Brief Introduction to Cluster Expansions
Single Association Site: Bond Renormalization
Single Association Site: Two‐Density Approach
A. The Monovalent Case
B. The Divalent Case
V. Multiple Association Sites: Multi‐Density Approach
VI. The Two‐Site AB Case
A. Steric Hindrance in Chain Formation
B. Ring Formation
C. Bond Cooperativity
VII. Spherically Symmetric and Directional Association Sites
VIII. Density Functional Theory
IX. Concluding Remarks
Acknowledgments
References

I.

INTRODUCTION

Since the time of van der Waals, scientists have sought to describe the m
­ acroscopic
behavior of fluids in terms of the microscopic interactions of the constituent
­molecules. By the early 1980s, accurate theories based on statistical mechanics had
primarily been developed for near-spherical molecules. Successes of the 1960s and
1970s particularly by Chandler, Weeks, and Andersen [1] and by Barker and Henderson
[2] produced perturbation theories for the properties of Lennard‐Jones (LJ) fluids.
Site–site theories such as reference interaction site model (RISM) [3] were developed,
Advances in Chemical Physics, Volume 160, First Edition. Stuart A. Rice and Aaron R. Dinner.
© 2016 John Wiley & Sons, Inc. Published 2016 by John Wiley & Sons, Inc.

1

2

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

in part, to provide reference fluid structure to extend these perturbation theories to
­polyatomic ­molecules. However, for certain classes of fluids, the accurate description
of the fluid phase in terms of the microscopic interactions has proven much more challenging. Hydrogen bonding interactions are strong, short‐ranged, highly directional
interactions that lie somewhere between a dipole/dipole attraction and a covalent
bond. The short range and directionality of hydrogen bonds result in the phenomena
of bond saturation, giving a limited valence of the hydrogen bonding attractions.
The same properties of the hydrogen bond, which complicate the theoretical
description of these fluids, also give rise to a number of macroscopic physical
properties that are unique to fluids that exhibit hydrogen bonding. Hydrogen
bonding is responsible for the remarkable properties of water [4], folding of
­proteins [5] and is commonly exploited in the self‐assembly [6] of advanced
materials. More recently patchy colloids, a new class of materials that shares
many qualities with hydrogen bonding fluids, have been developed. Patchy
­colloids are colloids with some number of attractive surface patches giving rise to
association like anisotropic inter‐colloid potentials [7]. For the purposes of this
review, patchy colloids and hydrogen bonding fluids are treated on equal footing
and will simply be termed “associating fluids.”
The first models used to describe hydrogen bonding fluids were developed
using a chemical approach, where each associated cluster is treated as a distinct
species created from the reaction of monomers and smaller associated clusters
[8, 9]. The “reactions” are governed by equilibrium constants that must be obtained
empirically. This type of approach has been incorporated into various equations of
state including a van‐der‐Waals‐type equation of state [10], the ­perturbed anisotropic chain theory equation of state (APACT) [11], and the Sanchez–Lacombe
[12] equation of state.
Alternatively, lattice theories may be employed to model hydrogen bonding
fluids. These approaches generally follow the method of Veytsman [13] who
showed how the free energy contribution due to hydrogen bonding could be
­calculated in the mean field by enumerating the number of hydrogen bonding
states on a lattice. Veytsman’s approach was incorporated into the Sanchez–
Lacombe equation of state by Panayiotou and Sanchez [14] who factored the
­partition function into a hydrogen bonding contribution and a non‐hydrogen
bonding contribution. The lattice approach has also been applied to hydrogen
bond cooperativity [15] and intramolecular [16] hydrogen bonds.
Both the chemical and lattice theory approaches to hydrogen bonding yield
semi‐empirical equations of state, which are useful for several hydrogen bonding
systems [8]. The drawback of these approaches is a result of their simplistic development. As discussed earlier, it is desired to describe the macroscopic behavior of
fluids through knowledge of the microscopic intermolecular interactions and
­distributions. This cannot be accomplished using a lattice or chemical theory.
To accomplish this goal, we must incorporate molecular details of the associating
fluid from the outset.

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

3

The starting place for any molecular theory of association is the definition of
the pair potential energy ϕ(12) between molecules (or colloids). Molecules are
treated as rigid bodies with no internal degrees of freedom. In total, six degrees of
freedom describe any single molecule: three translational coordinates represented

by the vector r1 and three orientation angles represented by Ω1. These six degrees

of freedom are represented as 1 r1 , 1 . It is assumed that the intermolecular
potential can be separated as
12

12

ref

as

12

(1)

where ϕas(12) contains the association portion of the potential and ϕref(12) is the
reference system potential, which contains all other contributions of the pair
potential including a harsh short‐ranged repulsive contribution.
Considering molecules (or colloids) that have a set of association sites
{A, B, C , , Z}, where association sites are represented by capital letters, the
association potential is decomposed into individual site–site contributions
as

12

AB
A

12

(2)

B

The potential ϕAB(12) represents the association interaction between site A on molecule 1 and site B on molecule 2. One of the challenges in developing theoretical
models for associating fluids stems from the short‐ranged and directional nature of
the association potential ϕAB, which results in the phenomena of bond saturation. For
instance, considering molecules which consist of a hard spherical core of diameter d

ref

12

HS

r12

r d
0 r d

(3)

and a single association site A (see Fig. 1), bond saturation arises as follows.
When spheres 1 and 2 are positioned and oriented correctly such that an association bond is formed, the hard cores of these two spheres may, depending on the

3

Figure 1. Illustration of bond saturation for hard spheres with a single monovalent ­association site.
(See insert for color representation of the figure.)

4

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

size and range of the association site, prevent sphere 3 from approaching and
­sharing in the association interaction. That is, if as (12) 0 and as (13) 0, then
, meaning that each association site is singly bondable (has a valence
HS (r23 )
of 1). In hydrogen bonding it is usually the case that each association bond site is
singly bondable, although there are exceptions. For the case of patchy colloids,
the patch size can be controlled to yield a defined valence controlling the type of
self‐assembled structures that form.
Conical square well (CSW) association sites are commonly used as a primitive
model for the association potential ϕAB. First introduced by Bol [17] and later
reintroduced by Chapman et al. [18, 19], CSWs consider association as a square
well interaction which depends on the position and orientation of each molecule.
Kern and Frenkel [20] later realized that this potential could describe the interaction between patchy colloids. For these CSWs the association potential is given by
AB

12

O AB 12
f AB 12

AB

O AB 12

1, r12
0
exp

rc ;

; B2
otherwise

AB

/k bT

A1

c

1 O AB 12

c

(4)
f AB O AB 12

where rc is the maximum separation between two colloids for which association

can occur, θA1 is the angle between r12 and the orientation vector passing through
the center of the patch on colloid 1, and θc is the critical angle beyond which
­association cannot occur. Equation (4) states that if the spheres are close enough
r12 rc , and both are oriented correctly A1 c and B 2 c, then an association
bond is formed and the energy of the system is decreased by εAB. Figure 2 gives an
illustration of two single‐site spheres interacting with this potential. The size of
the patch is dictated by the critical angle θc that defines the solid angle to be
2 (1 cos c ). The patch size determines the maximum number of other spheres
to which the patch can bond. Specifically considering a hard sphere reference
fluid with association occurring at hard sphere contact rc d , it is possible for a

θc
θA1

r12
θA2

Figure 2. Association parameters for conical association sites. (See insert for color representation of the figure.)

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

5

patch to associate at most once for 0 c 30, twice for 30 c 35.3, thrice for
35.3 c 45, and four times for 45 c 58.3 [21]. The advantage of the CSW
model is that it separates the radial and angular dependence of the potential and
allows for analytic calculations in the model while allowing for quick calculation
of association in a simulation since only two dot products are needed to determine
that the molecular orientation criteria is satisfied for association.
In the following sections we review some of the existing theories to model
associating fluids with potentials of the form of Eqs. (1)–(2). We focus mainly
on the multi‐density formalism of Wertheim [22, 23], which has been widely
applied across academia and industry. In Sections III and IV.A, only association
sites that are singly bondable are considered and steric hindrance between association sites is neglected. Extensions of Wertheim’s multi‐density approach for
the divalent case is described in Section IV.B. Section V addresses the case of
multiple association sites on a molecule within Wertheim’s multi‐density formalism. Section VI extends the theory to the case of a small angle between two
association sites, so that the sites cannot be assumed to be independent, and for
the case of cooperative hydrogen bonding. Section VII extends the theory to
account for association interactions between molecules with spherically symmetric and directional association sites (e.g., ion–water solvation). A brief
description of applying the density functional theory (DFT) approach for associating molecules is presented in Section VIII. Finally, Section IX gives concluding remarks. Prior to exploring the theory, a brief introduction to cluster
expansions is provided in Section II.
II. A BRIEF INTRODUCTION TO CLUSTER EXPANSIONS
In this section we give a very brief overview of cluster expansions. For a more
detailed introduction the reader is referred to the original work of Morita and
Hiroike [24] and also to the reviews by Stell [25] and Andersen [26]. Cluster
expansions were first introduced by Mayer [27] as a means to describe the
­structure and thermodynamics of non‐ideal gases. In cluster expansions Mayer f
functions are introduced:
f 12

exp

12
k BT

1

(5)

The replacement exp ( (12) / kBT ) f (12) 1 in the grand partition function
and the application of the lemmas developed by Morita and Hiroike [24] allows for
the pair correlation function g(12) and Helmholtz free energy A to be written as an
infinite series in density where each contribution is an integral represented
pictorially by a graph. A graph is a collection of black circles and white
­
­circles with bonds connecting some of these circles. The bonds are represented by

6

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

(a)
= ʃρ(2)ρ(3)ρ(4)f (12)f (23)d(2)d(3)d(4)
1

(b)

=

1
ʃρ(1)ρ(2)ρ(3)ρ(4)f (12)f (23)f (34)d(1)d(2)d(3)d(4)
2

=

1
ʃρ(1)ρ(2)ρ(3)ρ(4)f (12)f (23)f (34)f (14)d(1)d(2)d(3)d(4)
8

(c)

Figure 3.

Examples of integral representations of graphs. Arrows point towards articulation

circles.

two‐­molecule functions such as Mayer functions f(12) and the black circles are
called “field points” represented by single‐molecule functions such as fugacity
z(1) or density ρ(1) integrated over the coordinates (1). The white circles are called
“root points” and are not associated with a single‐molecule function, and the coordinates of a root point are not integrated. Root points are given labels 1, 2, 3,….
The value of the diagram is then obtained by integrating over all coordinates
­associated with field points and multiplying this integral by the inverse of the
symmetry number S of the graph. Figure 3 gives several examples. The volume

element d(1) is given by d (1) dr1 d 1, representing the differential position and
orientation of the molecule.
Before giving graphical representations of the pair correlation function g(12)
and the Helmholtz free energy A, a few definitions must be given as follows:
1. A graph is connected if there is at least one path between any two points.
Graph a in Fig. 3 is disconnected, and graphs b and c are both connected.
2. An articulation circle is a circle in a connected graph whose removal makes
the graph disconnected, where at least one part contains no root point and at
least one field point Arrows in Fig. 3 point to articulation circles.
3. An irreducible graph has no articulation circles. Graph c in Fig. 3 is an
example of an irreducible graph.
Using these definitions, the pair correlation function and Helmholtz free energy
are given as
g 12

sum of topologically different irreducible graphs that have two root
points labelled 1 and 2, any number of field points , and at most
one f bond between each pair of points
(6)

7

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

and
A
k BT

1 ln

1

3

1 d 1

c

o

where Λ is the de Broglie wavelength, ρ(1) is the density where
and c(o) is the graph sum given by

c

o

(7)

r

1 d ,

sum of topologically different irreducible graphs that have noo root
points any number of field points , and at most one f bond between
each pair of points
(8)

Equations (6)–(8) are rigorous and exact mathematical statements. Unfortunately,
the exact evaluation of these infinite sums cannot be performed and numerous
approximations must be made to obtain any usable result. In these approximations
only some subset of the original graph sum is evaluated.
Performing these partial summations in hydrogen bonding fluids is complicated by both the strength of the association interaction and the limited valence of
the interaction. Hydrogen bond strengths can be many times that of typical van
der Waals forces giving Mayer functions which are very large. If the entire cluster
series were evaluated for g(12) and c(o) many of these large terms would cancel;
however, when performing partial summations, care must be taken to eliminate
divergences if meaningful results are to be obtained. Similarly, in most hydrogen
bonding fluids, each hydrogen bonding group is singly bondable. Hence, any
theory for hydrogen bonding fluids must account for the limited valence of the
attractions. Again, if the full cluster series were evaluated for g(12) this condition
would be naturally accounted for; however, when performing partial summations
care must be taken to ensure this single bonding condition holds. There have been
three general methods to handle these strong association interactions using cluster
expansions. The first was the pioneering work of Andersen [28, 29] who ­developed
a cluster expansion for associating fluids in which the divergence was tamed by
the introduction of renormalized bonds, the second is the approach of Chandler
and Pratt [30] who used physical clusters to represent associated molecules, and
the third is the method of Wertheim [22, 23, 31–33] who used multiple densities.
Both Andersen and Wertheim took the approach of incorporating the effects of
steric hindrance early in the theoretical development in the form of mathematical
­clusters. In what follows, for brevity, we restrict our attention to the approaches of
Andersen and Wertheim, however, when possible we draw parallels between these
approaches and that of Chandler and Pratt.

8

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

III. SINGLE ASSOCIATION SITE: BOND RENORMALIZATION
Before discussing the more general case of associating fluids with multiple
­association sites, we will discuss the simpler case of molecules with a single association site A. For a single association site, the Mayer function is decomposed as
f 12

fref 12

FAA 12

(9)

where
FAA 12

eref 12 f AA 12

eref 12

exp

ref

f AA 12

exp

AA

12

kB T
12

kBT

1 fref 12

(10)

1

In Eq. (10) the fAA(12) accounts for the anisotropic/short‐ranged attraction of the
association interaction and the function eref(12) prevents the overlap of the cores
of the molecules. It is the functions eref(12) that give rise to the single bonding
condition. Now inserting Eq. (9) into Eq. (6) and simplifying

g 12

sum of topologically different irreducible graphs that have two
root points labelled 1 and 2, any number of field points , fref and
FAA bonds, and at most one bond between each pair of points
(11)

Andersen [28, 29] defines a renormalized association Mayer function FAA (12)
as the sum of the graphs in Eq. (11) which are most important in the determination
of g(12). Since the Mayer functions FAA may take on very large numerical values
in the bonding region, the most important graphs in the calculation of g(12) are
the ones whose root points are connected by an FAA bond. Hence, it is natural to
define FAA as
FAA 12

sum of graphs in 11 whose root points are connected
by an FAA bond

(12)

Andersen assumes that the intermolecular potential was such that the association
site was singly bondable. This single bonding condition was exploited in the

9

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

c­ luster expansion by use of the cancelation theorem as described by Andersen,
who was able to sum the diagrams in Eq. (12) as
FAA 12

FAA 12 Yp 12

1 2

1 4

AA

2

AA

where the term ΔCD is given by (where for a homogeneous fluid
CD

and

d
Yp 12

(13)

2
AA

Yp 12 FCD 12 d 2

(1) d

)

1

(14)

is the total number of orientations. The function Yp(12) is given by
sum of graphs in 11 which have no FAA bond attached
to either root, and no bond between the roots

(15)

It is easily shown that FAA (12) is bounded as follows:
FAA 12 d 2

0

1

(16)

Equation (16) shows that the renormalized association bond remains finite even
when the association potential ϕAA takes on infinitely large negative values. Using
this renormalized bond the average number of hydrogen bonds per molecule is
calculated as follows:
FAA 12 d 2

N HB

(17)

Comparing Eqs. (16) and (17) it is easy to see
0

N HB

1

(18)

Equation (18) demonstrates that the single bonding condition is satisfied and that
the method of Andersen was successful. Unfortunately, the function Yp(12) must
be obtained through the solution of a series of integral equations using approximate closures.
To the author’s knowledge, this approach has never been applied for numerical
calculations of the structure or thermodynamics of one‐site‐associating fluids.
Here we will show how a single simple approximation allows for the calculation

10

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

of NHB. To approximate Yp(12) we note that this function can be decomposed into
contributions from graphs that contain k association bonds FAA
Yp 12

Yp

k

12

(19)

k 0

k

The terms Yp give the contribution to Yp from graphs that contain k association
bonds. The simplest possible case is to keep only the first contribution k = 0 and
k
disregard all Yp for k > 0. For this simple case
Yp 12

yref 12

(20)

where yref is the cavity correlation function of the reference fluid, meaning association is treated as a perturbation to the reference fluid. This approximation is not
necessarily intuitive since the structure of a fluid is expected to be strongly
affected by association. Combining these results, the monomer fraction (fraction
of molecules that do not have an association bond) can be written as
Xo

1 N NB

1

1 4
2

AA

(21)

AA

where ΔAA is now given by
yref (12) FAA (12) d (2). Equation (21) gives
AA
a very simple relationship for the monomer fraction. This same equation was later
derived by Chandler and Pratt [30] and Wertheim [22] using very different cluster
expansions. Equation (21) has been shown to be highly accurate in comparison
to simulation data [19, 34, 35]. Now we will introduce Wertheim’s two‐density
formalism for one‐site‐associating fluids.
IV. SINGLE ASSOCIATION SITE: TWO‐DENSITY APPROACH
In the previous section it was shown that Andersen’s formalism can be applied to
derive a highly accurate and simple relationship for the monomer fraction. In
order to obtain this result the renormalized association Mayer functions FAA were
employed. The applicability of Andersen’s approach to more complex systems
(mixtures, multiple bonds per association site, etc.) is limited by the fact that for
each case the renormalized Mayer functions must be obtained by solving a rather
complex combinatorial problem. A more natural formalism for describing association interactions in one‐site‐associating fluids is the two‐density formalism of
Wertheim [22, 31].
Instead of using the density expansion of the pair correlation function g(12) or
Helmholtz free energy A, Wertheim uses the fugacity expansion of ln Ξ, where Ξ

11

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

is the grand partition function, as the starting point. Building on the ideas of
Lockett [36], Wertheim then regroups the expansion such that individual graphs
are composed of s‐mer graphs. An s‐mer represents a cluster of points that are
connected by paths of FAA bonds; each pair of points in an s‐mer, which are not
directly connected by a FAA bond, receives an eref(12) bond. This regrouping serves
to include the geometry of association with the eref(12) bonds enforcing the limited valence of the association interaction. In the s‐mer representation, graphs that
include unphysical core overlap are identically zero. That is, if the association site
is singly bondable all graphs composed of s‐mers of size s > 2 immediately vanish
due to steric hindrance. This is not the case in Andersen’s approach where these
unphysical contributions are allowed in individual graphs, with the single bonding
condition being exploited with the cancelation theorem.
This regrouping of the fugacity expansion allows for the easy incorporation of
steric effects. Now, unlike Andersen who tamed the arbitrarily large FAA bonds
through the introduction of a renormalized FAA , Wertheim uses the idea of multiple
densities, splitting the total density of the fluid as
1

o

1

b

1

(22)

where ρo(1) is the density of monomers (molecules not bonded) and ρb(1) is the
density of molecules that are bonded. The density ρo(1) is composed of all graphs
in ρ(1) which do not have an incident FAA bond, and ρb(1) contains all graphs
which have one or more incident FAA bonds. Performing a topological reduction
from fugacity graphs to graphs which contain ρo(1) and ρ(1) field points, allowed
Wertheim to arrive at the following exact free energy
A
kBT

1 ln

o

1

3
o

1 d 1

c

o

(23)

where for this case the graph sum c(o) is given as follows:

c

o

sum of all irreducible graphs consisting of monomer points carrrying
factors of , s -mer graphs with s 2 and every point carrying a factor
of o , and fref -bonds between some sets of points in distinct s -mers.
(24)

The first few graphs in the infinite series for c(o) are given in Fig. 4. In Fig. 4
crossed lines
represent FAA bonds, dashed lines represent eref bonds, and
solid lines represent fref bonds. All points with one or more incident FAA bonds
carry a factor ρo(1), and each point with no incident FAA bonds carries a factor
ρ(1). All graphs without any FAA bonds (graphs a, c, g, h, and i in Fig. 4) represent

12

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

(a)

(b)

(c)

+

c(o) =

(d)

+

(e)

+

(f)

(g)

(h)

+

+

(i)

+

(j)

(k)

+

(l)
+

+

(m)

(n)

(o)
+

+

+ •••

Figure 4. Graphical representation of Eq. (24) where crossed lines
FAA bonds, dashed lines represent eref bonds, and solid lines represent fref bonds.

represent

o

the ­reference system contribution cref . Any point that has two incident FAA bonds
(graphs e and f in fig. 4 are s = 3‐mers) represents a molecule with an association
site which is bonded to two other molecules.
A. The Monovalent Case
If ϕ(12) is chosen such that the single bonding condition holds, then all s‐mer
graphs with s > 2 vanish (e.g., graphs e and f in Fig. 4 are zero) and Eq. (24) can
be summed exactly to yield the following:
c

o

o

cref

1
2

o

1 f AA 12 goo 12

o

2 d 1 d 2

(25)

Note that Eq. (25) contains monomer densities since only monomers can associate. The use of monomer densities bounds the association term. The quantity
goo(12) is the monomer/monomer pair correlation function which can be ordered
by graphs that contain k FAA bonds:
k

goo 12

goo 12
k 0

(26)

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

13

Similar to the approximation made for Yp in the formalism of Andersen (Eq. 19),
only the lowest‐order contribution is retained, and all contributions with k > 0 are
neglected. This is the single chain approximation, that yields the following:
goo 12

gref 12

yref 12 eref 12

(27)

Equation (27) forms the basis of Wertheim’s TPT, which assumes the monomer–
monomer correlation function is the same as that of the reference fluid. This
amounts to neglecting all graphs in Eq. (24) which contain more than a single FAA
bond (e.g., neglecting graphs n and o in Fig. 4). Although this is the same approximation that we introduced in Eq. (20) for Andersen’s theory, the approximation is
more intuitive in terms of the monomer–monomer distribution. Considering a
dense fluid of hard spheres associating at contact, the packing ­fraction of the fluid
does not change with extent of association. Therefore, we might expect that the
monomer–monomer distribution function would remain relatively unchanged
with association. For association near hard sphere contact (or sigma for LJ
­molecules), molecular simulation results show this to be a r­easonably accurate
approximation [18, 19, 35, 37–39].
To obtain an equation for ρo(1), Eq. (23) is minimized:
A / kB T
o 1

1
o

1

1

f AA 12 gref 12

o

2 d 2

0

(28)

The operator δ/δρo(1) represents the functional derivative. Combining Eqs. (23),
(25), and (28) the free energy is simplified as
A Aref
kB T

1 ln X o 1

Xo 1
2

1
d 1
2

(29)

where Aref is the Helmholtz free energy of the reference system and
X o (1)
(1) is the fraction of monomers. Now, assuming a homogeneous
o (1) /
fluid
(1) and solving Eq. (28), the monomer fraction Eq. (21) is obtained.
As can be seen, under the single bonding condition when treated as a perturbation theory, both Andersen’s and Wertheim’s approaches give the same result for
homogeneous fluids. Indeed, Eq. (13) can be rewritten in terms of monomer
fractions
FAA 12

FAA 12 yref 12 X 02

f AA 12 gref 12 X 02

(30)

14

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

and for a homogeneous fluid the renormalized Mayer functions can be used to
represent c(o).
c

o

o

cref

1V
2

2

FAA 12 d 2

(31)

Note that the monomer fraction provides the scaling that keeps the perturbation
bounded even for large association energies. While Andersen’s and Wertheim’s
approaches produce identical results for singly bondable sites in the single chain
approximation (perturbation theory), the two‐density approach of Wertheim is
much more versatile than the approach of Andersen. For instance, for the case that
the association site can bond a maximum of n times, there is a clear path forward
in the development of a perturbation theory using Wertheim’s approach (keep all
s‐mer graphs with s n). Attempting to apply Andersen’s formalism to this case
would be hopelessly complex. Also, Eqs. (28) and (29) are generally valid for
inhomogeneous fluids where the density and monomer fraction vary with position
and orientation. In fact, DFTs based on Wertheim’s TPT have proven to be very
accurate in the description of inhomogeneous one‐site‐associating fluids [40, 41].
It seems unlikely the approach of Andersen could be utilized to derive the inhomogeneous form of the theory.
The accuracy of the theory for hard spheres and LJ spheres with a single association site is remarkable in comparison with molecular simulation results for the
extent of association, fluid pressure, and internal energy to high association energy
[18, 35, 37]. In the limit of infinite association energy, the theory accurately
­predicts the equation of state for a fluid of diatomic hard spheres or diatomic
LJ molecules [18, 33, 35, 37]. Interestingly, in the limit of infinite association
energy, the residual free energy in the theory predicts a correction to the ideal gas
term to convert from an ideal gas of spheres to an ideal gas of diatomics. For LJ
diatomics, the theory accurately predicts the change in potential energy from a
fluid of independent LJ spheres to a fluid of LJ diatomics. Accurately predicting
the equation of state of the species produced by association is necessary to
­accurately model the association equilibrium.
B. The Divalent Case
One of the main assumptions in the development of Wertheim’s first‐order
­thermodynamic perturbation theory (TPT1) is that association sites are singly
bondable. Indeed, the entire multi‐density formalism of Wertheim is constructed
to enforce this condition. For the case of hydrogen bonding, the assumption of
singly bondable sites is justified; however for patchy colloids (see Section I for a
background on patchy colloids), it has been shown experimentally [42, 43] that
the number of bonds per patch (association site) is dependent on the patch size.
It has been 30 years since Wertheim first published his two‐density cluster

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

15

(a)

(b)

(c)

Figure 5. Associated clusters for patchy colloids with a single double bondable patch: (a) dimers,
(b) chains with double bonded sites, and (c) triatomic rings of double‐bonded sites. (See insert for color
representation of the figure.)

e­ xpansion for associating fluids, and only very recently have researchers applied
his approach (or a similar approach developed for associating fluids with spherically symmetric association potentials [44]) to the case that association sites are
divalent [21, 45, 46].
Application of TPT to divalent association sites is complicated by the fact that
three‐body terms must be included to account for blocking effects caused between
two molecules attempting to bond to an association site on a third molecule. For a
pure component fluid of associating spheres with a single divalent association site
the dominant types of associated clusters are chains and triatomic rings of doubly
bonded sites as shown in Fig. 5. The application of TPT to this divalent case is an
excellent teaching example of how to extend TPT beyond first‐order (­monovalent
sites). For clarity we consider the specific case of a homogeneous fluid of patchy
hard spheres (PHS) whose potential model is defined with a hard sphere reference
(Eq. 3) and a single conical square well association site (Eq. 4).
To begin we first separate Δc(o) into contributions for chain and ring formation
as follows:
c

o

o
cˆ1

o

cchain

o

cring

(32)

o

The contribution cchain is further decomposed into contributions from chains of n
bonds and n + 1 colloids as follows:
o

cnchain

cchain
n 1

(33)

16

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

For instance, graph e in Fig. 4 belongs to the contribution c2chain and graph f
o
belongs to cring. Each of these contributions consists of an infinite series of graphs
with a single associated cluster interacting with the reference fluid. These series
can be summed as follows:
cnchain
V

1
2

d k 1

n

n 1 n
o
AA

gHS 1

f

n 1

AA

k, k 1

(34)

k 1

and
cring

1
6

V

3 3
o AA

gHS 123

f

AA

12

AA

23

AA

13

d 2 d 3

(35)

Here
4 and the functions gHS(1 … k) are the k body correlation functions of
the hard sphere reference system. Since little is known about the correlation
­functions gHS(1 … k) for k > 3, we must approximate the higher order gHS(1 … k) in
superposition. For the current case, a particularly convenient approximation for
the chain contributions will be the following:
k 1

gHS 1

k

k 2

gHS rj , j

eHS ri ,i

1

j 1

(36)

2

i 1

The superposition given by Eq. (36) prevents overlap between nearest and next
nearest neighbors in the chain and should be most accurate at low densities. We
note that the probability that an isolated associated chain of n + 1 colloids has a
configuration (1 n 1) is given by the following equation:
n 1

n
n
chain

P

AA

1

k, k 1 eHS rk , k

eHS ri ,i

1

k 1

n 1

2

i 1

(37)

n

Z chain

The probability in Eq. (37) accounts for steric interactions between nearest and
n
next nearest neighbors in the chain, and the term Z chain is the chain partition
­function given by the following equation:
n

n 1

Z chain

n

eHS ri ,i
i 1

2

AA

k, k 1 eHS rk , k

1

d k 1

(38)

k 1

Combining Eqs. (34) and (36)–(38) we obtain the following:
cnchain
V

1
2

n

n 1 n
o
AA

f

Z ch

n

yHS rj , j

n
j 1

(39)

1
Pchain

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

17

The cavity correlation function yHS (rj , j 1 ) gHS (rj , j 1 ) / eHS (rj , j 1 ). The brackets
in Eq. (39) represent an average over the distribution function Eq. (37). To an
excellent approximation this average can be evaluated as a product of individual
­averages over the bonding range
n

yHS rk , k

yHS r

1

j 1

n
br

(40)

Pchain

where
rc

yHS r r 2 dr

4
yHS r

d

(41)

rc

br

b

2

4

r dr
d

The constant νb is the volume of a spherical shell defined by the denominator of
the second term in Eq. (41) and ξ is defined by the numerator. As has been shown
n
[21], integrals of similar form to Z chain can be very accurately factored as
1

n

Z chain

Z chain

n

n 1
chain

(42)

where
2

Z chain
chain

Z

1
chain

1
2

2

AA

12

AA

23 eHS r12 eHS r23 eHS r13

d 2 d 3

b

(43)
The accuracy of the factorization in Eq. (42) stems from the fact that double
­bonding of a patch is dominated by two‐ and three‐body effects. When chain 0,
­multiple bonding of an association site (patch) is impossible; while for the case
1, there is no steric hindrance between two PHS bonding to the same
chain
­association site on a third. Indeed, the geometric integral Φchain encodes the effect
of steric hindrance for doubly bonded sites. Combining the previous results, we
obtain the following:
cnchain
V

1
2

n 1
o

f AA

n

n 1
chain

(44)

18

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

The constant κ represents the probability that two colloids are in mutual bonding
orientations and is given by the following equation:
2

1 cos

c

(45)

4

Now using Eq. (44) to evaluate the infinite sum in Eq. (33) we obtain the following:
o

2
o AA

f

1
2

cchain
V

1 f AA

(46)

o

chain

When multiple bonding of a patch is impossible chain 0, and we recover the
TPT1 result Eq. (25) for the case goo = gHS.
Now we turn our attention to the ring contribution cring Eq. (35). For this case
we approximate the triplet correlation function using Kirkwood superposition.
Following a similar process to the one desribed above in the development of
Eq. (46), we obtain the following result:
o

cring

1
6 b

V

3

f

2

o AA

(47)

ring

In Eq. (47) Φring is given by
1
ring

AA

2

12

23

AA

AA

d 2 d 3

13 eHS r12 eHS r23 eHS r13

b

(48)
When multiple bonding of a site becomes impossible, ring 0, resulting in
o
cring 0. Now that Δc(o) has been completely specified the free energy is
­minimized with respect to ρo giving the following relation:
2
o AA

f

o

1 f AA

o

chain

1
2

3
o

chain

2

f AA
1 f AA

o

chain

1
2 b

f

o AA

3

2
ring

(49)
Equation (49) is simply conservation of mass. From Eq. (49) we identify the
­density of colloids bonded twice in rings 2ring , bonded once in a chain 1chain , and
bonded twice in a chain 2chain as follows:
ring
2

1
2 b

f

o AA

3

2
ring

(50)

19

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

chain
1

chain
2

1
2

3
o

2
o AA

f

1 f AA

(51)
o

chain

2

f AA
chain

1 f AA

(52)
o

chain

Using these density definitions the free energy can be simplified to
A AHS
VkBT

ln

chain
1

o
o

2

ring
2

3

(53)

Equation (53) completes the theory for molecules/colloids with a single doubly
bondable association site. To obtain the free energy, Eq. (49) is first evaluated for
ρo which allows the free energy to be calculated through Eq. (53).
All that remains to be done is to calculate the integrals Φchain (43) and Φring (48).
To evaluate these integrals we exploit the mean value theorem and employ Monte
Carlo integration [47] to obtain the following:

chain

The probability that if the positions of two colloids are geneerated
such that they are correctly positioned to associate with a thhird colloid,
that there is no core overlap between the two generatedd colloids
(54)

ring

The probability that if the positions and orientations of two colloids
are generated such that they are positioned and oriented coorrectly to
bond to a third colloid, that there is no core overlap betw
ween the
two generated colloids and that these two generated colloid
ds are
oriented and positioned such that they share an association boond
(55)

Equations (54) and (55) are easily evaluated on a computer; the calculations are
independent of temperature and density—they depend only on the potential
parameters rc and θc. Table I gives calculations for a critical radius of rc 1.1d.
Numerical results are given in Fig. 6, for theoretical predictions of the reduced
excess internal energy E * E AS / NkBT as well as the fraction of PHS that are
bonded twice in chains X 2chain and rings X 2ring. Results are plotted against reduced
association energy * AA / kBT at a packing fraction of
0.2. For comparison

20

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

Table I
Numerical Calculations of Integrals Φchain and Φring for a
Critical Radius rc 1.1d
θc
27°
28°
29°
30°
31°
32°
33°
34°
35°
40°
45°
50°
55°

Φchain

Φring

0
2.89 × 10−5
5.91 × 10−4
2.82 × 10−3
7.42 × 10−3
1.44 × 10−2
2.35 × 10−2
3.45 × 10−2
4.70 × 10−2
0.123
0.207
0.285
0.355

0
0
0
0
5.10 × 10−8
2.41 × 10−6
1.79 × 10−5
6.20 × 10−5
1.51 × 10−4
1.52 × 10−3
4.37 × 10−3
8.03 × 10−3
1.18 × 10−2

(a)

(b)

0

1

–2

0.8

TP

T1

–4
–6

E*

0.4

–8

0.2

–10
–12

X2ring

0.6

0

2

4

6
ε*

8

10

12

0

X2chain
0

2

4

6

8

10

12

ε*

Figure 6. Left: Reduced internal energy versus reduced association energy. Dashed curve
gives theory predictions assuming a monovalent association site and solid curve gives theory predictions when double bonding of a site is accounted for and symbols give Monte Carlo simulation [45]
results. Right: Fraction of spheres bonded twice in chains and spheres. Squares and circles give
respective Monte Carlo simulations [45], and curves are from divalent theory.

we include the simulation results of Marshall et al. [45] and predictions of TPT1.
As can be seen, TPT1 significantly under predicts the magnitude of E* due to the
fact that the possibility of two bonds per patch is not accounted for. The theory
derived here (solid curve) is in excellent agreement with the simulation data over
the full range of ε*. In addition to the internal energy, the theory also accurately
predicts the structure of the fluid. In agreement with simulation the theory shows
that triatomic rings dominate at strong association.

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

21

V. MULTIPLE ASSOCIATION SITES: MULTI‐DENSITY APPROACH
{A, B, C , , Z} will
Now the case of molecules with a set of association sites
be considered. The majority of hydrogen bonding molecules contain multiple
association sites: water, alcohols, proteins, hydrogen fluoride, etc. Theoretically,
this case is more difficult to model than the single‐site case due to the fact that
these molecules can form extended hydrogen bonded structures.
The two‐density approach of Wertheim allows the development of accurate and
simple theories for molecules with a single association site. To extend this idea to
the case of multiple association sites n ( ) 1, Wertheim again begins with the
fugacity expansion of ln Ξ which he regroups into the s‐mer representation.
Where, as for the one‐ site case, an s‐mer represents a cluster of s points
(­hyperpoints here) connected by association bonds fij. However, in contrast to the
two‐density case, all points in an s‐mer are not connected by eref bonds. Only
points with bond‐connected association sites within an s‐mer are connected by eref
bonds. Wertheim defines two association sites as bond connected if there is a
­continuous path of association sites and bonds between these two association
sites. Figure 7 demonstrates this for the case of a two‐site AB molecule. The wavy
lines represent fAB bonds and the dashed lines represent eref bonds, remembering
f AB (12)eref (12) FAB (12). All molecules that share fAB bonds are bond connected
receiving eref bonds. Molecules that do not share association bonds (e.g., molecules 1 and 3 and 3 and 5 in Fig. 7) can only be bond connected if an association
site is bonded more than once. This is the only way two association sites not
directly connected by a fAB bond can be connected by a continuous path of sites
and fAB bonds. For this reason, molecules 1 and 3 receive an eref bond and molecules 3 and 5 do not. This choice to only fill with eref bonds between bond connected sites greatly facilitates the formulation of approximation methods.
2

3
5
4

1

Figure 7. Representation of graph for two‐site‐associating fluids, where wavy lines represent
association bonds and dashed lines represent reference system e bonds. (See insert for color representation of the figure.)

22

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

In the two‐density formalism for one‐site‐associating molecules, separate
d­ ensities were assigned to molecules that were bonded and those that were not
bonded. For multiple association sites this choice would result in the loss of
­information on site–site‐level interactions. For this reason, Wertheim expresses
the total density as the sum over densities of individual bonding states of the
molecules
1

1

(56)

where ρα(1) is the density of molecules bonded at the set of sites α. For example,
ρAB(1) is the density of molecules with sites A and B bonded. To aid in the reduct­
ion to irreducible graphs Wertheim defines the density parameters:
1

1

(57)

Two important cases of Eq. (57) are o
. Using these density
o and
parameters, Wertheim transforms the theory from a fugacity expansion to an
expansion in σγ through the use of topological reduction, ultimately arriving at the
following exact free energy.
A
kBT

1 ln

o

1

3

Q 1 d 1

c

o

(58)

The graph sum in Eq. (58) is now defined as follows:

c

o

sum of all irreducible graphs consisting of s -mer graphs (including
monomer hyperpoints) and fR bonds. Points which are bonded
at a set of sites carry a factor
1
(59)

The term Q(1) is given by
Q 1

1

1 c 1

(60)

with
c 1

c

o

1

(61)

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

23

The densities are related to the cα(1) by the relation
1

o

1

c 1

(62)

P

where P( ) { } is the partition of α into subsets indexed by {γ}. For instance,
the density ρAB(1) is given by AB (1)
c A (1)cB (1)) .
o (1)(c AB (1)
The reference system Helmholtz free energy is given as follows:
Aref
kB T

1 ln

1

3

1 d 1

o

(63)

cref

(o)
The reference system cref
contains all of the graphs in Eq. (59) which are devoid
of association bonds. Subtracting Eq. (63) from (58) we obtain the following:

A Aref
kB T

o

1 ln

1

Q 1

1

1 d 1

c

o

(64)

(o)
The association graph sum c ( o ) c ( o ) cref
contains all the graphs in Eq. (59)
which contain association bonds. Figure 8 gives examples of graphs in the sum
Δc(o) for the two‐site AB case discussed in Fig. 7.
Equation (64) provides a very general starting point for the statistical mechanics of associating fluids and is exact so long as the system is pairwise additive, and
the intermolecular potential can be separated into a reference and association
­portion as in Eq. (1). The challenge is to approximate the graph sum Δc(o) (here we
assume the properties of the reference fluid are known).
A simple and widely used approximation of this sum forms the foundation of
Wertheim’s TPT [32, 33]. To start we decompose Δc(o) as

c

o

o
cˆk

(65)

k 1

where cˆk( o ) is the contribution for graphs which contain k associated clusters. For
example, graph d in Fig. 8 belongs to the sum cˆ2( o ), while the remaining graphs
all belong to cˆ1( o ). TPT can be defined as the neglect of all cˆk( o ) for k > 1 giving
c

o

o
cˆ1

(66)

The approximation Eq. (66) accounts for interactions within the associated
­cluster and between the clusters and the reference fluid, but not the interactions
between associated clusters. This approximation is accurate for cases in which
the pair correlation function between associating molecules is similar to that of
the reference fluid.

24

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 8. Examples of graphs in Δc(o) for the two site AB case. Wavy lines represent fAB bonds,
dashed lines represent eref bonds, and solid lines represent fref bonds.

Equation (66) accounts for multiple bonded association sites (graph f in Fig. 8
belongs to this class), cycles of association bonds (graph g in Fig. 8), multiple
bonds between two molecules (graph c in Fig. 8), and chains (trees) of association
bonds (graphs a, b, and e in Fig. 8). For the time being, it will be assumed that the
intermolecular potential and placement of association sites is such that contributions from cycle formation, multiple bonded association sites, and multiple bonds
between molecules can all be ignored, leaving only contributions for the ­formation
o
of chain and tree like clusters. With these restrictions ĉ1 can now be written as
ĉ1
o

o

o

cTPT 1

o

cTPT 2

o

cTPT 3 

(67)

where cTPT 1 is the first‐order contribution that contains all contributions for
­association between a pair of molecules (graph a and b in Fig. 8 are examples),
o
cTPT 2 is the second‐order contribution that contains information about the
­simultaneous association of three molecules (graph e in Fig. 8 belongs to this
o
o
class) etc. For the case of molecules with a single association site ĉ1
cTPT 1 .

25

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

If it can be assumed that steric hindrance between association sites is small, the
sum in Eq. (67) can be truncated at first order (TPT1) giving the simple result [32]
c

o

o
cˆ1

o

cTPT 1

1
2A

A

1 gref 12 f AB 12

2 d 1 d 2

B

B

(68)
(o)
It should be noted that while cTPT
1 only accounts for interactions between pairs
of molecules, all possible trees of associated clusters can be reproduced. However,
in first order, the only steric interactions are between nearest neighbors in the
cluster giving a theory that is independent of bond angles.
The combination of Eqs. (64) and (68) summarizes Wertheim’s first‐order
­theory for multiple association sites. The density parameters in Eq. (57) are determined by minimizing the free energy with respect to variations in these parameters.
Once an association model is specified, the free energy function and equations for
the density parameters can be derived from Eqs. (64) and (68).
Since the association sites on a molecule are independent at first order in
­perturbation, Chapman derived a closed‐form solution for mixtures of molecules
with any number of association sites [18, 34].

i

A Aref

i

k BT

i A

1

XA 1

i

ln X A 1

2

i

1
d 1
2

(69)

i

The fractions of molecules of species i that are not bonded at site A, X A , are
solved for self‐consistently by minimizing the free energy with respect to the
density parameters in Eq. (57) as
i
i

XA 1

i

1

A
i

1

1
j

1
j

1 gref 12

f AB 12 X B
B

j

j

2 d 2
(70)

For homogeneous fluids, the densities and X A(i ) parameters are no longer functions of position. The primary approximation in Wertheim’s perturbation theory is
that the unbonded site–unbonded site distribution function can be approximated by
the pair correlation function of the spherical reference fluid. Since the structure of
an associating fluid is known to be much different from that of a simple fluid, this
approximation might seem surprising. In the section on single association sites, we
stated that, for association near sphere contact, the packing fraction of the fluid
does not change with extent of association, and therefore the monomer–monomer
pair correlation function is assumed independent of the extent of association.
In this case, approximating the unbonded site–unbonded site correlation function

26

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

with the reference fluid pair correlation function is in reasonably good agreement
with molecular simulation results [18, 19, 35, 37, 38, 48]. Since association is short
ranged and commonly modeled as a square well, it is common to approximate the
integral in Eq. (70) using the pair correlation function at contact [18, 19, 34, 39].
Numerical solution of Eq. (70) is simplified since the X A(i ) parameters are the s­ olution
to an unconstrained minimization of the free energy [49]. Given that the association
sites are independent of one another at first order in perturbation, the monomer
­fraction of component i can be calculated from the product of the X A(i )’s.
This simple TPT1 result has proven to be a very powerful tool in the theoretical
description of associating fluids. Equations (69) and (70) have been shown to be
accurate for bulk fluids composed of hard spheres with one [19], two [19, 50], and
four association sites [38, 51]. The theory has also been shown to be accurate for
associating LJ spheres with one, two [37], and four [38] association sites and
associating molecules with a square well reference fluid [52]. Also, TPT1 has
been shown to accurately describe novel phase behavior in patchy colloid fluids
[50, 53–57]. In addition to spherical molecules, TPT1 has been shown to be
­accurate for associating chains [18] of tangentially bonded hard [48], LJ [37, 38,
48, 58, 59] and square well [52] spheres. When applying TPT1 to associating
chains, a chain reference fluid must be used, which is obtained using TPT1in the
complete association limit [18].
Given the accuracy of the theory in describing the properties of networking
forming fluids, one might imagine using the associating spheres as molecular
building blocks to build specific structures. By defining a mixture of molecules
with one, two, or more association sites that can only bond to specific sites on
other molecules, it is possible to define the structure that will form in the limit of
strong association. For example, if molecule 1 can only bond to the A site on
­molecule 2, and the B site on molecule 2 can only bond to the A site on molecules 3,
etc. until the A site on molecule (m − 1) can only bond to the B site on molecule m,
a linear chain of length m segments can form given a stoichiometric ratio of components and large enough association energy. In this way, by taking the limit of
infinite association energy, the perturbation theory produces an equation of state
for polyatomic molecules [18, 34] that is in reasonable agreement with molecular
simulation results for polyatomic molecules made of chains of hard spheres, LJ
chains, and square well chains [18, 33, 34, 37, 58–61].
The theory was further extended to mixtures of associating polyatomic molecules [18, 34, 37, 48, 62–64]. This requires the pair correlation function between
unbonded sites on the molecules. The simplest approximation is that the pair
­correlation function between unbonded sites is the same as that in the reference
fluid. This is valid as long as the angle between association sites is large enough
that the sites do not affect each other. One might question using the pair correlation function for spheres as an estimate of the pair correlation function between
association sites on chain‐like molecules. Some may suggest that this contradicts
the correlation hole effect seen in site–site distribution functions; however, the
correlation hole effect is largely a result of averaging over spherically symmetric

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

27

potentials. Over the small range of angles for which molecules can associate,
the approximation of the unbonded–unbonded distribution function with that
for spheres (particularly for hard sphere chains) has been shown to be a good
approximation [37, 38, 48, 58, 59].
When LJ spheres (or square well spheres) associate with each other, we expect
that the LJ (or square well) contribution to the internal energy (and free energy)
should be a function of the extent of association due to shielding of molecules in
the associated cluster. In other words, the dispersion contribution to the internal
energy should be different from that of a fluid of independent LJ spheres. The
TPT1 predicts the change in internal energy due to association or chain formation
in good agreement with molecular simulation results particularly for dense fluids
[37, 38, 48, 58, 59].
In addition to model systems, TPT1 has been widely applied in both academia
and industry as an engineering equation of state for hydrogen bonding fluids. The
equation of state resulting from the extension of Wertheim’s TPT1 to mixtures of
associating polyatomic molecules is called the statistical associating fluid theory
(SAFT) [18, 34, 37, 59, 62, 63]. Many of the SAFT forms [63, 65–67] have found
wide application among scientists and engineers in academia and industry to
describe the properties and phase behavior of solvents to associating polymers as
well as patchy colloids.
Though widely applied, TPT1 is far from perfect with a number of limitations
resulting from the simplifying approximations employed. These approximations
are summarized as follows:
1. Single chain approximation—Neglects all graphs with more than one associated cluster. This is TPT, which assumes that the structure of unbonded
sites in the fluid is similar to that of the reference fluid. The single chain
approximation will fail, for instance, for fluids with a nematic phase [68].
2. Singly bondable association sites—Assumes each association site saturates
after sharing in a single association bond. This approximation is not valid
for patchy colloids with large patch sizes [43].
3. No multiple bonding of molecules—Assumes that any two molecules can
share at most one association bond. Carboxylic acids [69] and water [70]
are known to violate this condition.
4. No cycles of association bonds—Only chains and trees of association
bonds are accounted for. Cycles are irreducible and cannot be reproduced
in TPT1. It is well known that hydrogen fluoride exhibits significant ring
formation [71].
5. No steric hindrance between association sites—All contributions to the
irreducible graph sum with more than a single association bond were
­
neglected. For this reason, association at one site is independent of association at all other sites. Most polyfunctional associating molecules will exhibit
some degree of steric hindrance between association sites.

28

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

6. Association is independent of bond angles—There is no information in
TPT1 on location of association sites. This is intimately related to approximations 3–5 above.
7. In Wertheim’s multi‐density formalism, pairwise additivity of the pair
potential was assumed. Most polyfunctional hydrogen bonding molecules
exhibit some degree of bond cooperativity (non‐pairwise additivity) [72].
Hydrogen bond cooperativity is particularly important for hydrogen
­fluoride [71].
In recent developments, a number of these approximations have been relaxed.
The remainder of this chapter will review some of these extensions of the theory
and development of a molecular DFT for inhomogeneous fluids.
VI. THE TWO‐SITE AB CASE
In Section V it was shown how Wertheim’s multi‐density approach could be used
to develop an equation for associating fluids with an arbitrary number of association sites provided a number of assumptions were satisfied. The simplicity of
the TPT1 solution results from the fact that the solution is that of an effective
two‐body problem. Only one bond at a time is considered. This allows the theory
to be written in terms of pair correlation functions only, as well as obtain analytical solutions for the bond volume. Moving beyond TPT1 is defined as considering
irreducible graphs that contain more than one association bond.
The simplest case to illustrate this extension is for molecules with a single type
A and type B association site, where the center of these two sites is separated by
the angle αAB. We loosely call αAB the bond angle. Figure 9 gives examples of
association into linear chains for both cases of large (case a), and small (case b)
αAB. In case a, since the sites are widely separated, association at each site will be
(a)

αAB

(b)

Figure 9. Examples of linear triatomic clusters for two angles αAB. (See insert for color representation of the figure.)

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

29

independent. For case b this is not the case. Due to the small bond angle, association
at one site may block a third molecule coming to associate with the other site.
There is steric hindrance. Further, for simplicity, we will assume that there is only
attraction between unlike sites. That is there are AB attractions, but no AA or AB
attractions. This could serve as a primitive model for a bifunctional hydrogen
bonding molecule such as hydrogen fluoride, monomer in a supramolecular
­polymer, globular protein, or patchy colloid.
To extend the TPT1 results, such that systems with small αAB can be considered, we need to introduce the effect of steric hindrance. However, inclusion of
steric hindrance alone is not sufficient due to the fact that as αAB is decreased, the
probability of forming small rings increases. Indeed, as has been shown [39, 73],
for small bond angles rings become the dominant type of associated cluster. In
what follows we outline the methodology to include these types of higher‐order
interactions into an accurate equation of state which explicitly depends on αAB.
A. Steric Hindrance in Chain Formation
To incorporate steric effects in chain formation we must employ Wertheim’s
TPTM. In TPTM, Δc(o) is approximated by considering all chain diagrams that
contain a single chain of M or less association bonds and is given as [33]
c

o

M

o
cˆ1

o

(71)

cTPTn
n 1

(o)
where cTPTn
is the nth‐order contribution (involves chains of n association bonds)
and is given by (assuming a homogeneous fluid)
o

cTPTn
V

A

B

n 1
o
n

(72)

I

The integrals In are given by
In

1
n

f AB 12

f AB n, n 1 Gref 1 n 1 d 2

d n 1

(73)

where Ω = 8π2 and the fAB(ij) are the association Mayer functions. Wertheim
defines the functions Gref (1 n 1) as “the subset of graphs in gref (1 n 1) such
that combining them with the chain produces an irreducible graph; gref(1 … s)
denotes the s particle correlation function of the reference system” [33]. This
means, for instance, that in a TPT2 the contribution ΔcTPT2 will include the triplet
correlation function gref(123), but one must subtract off the contribution from the
(o)
first‐order term cTPT
1 to keep from double counting. We then obtain the Gref(1 … s)
by summing gref(1 … s) and all products of gref’s obtained by partitioning 1…s into

30

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

subsequences which share the switching point and associating a negative 1 with
each switching point [33]. A few examples include the following:
Gref 12

gref 12

Gref 123

gref 123

Gref 12234

gref 12 gref 23

gref 1234 gref 123 gref 34
gref 12 gref 23 gref 34

(74)

gref 12 gref 234

The general idea of TPTM is then to build up chains by adding in higher‐order
contributions and subtracting off lower‐order contributions.
(o)
Returning our attention to case a in Fig. 9, all cTPTn
with n > 1 can be neglected
due to the fact that there will be little steric hindrance between sites. For case b the
situation is more complex as there are steric effects between association sites. To
include these steric effects we must go to a minimum of TPT2. For clarity we will
now assume a hard sphere reference system with CSW association sites as given
in Eq. (4). In TPT2 we keep all contributions in Eq. (71) up to and including M = 2.
The integrals I1 and I2 are evaluated as
I1
I2

1

f AB 12 gHS 12 d 2

1
2

f AB

AB

2
AB

f AB 12 f AB 23 GHS 1223 d 2 d 3

2
AB

yo 123

1
(75)

In Eq. (75) the triplet function yo(123) is defined as
yo 123

yHS 123
yHS 12 yHS 23

(76)

and
represents the average over all states where both sites on the molecule are
bonded
yo 123

AB

12

AB

23 eHS 12 eHS 23 eHS 13 yo 123 d 2 d 3
2
b

(77)
The steric effects in TPT2 are now wholly included in this average. Employing
the mean value theorem Eq. (77) can be further simplified as
yo 123

ch

yo 123

AS

(78)

31

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

where Ψch is given by

ch

The probability that if the positions of two molecules are generaated
such that site A on molecule 1 is bonded to site B on molecule 2 and
that site B on molecule 3 is bonded to site A on molecule 2, that there
is no overlap between molecules 1 and 3
(79)

and ⟨ ⟩AS represents the average over all associated states of the cluster in which
there is no hard sphere overlap. In short, the average in Eq. (77) includes states
with hard sphere overlap in the normalization (κνbΩ)2, while the average ⟨ ⟩AS does
not with the normalization Ψch(κνbΩ)2.
Equation (78) splits the TPT2 contribution into a geometric contribution Ψch,
which accounts for the decrease in bond volume due to steric hindrance, and a
density dependent contribution ⟨yo(123)⟩AS which accounts for the effect of bond
angle on the fluid packing. At this point ⟨yo(123)⟩AS has not been evaluated;
­however, this average should be relatively easy to evaluate using the fitting f­ unction
of Müller and Gubbins [74] in combination with Monte Carlo integration.
A desirable limit of this approach is that the density of molecules bonded at
both sites ρAB vanishes as the angle AB 0. In this limit it becomes impossible
for a molecule to be bonded at both sites. To check this limit we set ch 0 and
evaluate ρAB through Eq. (62) as (enforcing A
B).
AB
o

c AB
AB

c A cB

A

o

2
AB

2

2

2
A

o

3
AB

0

(80)

0

As can be seen, TPT2 does not satisfy this limit showing that TPT2 does not
include full steric information between the two association sites. This deficiency
results from the way in which reducible graphs are created from irreducible graphs
through elimination of the density parameters σA and σB.
To include full steric information it is necessary to consider perturbation to
infinite order. Unfortunately, we know very little of correlation functions
beyond the triplet case. To proceed further we note that the TPT2 contribution
given in Eq. (78) is split into a purely geometric contribution and a density
dependent contribution which depends on knowledge of the triplet correlation
function. A simplification of Eq. (78) would be to assume that the average
⟨yo(123)⟩AS was approximately unity giving
yo 123

ch

(81)

Equation (81) will account for blocking effects between the two association sites
by correcting the decrease in the total number of associated states of the cluster.

32

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

Equation (81) will be exact in the low‐density limit, and should be a fair approximation at higher densities, becoming more accurate as αAB increases.
To apply the theory as an infinite‐order perturbation theory we must approximate gHS (1 n 1) in such a way that the infinite sum over all chain graphs can be
performed. For this we take the same approach as in Section IV.B and approximate
the multi‐body correlation functions with Eq. (36). Like Eq. (81), the superposition
Eq. (36) treats higher‐order effects in a density‐independent way by incorporating
purely geometric constraints in the association model. With Eq. (36), the infinite
sum in Eq. (71) for M
can be approximated as follows:
o

cch
V

A

1

f AB
f AB

B

1

(82)
o

Now we use Eq. (82) to evaluate the density of molecules bonded at both sites as
follows:
AB

A
ch

o

1

1

ch

f AB
f AB

2

(83)
o

It is easy to see from Eq. (83) that as ch
0 the density of molecules bonded
twice vanishes. This shows that, in contrast to TPT2, the theory properly accounts
for blocking effects between the two association sites.
B.

Ring Formation

In addition to steric effects, when bond angles are small as in case b, ring formation can become significant. Where now we are assuming association sites are
singly bondable and rings are composed of cycles of AB bonds as depicted in
Fig. 10. It was Sear and Jackson (SJ) [75] who were the first to introduce contributions for association into rings. In this approach the associated rings were treated

Figure 10. Examples of associated rings. (See insert for color representation of the figure.)

33

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

ideally such that non‐adjacent neighbors along the ring can overlap. The probability
that a chain of colloids was in a valid ring state was approximated by the expression of Treolar [76] for the distribution of the end‐to‐end vector in a polymer
chain. In this approach any dependence on αAB is neglected when in reality αAB
plays a dominant role in determining if association into rings will occur. A recent
study using lattice simulations has shown that ring formation is strongly dependent on αAB [73]. For instance, it is impossible to form 4‐mer rings (and satisfy the
one‐bond‐per‐site condition) when AB 180; however, decreasing αAB to 90° this
type of ring would be possible. Tavares et al. [77] explored the possibility of ring
formation in two patch colloid fluids with AB 180 by extending the approach
of SJ [75] and found that to achieve appreciable ring formation the parameters of
the association potential had to be chosen such that the one‐bond‐per‐patch condition would be violated. To correct for this in the simulations they used a model
that restricts bonding to at most one bond per patch [78].
Recently, Marshall and Chapman [39] extended TPT to account for the effect
of αAB on ring formation. To allow for rings of all sizes the ring contribution to the
graph sum is split into contributions from rings of size m
o

cmring

cring

(84)

m 3

where cmring is the contribution for rings of size m. The contributions cmring are
obtained by summing over all graphs that contain a single ring of m association
bonds.
cmring
V

m

f

o AB

m 8

2

m 1

AB

12

AB

m 1, m

AB


1, m gHS r1

 
rm dr2 d

2


drm d

m

(85)
Pictorially, graph g in Fig. 8 represents the low‐density limit of c4ring . Evaluation
of Eq. (85) is much more difficult than for the chain contribution due to the presence of the additional association bond that closes the ring, rendering the graph
 
irreducible irrespective of the superposition used to approximate gHS (r1 rm ). To
 
proceed, we consider the following simple superposition of gHS (r1 rm ):

gHS r1


rm

bonded pairs
i,j

yHS rij

eHS rlk

(86)

all pairs
l ,k

In Eq. (86) each pair of spheres in the ring receives an eHS(rlk), meaning the ring is
fully self‐avoiding. This is in contrast to a similar approximation for chain
­formation in Eq. (36), for which the steric effects are limited to first and second

34

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

neighbors along the chain, which ultimately allows for the factorization of the
chain into dimer and triplet contributions. Here, since the ring graph cannot be
factorized regardless of the superposition used, the additional connectivity of
Eq. (86) imposes no penalty. Using Eqs. (85) and (86) Marshall and Chapman
[39] arrived at the following form for cmring
f AB

cmring
V

o

gˆ HS K

md

m
m

(87)

3

where
2 p gHS d

gˆ HS

rc / d 1

(88)

p

and Γ(m) is proportional to the partition function of an isolated ring of size m.
Both Γ(m) for m = 3–10 and ψ have been evaluated numerically for CSW association sites with potential parameters c 27, rc 1.1d over the full bond angle
range [39]. Figure 11 gives results for Ψ and Γ(m) for rings of size m = 3 and 4.
As expected, Ψ vanishes for small αAB due to steric hindrance, and becomes unity
for large αAB when association at one site no longer interferes with the ability of
the other site to bond. The ring integrals Γ(m) are peaked around an optimum bond
angle for ring formation, and the maximums of Γ(m) decrease and shift to larger
bond angles as m increases. It is the relative magnitudes of these geometric
­integrals which control the types of associated structures that exist.
1

Ψ

0.8

Γ (3)

0.6
0.4

Γ (4)

0.2
0
0

30

60

90
αAB

120

150

180

Figure 11. Numerical evaluation of geometric integrals for the two‐site case [79]. Ψ represent
the probability that there is no core overlap in a triatomic chain and ring integrals Γ(m) are proportional
to the partition function of an isolated ring of size m.

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

35

The Helmholtz free energy for the two‐site case with steric hindrance and ring
formation can be developed from the results of Sections V and VI as follows:
A A HS
NkBT

A A HS
NkBT

(89)

aHO
TPT 1

Where the first term on the right‐hand side represents the TPT1 solution obtained
from Eq. (69) and ΔaHO represents the higher‐order correction given by
aHO

ln

Xo
XA XB

m

X mring
3 m

(90)

Here, X mring is the fraction of molecules in rings of size m given by [39]
X

ring
m

f AB

o

gˆ HS K

m
m

(91)

d3

From Eq. (90), we see the first‐order limit will be obtained when αAB is large
enough that ring formation becomes improbable X mring 0, and steric hindrance
between sites is small enough that association at one site is essentially independent of the other. When association sites are independent due to the lack of ring
formation and steric hindrance, the relation X o X A X B holds. Under these
­conditions, ΔaHO will vanish and a treatment in TPT1 will be justified.
Figure 12 compares theory predictions and Monte Carlo simulations for the
fraction of spheres bonded in rings of size m = 3–4, and the fraction of spheres that
1

Fraction

0.8

η = 0.3

0.6

X3ring

αAB = 60°

0.4
X4ring

X2c

0.2
0
2

4

6

8

10

12

ε*

Figure 12. Comparsion of TPT (curves) and simulation results (symbols) for the fraction of
two‐site molecules in trimer rings, 4‐mer rings, and chain centers X2c (chains of all sizes) versus reduced
association energy. Potential parameters are c 27 and rc 1.1 d. The figure is modified from Ref. 39.

36

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

are bonded at both sites in chains of any size X2c (spheres in the middle of a chain)
[39]. The comparison is made for an angle AB 60 and a liquid‐like packing
fraction of η = 0.3. For low *
AB / k BT , the majority of fully bonded colloids are
in the center of triatomic chains resulting in X2c being the dominant contribution.
Increasing ε*, X3ring rapidly becomes the dominant type of associated cluster in the
fluid forcing a maximum in X2c which becomes very small in strongly associating
systems. The fraction X 4ring shows a nearly linear increase with ε*, overtaking X2c
near *~ 9.5.
As can be seen, theory and simulation are in excellent agreement. This may
come at a surprise to many. After all, the superpositions Eqs. (36) and (86) neglect
any density‐dependent packing effects beyond nearest neighbors in the cluster.
What this shows is that the accuracy of the theory is largely (not completely)
determined by getting the geometry right. That is, if the number of associated
states of an isolated cluster can be calculated, TPT can give good predictions over
a wide range of densities. Of course, at some density, packing effects must become
important. Marshall and Chapman [39] showed that the theory loses accuracy
somewhat for η = 0.4; however, even at this high density, the theory still performed well.
C. Bond Cooperativity
In the previous sections we have shown how TPT can be extended to describe
a wide variety of associating fluids. In each case, the distribution of associated
­clusters and the resulting equation of state were strongly dependent on a delicate balance between the energetic benefits of association and the resulting
entropic penalty. In each case it was assumed that the total system energy is
pairwise additive, there is no bond cooperativity. In nature, hydrogen bond
cooperativity arises from the fact that when a multi‐functional hydrogen
­bonding molecule forms hydrogen bonds, the polarization of the molecule is
necessarily increased [72]. As has become increasingly apparent in recent
years, hydrogen bond cooperativity plays a significant role in many physical
processes. Both hydrogen fluoride (HF) [71] and alcohols [80] have been
shown to exhibit strong hydrogen bond cooperativity. In addition, hydrogen
bond cooperativity has been shown to stabilize peptide hydrogen bonds [81].
Indeed, it is believed that all polyfunctional hydrogen bonding molecules
exhibit some degree of bond cooperativity [72].
The assumption of a pairwise additive association potential, Eq. (2), forms the
foundation from which Wertheim’s multi‐density formalism is built. Strictly
speaking, TPT cannot be applied to non‐pairwise additive association potentials
in a rigorous way. However, it was recently [79, 82] demonstrated how the same
ideas used to develop Eq. (82) can be applied to develop an equation of state
for hydrogen bonding fluids that exhibit bond cooperativity. The approach
employs the potential model of SJ [83] for a fluid composed of Np hard spheres of

37

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

diameter d with two association sites A and B (only AB attractions and singly bondable sites) with a total energy composed of pairwise and triplet contributions [83]
U 1

1
2

Np

2
as

rij

HS

1
6 i, j ,k

ij

i, j

3
as

(92)

ijk

The term as( 2 ) (ij ) is the pairwise contribution given by Eqs. (2), (4) and
is the triplet association contribution:
2
as
3
as

ijk

1

ij
2

1

ij

AB
AB
AB

ij

AB

ji
ki

BA
BA

jk
kj

BA

BA
BA

BA

ji
ki

AB
AB

(ijk )

(93)

ij

ik

BA

(3)
as

ij

AB

ik

jk
kj

The triplet contribution as(3) serves to add a correction ( ( 2 ) (1) ) for each sphere
bonded twice. With this potential an associated chain of n spheres will have a
cluster energy:
n
ch

1

2

n 2

.

(94)

Strictly speaking, Wertheim’s formalism cannot be rigorously applied to this
potential; however, employing the same ideas used to develop Eq. (82), a very
(o)
accurate approximation of cch
can be derived for this potential model. The
­derivation begins with the observation that the energy given by Eq. (94) can be
partitioned such that the first association bond in a chain receives an energy ε(1),
while each remaining bond in the chain receives an energy ε(2). This partitioning is
illustrated in Fig. 13. Partitioning the bond energies in this manner allows one to
replace the product of association Mayer functions in Eq. (73) with the following:
f AB 12

f AB n, n 1

1

2

f AB 12 f AB 23

– ε (2)

– ε (1)

Figure 13.
given by Eq. (93).

– ε (2)

2

f AB n, n 1

(95)

– ε (2)

– ε (2)

Diagram of bond energy distribution in associated chain with bond cooperativity as

38

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

1
X1

Fractions

0.8

X2

0.6
0.4
X0
0.2
0

0

2

4

6

8

10

ε (2) / kBT

Figure 14. Comparison of theory and simulations for the bonding fractions (fraction of spheres
bonded k times) of two‐site‐associating spheres with bond cooperativity. The pairwise association
energy and density are held constant at ε(1) = 7 kBT and ρ = 0.6d3. Curves give theory predictions and
symbols give Monte Carlo simulation results. The figure is modified from Ref. 79.
(k )
(k )
/ kBT ) 1) AB (12) . Combining
where the Mayer functions f AB (12) (exp(
Eqs. (73) and (95) effectively map the non‐pairwise additive association potential
onto Wertheim’s multi‐density formalism. Following this approach it was shown
that when the transformation Eq. (95) is combined with the superposition
­approximation (36), Eq. (71) can then be summed in the limit of an infinite‐order
perturbation theory as [79, 82] follows:

o

cch
V

A

1

f

1
AB

1
B AB

f

(96)

2

f AB

o

Equation (96) is a very simple result which accounts for both steric hindrance
and hydrogen bond cooperativity in two‐site‐associating fluids. When hydrogen
(1)
(2)
bonding is non‐cooperative f AB
f AB
, and Eq. (96) simplifies to Eq. (82). To
demonstrate the accuracy of this approach, Fig. 14 compares theory predictions to
Monte Carlo simulations for the fraction of molecules which are m
­ onomers Xo,
bonded once X1 and bonded twice X2. For simplicity we consider the case of a
bond angle αAB = 180°, pairwise association energy ε(1) = 7 kBT and density ρ = 0.6d3.
For ( 2 ) 0, there is no energetic benefit for a sphere to bond twice which results
in X 2 0. Increasing ε(2) results in a steady increase in X2 and the fractions X1 and
Xo remain nearly constant until ε(2) ~ 5kBT at which point they decline sharply.

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

39

The lack of fully bonded molecules for small ε(2) is due to the fact that there is little
energetic benefit to a molecule to becoming fully bonded, while the entropic penalty must still be paid in full. For large ε(2) the energetic push for molecules to become
fully bonded overpowers the entropic penalty resulting in X 2 1 for large ε(2). Theory
and simulation are in near‐perfect agreement.
In addition to chains, rings can also be included in the bond cooperative
(2)
­perturbation theory through the simple transformation f AB
f AB
in Eq. (85).
Including the effect of bond angle, bond cooperativity, and ring formation,
Marshall et al. [82] were able to show that both bond angle and bond cooperativity
play a huge role in the types of associated clusters which are formed. In agreement
with detailed quantum calculations [71], it was shown that bond cooperativity
favors ring formation.
VII. SPHERICALLY SYMMETRIC AND DIRECTIONAL
ASSOCIATION SITES
Throughout this chapter association has been defined as being between two
­molecules that must be positioned and oriented correctly for association to occur.
That is, both molecules participating in the association bond have directional
association sites. Another common case would be an association interaction
between two molecules where one has a directional association site, while the
other has a spherically symmetric association site. This type of interaction could
describe ion–water solvation or mixtures [42] of patchy and spherically symmetric colloids.
The single‐component version of this type of interaction is the Smith and
Nezbeda [84] model of associating fluids. This model considers a spherical core
with a single directional bonding site. The directional association sites are singly
bondable and the spherical cores are treated as spherically symmetric association
sites, with the maximum number of bonds the spherical core can receive being
determined by steric constraints. Wertheim [85] developed an integral equation
theory for this model of associating fluids, which was later solved analytically by
Kalyuzhnyi and Nezbeda [86]. Also, we may draw parallels with the study of
highly asymmetric electrolyte solutions. These solutions contain large polyions
and small single‐charge counterions. Previous multi‐density integral equation
theory studies of these solutions [87–89] have treated the counterions as singly
bondable and the maximum number of times the polyion ion can bond is unrestricted and determined by steric constraints.
Very recently Marshall and Chapman [90, 91] developed a new TPT to model
mixtures of these types. Specifically, they considered a mixture of molecules with
directional association sites (d molecules) and molecules with a single spherically
symmetric associations site (s molecules). The d molecules have CSW association
( d ,d )
sites as given in Eq. (4) with an association energy between d molecules AB
,

40

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

the s molecules do not attract other s molecules, and the association between s
­molecules and d molecules is governed by the following association potential
s ,d
A

s ,d
A

12

, r12

rc and A
otherwise

0

(97)

c

which states that if s molecule1 and d molecule 2 are within a distance rc of each
other, and the d molecule is oriented such that the angle between the site A orientation vector and the vector connecting the two segments θA is less than the critical
angle θc, the two molecules are considered bonded and the energy of the system is
decreased by a factor A( s , d ).
In what follows, attention will be restricted to the case that both s and d molecules have a hard spherical core Eq. (3), equal diameters, and the d molecules will
be restricted to having a single association site. For this case, this mixture can be
treated as a binary mixture of associating molecules in Wertheim’s two‐density
formalism outlined in Section IV. As done throughout this chapter, we will consider a perturbation treatment with a hard sphere reference fluid. Like previous
cases, the challenge is determining the graph sum Δc(o).
For the current case, the s molecule is a single spherical association site
which can clearly not be modeled in the single bonding condition. The maximum number of bonds is simply the maximum number of d molecules nmax that
can pack in the s molecule’s bonding shell d r rc . To account for all association possibilities, we will have to include contributions for each association
­possibility explicitly (one s molecule with one d molecule, two d molecules,
three d molecules, etc.). To accomplish this, we decompose Δc(o) as follows:
c

o

nmax

o

(98)

cn
n 1

where cn( o ) is the contribution for n directional molecules bonded to a s molecule.
Figure 15 gives a pictorial representation of Eq. (98).

Δc1(o)
Δc(o) =

Δc2(o)
+

Δc3(o)
+

+ ... + Δc

(o)
nmax

Figure 15. Illustration of graph sum contributions Eqs. (98) and (99) for a binary mixture of s
molecules and d molecules.

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

41

The development of cn( o ) follows a similar path to the chain contribution for
doubly bonded association sites Eq. (39). The final result is (in a slightly different
form) [91]
o

cn
V

1
n!

s
o

n

n

Zn

s ,d

(99)

where


1

s, d

exp

kB T

d

1

o

yHS d

(100)

with o( s ) being the monomer density of s molecules, o( d ) the monomer density
of directional molecules and δ(n) being a second‐order correction to the superposition approximation Eq. (86). Here Z n( s ,d ) is the cluster partition function for an
isolated cluster of n directional molecules bonded to a single s molecule.
As with the other perturbation theories discussed in this chapter, the key quantity in the theory is the cluster partition function which enumerates the number of
associated states for which an isolated associated cluster can exist. For this case
the cluster partition function is simplified as follows:
Zn

s ,d

vb

n

P

n

(101)

The term P(n) is the probability that if n directional molecules are randomly generated in the bonding shell of the s molecule that there is no hard sphere overlap. An
interesting feature of this result is that the square root
1 cos c / 2 appears
in the equations, as opposed to κ as in Eq. (46). The reason for this is that when
one d molecule forms an association bond to an s molecule, there is a total orientation entropic penalty of kB ln
. On the other hand, when two d molecules form
association bond, there is a total orientation entropic penalty of kB ln κ. The difference in the two cases lies in the fact that the s molecules do not pay a penalty in
decreased orientational entropy when an association bond is formed. The number
nmax is defined as the largest integer n for which the probability P(n) is nonzero. For
bonding at hard sphere contact nmax = 12.
Figure 16 compares theoretical predictions and Monte Carlo simulation results
for a mixture d and s molecules with an association energy ( s ,d ) 7kBT at both
low and high densities. The average “solvation number” of s molecules n , or average number of bonds per s molecule, is plotted against the number fraction of s
molecules x(s). For each case, n increases with decreasing x(s) reaching a maximum
when x ( s ) 0. This is due to the fact that when x(s) is small, there is an abundance
of d molecules available to bond to the s molecules. As x(s) is increased, n decreases

42

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

10
8
6
n–

ρ* = 0.7

4
2
0

ρ* = 0.2
0

0.2

0.4

0.6

0.8

1

x (s)

Figure 16. Comparison of theory (curves) and Monte Carlo simulation (symbols) results for
the average solvation number of s molecules versus number fraction of s molecules at an association
energy ( s ,d ) 7kBT and densities ρ* = ρd = 0.2 (circles, dashed curve) and 0 (triangles, solid curves).
The figure is modified from Ref. 91.

because there are less d molecules available for association due to a decreased
fraction of d molecules and competition with other s molecules. Overall, theory
and simulation are in good agreement. In addition to average solvation numbers,
it was shown that the theory accurately predicts the distribution of s molecules
bonded n times [91].
Going beyond the single site case, the theory was recently extended such that
the d molecules can have an arbitrary number of association sites [90]. In this
approach the interaction between s molecules was also that of the hard sphere
reference fluid. To add spherically symmetric attractions (square well, LJ, etc.)
between s molecules, one simply needs to employ the appropriate reference
­system (square well, LJ, etc.). Work is currently under way to employ this association theory as a model for ion–water solvation.
VIII.

DENSITY FUNCTIONAL THEORY

Thus far, the focus of this review has been homogeneous fluids. For many interesting phenomena observed in biological and soft material systems, the micro‐ or
mescoscale structure determines the properties of the system. DFT provides
a ­valuable tool to predict mesoscale structure and interfacial properties assuming
a suitable free energy functional can be developed. Excellent reviews of DFT for
associating molecules have been written [92–94], so only a brief introduction will
be provided here.
Calculating the inhomogeneous fluid structure of associating molecules based
on Wertheim’s perturbation theory was first proposed by Chapman [18]. At that
time, accurate molecular DFTs for non‐polar spherical molecules (e.g., Tarazona’s

43

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

weighted density approximation [95] for hard spheres) had been developed. Two
methods were suggested to include a perturbation for association in a DFT for
nonpolar spheres. One approach was to include association using the local density
approximation [52, 96] or using weighted densities of the bulk free energy due to
association in a similar way to that used to create a hard sphere free energy
­functional by Tarazona or later Rosenfeld [97]. The second approach was to take
advantage of the fact that Wertheim’s TPT for associating molecules was already
written in the form of a free energy functional (see Eqs. 69 and 70). A challenge
with the associating free energy functional is to approximate the inhomogeneous
cavity correlation function required by the theory.
Kierlik and Rosinberg [68, 98, 99] were the first to apply Wertheim’s theory in
the form of a free energy functional to produce a DFT for non‐associating polyatomic molecules. As input to the theory, they estimated the cavity correlation
function from a first‐order functional Taylor series around the homogeneous result
[99]. Results were in good agreement with molecular simulation results for hard
sphere chains.
Segura et al. [51, 100] used two perturbation approaches to produce DFTs for
associating molecules with one, two, and four association sites. In the first
approach, they calculated the bulk free energy due to association at the same
weighted densities as from Tarazona’s DFT [95] for hard spheres. The theory
results were in good agreement with molecular simulation for associating hard
spheres near a hydrophobic surface. Further studies have showed accurate results
in comparisons with molecular simulation for mixtures and confined associating
molecules with hydrophobic and hydrophilic surfaces as described in the reviews
[92, 93]. The same approach has been applied using various weighted densities or
fundamental measure theory with similar accurate agreement with molecular simulation. Still other studies have used gradient theory or a local density approximation for vapor–liquid interfaces and shown good agreement with interfacial
tension data [92, 93, 101, 102]. Extensions of the weighted bulk association free
energy approach of Segura et al. have resulted in an accurate DFT for polyatomic
molecules [103–106] by taking the limit of complete association in the bulk as
described earlier.
The second approach of Segura et al. [51] was to use the free energy functional
of Eqs. (69) and (70) as a perturbation to a hard sphere DFT. To minimize the
system free energy requires the functional derivative of Eq. (69) with respect to
the singlet density. The result is
Aex ,assoc
j r

ln

j
A

r

j

A

1
2

dr1dr2

N

N
i

i 1 k 1

r1

k

i
A

r2
A

i

B

k

r1

k
B

r2

ik
AB

r1,r2
j r

(102)

44

BENNETT D. MARSHALL AND WALTER G. CHAPMAN

Based on Eq. (102), Jain et al. [107] proposed a DFT for heteronuclear polyatomic molecules by taking the limit of the free energy functional for complete
association of a fluid of spheres. Interestingly, in this limit, the theory automatically corrects the ideal gas free energy functional from an ideal gas of spheres to
predict the exact free energy and density distribution for an ideal gas chain.
Bymaster and Chapman [108] have shown that, in addition to associating spheres,
Eq. (102) is applicable to associating polyatomic molecules. Results based on this
DFT for associating molecules and non‐associating polyatomics are in good
agreement with molecular simulation results for associating molecules near a
hydrophobic surface [109], associating grafted polymers [110], surfactants [111],
and mixtures of associating polyatomics with intermolecular and intramolecular
association [112]. For further information, we recommend reviews available in
the literature as well as more recent literature [92, 94, 110, 113–116].

IX.

CONCLUDING REMARKS

In this chapter, the basics of Wertheim’s and Andersen’s cluster expansions for
associating fluids have been reviewed, specifically focusing on thermodynamic
perturbation theory (TPT). Despite the severe approximations made in TPT, the
approach yields equations of state that accurately reproduce simulation data, and
are widely used for modeling the properties and phase behavior of solvents to
polymers in the engineering community. As was shown throughout the chapter,
TPT can be applied to develop theories for complex associated structures with
steric effects, so long as the number of associated states of an isolated cluster can
be calculated. Of course, this TPT approach will fail once the correlation f­ unctions
of the real fluid begin to deviate from that of the reference fluid; for instance,
association of monomers into rigid chains resulting in a nematic transition.

Acknowledgments
The authors are grateful for the financial support of The Robert A. Welch
Foundation (Grant No. C-1241). Amin Haghmoradi is acknowledged for helpful
discussions.
REFERENCES
1. Weeks, J. D.; Chandler, D.; Andersen, H. C. The Journal of Chemical Physics 1971, 54, (12),
5237–5247.
2. Barker, J. A.; Henderson, D. The Journal of Chemical Physics 1967, 47, (11), 4714–4721.
3. Chandler, D.; Andersen, H. C. The Journal of Chemical Physics 1972, 57, (5), 1930–1937.

THERMODYNAMIC PERTURBATION THEORY FOR ASSOCIATING MOLECULES

45

4. Jeffrey, G. A. An Introduction to Hydrogen Bonding. Oxford University Press New York: 1997;
Vol. 12.
5. Dill, K. A. Biochemistry 1990, 29, (31), 7133–7155.
6. Whitesides, G. M.; Grzybowski, B. Science 2002, 295, (5564), 2418.
7. Bianchi, E.; Blaak, R.; Likos, C. N. Physical Chemistry Chemical Physics 2011, 13, (14),
6397–6410.
8. Economou, I. G.; Donohue, M. D. AIChE Journal 1991, 37, (12), 1875–1894.
9. Dolzalek, F. Zur Theorie der Biniiren Gemische und Konzentrierten Loungen. Z. Phys. Chem.
1908, (64), 727–747.
10. Heidemann, R. A.; Prausnitz, J. Proceedings of the National Academy of Sciences 1976, 73, (6),
1773–1776.
11. Ikonomou, G.; Donohue, M. AIChE Journal 1986, 32, (10), 1716–1725.
12. Panayiotou, C. G. The Journal of Physical Chemistry 1988, 92, (10), 2960–2969.
13. Veytsman, B. Journal of Physical Chemistry 1990, 94, (23), 8499–8500.
14. Panayi