Main
Advances in Chemical Physics, Volume 160
Advances in Chemical Physics, Volume 160
Stuart A. Rice, Aaron R. Dinner
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
 MeanField Electrostatics Beyond the PointCharge Description
 First Passage Processes in Cellular Biology
 Theoretical Modeling of Vibrational Spectra and Proton Tunneling in HydroenBonded Systems
Categories:
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
Your tags:
1

2

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: 589935 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 nearspherical 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 FAA (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 FAA as FAA 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 FAA 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 FAA (12) is bounded as follows: FAA 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: FAA 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 FAA 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 FAA , 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 reasonably 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 FAA 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 FAA 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. C1241). 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