We extend Wertheim’s thermodynamic perturbation theory to derive the association free energy of a multicomponent mixture for which double bonds can form between any two pairs of the molecules’ arbitrary number of bonding sites. This generalization reduces in limiting cases to prior theories that restrict double bonding to at most one pair of sites per molecule. We apply the new theory to an associating mixture of colloidal particles (“colloids”) and flexible chain molecules (“linkers”). The linkers have two functional end groups, each of which may bond to one of several sites on the colloids. Due to their flexibility, a significant fraction of linkers can “loop” with both ends bonding to sites on the same colloid instead of bridging sites on different colloids. We use the theory to show that the fraction of linkers in loops depends sensitively on the linker end-to-end distance relative to the colloid bonding-site distance, which suggests strategies for mitigating the loop formation that may otherwise hinder linker-mediated colloidal assembly.

It has been over 35 years since Wertheim developed a thermodynamic perturbation theory (TPT) for fluids of molecules that associate through strong directional attractions.1–4 Wertheim’s TPT has proven to be powerful and broadly applicable, capturing the thermodynamic behavior of fluids comprising, e.g., small molecules with chemically reactive sites, associating macromolecules, such as DNA nanostars,5,6 and patchy colloidal particles (“colloids”).7–10 Indeed, Wertheim’s TPT is a foundation for statistical associating fluid theory (SAFT),11–14 which has been widely used to develop equations of state for complex fluid mixtures.15 

In practice, many applications of Wertheim’s TPT, including SAFT, adopt a first-order approximation that treats association as pairwise bonding into tree-like networks and is often referred to as TPT1.2,4 We previously used TPT1 to predict the phase behavior of mixtures of colloids and difunctional chain molecules (“linkers”) in which linker ends associate with sites on the colloids through strong reversible bonds (Fig. 1).16 Homogeneous colloid–linker mixtures of this type lose thermodynamic stability for a range of colloid or linker concentrations,16,17 phase separating into colloid-rich and colloid-lean phases and potentially gelling.18,19 TPT1 qualitatively predicts this transition, but we have found that it underestimates the amount of linker required for phase separation.16,19 We hypothesized that this underestimation was due, in part, to TPT1’s neglect of other prevalent bonding motifs such as linkers that “loop” so that both ends associate to the same colloid [Fig. 1(a)], i.e., form a double bond between a linker and a colloid.

FIG. 1.

Snapshots of associating colloids (red) and linkers (blue), rendered using OVITO 2.9.0.20 Free bonding sites on the colloids are shown in gray. (a) Two linkers “loop” to form a double bond with a colloid, while the other two linkers bridge between different colloids. (b) Two colloids and two linkers form a “ring.” Loop and ring bonding structures are both neglected in TPT1.

FIG. 1.

Snapshots of associating colloids (red) and linkers (blue), rendered using OVITO 2.9.0.20 Free bonding sites on the colloids are shown in gray. (a) Two linkers “loop” to form a double bond with a colloid, while the other two linkers bridge between different colloids. (b) Two colloids and two linkers form a “ring.” Loop and ring bonding structures are both neglected in TPT1.

Close modal

TPT1 has been extended to include double bonds for restricted numbers of components, bonding sites, and/or pairs of sites that can participate in the bonds.21–28 The most general result to date is the work of Marshall,29 which includes multiple components with multiple bonding sites but permits only one pair of sites per component to double bond. This restriction is well justified when the geometry of the sites limits potential double bonding pairs, e.g., direct association of colloids with small numbers of attractive patches. However, it can break down for common linking schemes for self-assembly. For example, nanoparticles are often coated with flexible ligands that are functionalized to bond with a complementary functional group on another nanoparticle30–35 or a linking molecule.19,36–38 Control over the specific number of functionalized ligands on each nanoparticle is infeasible or produces low yields,39 and due to their flexibility and surface coverage, multiple ligands on the same nanoparticle might readily double bond with another nanoparticle or with a linker.19 Given the success of TPT1 for other associating fluids, we aim at developing the theory to handle such cases.

Building on previous TPTs, we derive a general expression for the association free energy of a mixture when a double bond can form between any two pairs of sites on any of the components; this expression reduces to TPT114 and Marshall’s result29 as special cases. We apply the theory to the colloid–linker mixture we previously studied, where the linker readily double bonded with nearly any pair of sites on the colloids.16 We compute the fraction of linkers forming double-bond loops. In the strong association limit, this fraction is strongly dependent on the distribution of linker end-to-end distances. Modifications to linker length or flexibility that make the linker end-to-end distance incompatible with the distance between colloid bonding sites inhibit loop formation; this consideration may prove helpful in the design of molecular linkers for self-assembly. We also show that the theory readily incorporates other prevalent bonding motifs with two-site association, such as “rings” made from multiple linkages between colloids [Fig. 1(b)], and that incorporating these motifs improves its accuracy.

The rest of this article is organized as follows. In Sec. II, we review the fundamentals of Wertheim’s theory, which we hope will serve as a pedagogical reference complementary to and expanding on excellent recent reviews.40,41 In Sec. III, we derive the new theory for double bonding between multiple components with multiple sites. We apply the theory to the colloid–linker mixture in Sec. IV before summarizing in Sec. V.

We begin by outlining Wertheim’s theory for associating fluids,1–4 which is based on a graphical expansion in the grand canonical ensemble using the densities of various bonding states. Our aim is to provide a useful entry point for those seeking to understand and build upon this work, as well as for those who wish to similarly apply or extend Wertheim’s theory. We will follow Wertheim’s derivation for molecules with multiple association sites3 with small modifications to include multiple components; any subscripts or superscripts denoting the component type in this article should be neglected when comparing to the notation from his work.

Consider a mixture of k components that are able to associate with each other through short-ranged attractions. We will adopt a physical model for the mixture that refers to the constituents of a component as “molecules” that are covered with bonding “sites,” but here, a molecule should be interpreted liberally to include not only small molecules but also macromolecules and colloids. A molecule of component i has a set of sites Γi={A1,,Ani}, which we will assume for now are rigidly fixed within the molecule so that its conformation can be specified by its position r and orientation Ω. The total potential energy of the mixture U is taken to be the sum of a one-body potential u1(i) acting on each molecule of component i and a pairwise interaction u2(ij) between two molecules of components i and j. The pairwise interaction is split into a repulsive core uR(ij) and a short-ranged attraction uAB(ij) between sites,

u2(ij)(1,2)=uR(ij)(1,2)+AΓiBΓjuAB(ij)(1,2),
(1)

where 1 = (r1, Ω1) is a shorthand notation for the position and orientation of molecule 1. The site attractions are functions of 1 and 2 through the positions of sites A and B, rA(i)(1) and rB(j)(2). The summation is taken over all pairs of sites in the two molecules, but some sites may be noninteracting.

In the grand canonical ensemble, the temperature T, volume V, and the chemical potentials {μi} of all components are held constant, and the numbers of molecules {Ni} fluctuate. The grand partition function is

Ξ=N1=0Nk=0i=1kdNiΛiNieβμiNiNi!eβU({Ni}),
(2)

where Ni indicates all the positions and orientations of the Ni molecules of component i, Λi accounts for integrals over translational and rotational momenta (i.e., is related to the thermal wavelength), Ni! is a factor for permutating labels of the molecules, and β = 1/(kBT), with kB being Boltzmann’s constant. The one-body contributions u1(i) to the potential energy U can be grouped with the chemical potential μi to define a modified fugacity

zi(1)=Λiexpβ(μiu1(i)(1))
(3)

so that for our model,

Ξ=N1=0Nk=0i=1kdNi1Ni!1zi(1)×expβ1,2u2(ij)(1,2).
(4)

The product of zi is taken over all molecules of component i (indexed 1), and the sum of u2(ij) is taken over all unique pairs of molecules (indexed 1 and 2).

We now define e(ij)(1,2)=exp[βu2(ij)(1,2)] as the e-function (Boltzmann weight) associated with an interaction between two molecules of components i and j, respectively, and f(ij)(1, 2) = e(ij)(1, 2) − 1 as the corresponding Mayer f-function. The integrand in Eq. (4) can be written as a product of f-functions using

expβ1,2u2(ij)(1,2)=1,21+f(ij)(1,2).
(5)

The integrals in Eq. (4) can then be expanded into multiple integrals over different products of f-functions, and it is useful to express these integrals as graphs (Fig. 2).42 Each integral is represented by a “z-graph” whose nodes, called field points, are the integration variables and are each assigned a factor of the fugacity zi. Two field points in a graph may be connected by at most one f-bond whose presence is represented by an edge between the points that carries a weight f(ij). There is an infinite number of z-graphs in Ξ of which only a subset are connected graphs, i.e., all points have joined to each other by f-bonds. It can be shown that ln Ξ is the sum of only these connected graphs, and the mixture thermodynamics can be computed using the grand potential F = −β−1 ln Ξ. However, using only f-bonds to represent association can give a poorly convergent summation for ln Ξ, which motivated Wertheim to further decompose the graphs.1 

FIG. 2.

Graphical representations of the integral ∫d1d2d3z1(1)z2(2)z2(3)f(12)(1, 2)f(12)(1, 3). The first graph is a connected z-graph with f-bonds between molecule 1 of component 1 (red) and molecules 2 and 3 of component 2 (blue). The other graphs show the expansion of f using Eq. (6) into hyperpoints (larger open circles) and bonding sites (smaller filled circles). Hyperpoints connected by an fAB-bond also have an eR-bond (dashed line) between them, while those without an fAB-bond have an fR-bond (dotted line). In this example, a molecule of component 1 has two bonding sites, and a molecule of component 2 has one bonding site. We omit the graphs having two fAB-bonds at a site because these can usually be neglected due to steric hindrance, and we do not show the label-permutation prefactor (1/2) for the z-graph from Eq. (4).

FIG. 2.

Graphical representations of the integral ∫d1d2d3z1(1)z2(2)z2(3)f(12)(1, 2)f(12)(1, 3). The first graph is a connected z-graph with f-bonds between molecule 1 of component 1 (red) and molecules 2 and 3 of component 2 (blue). The other graphs show the expansion of f using Eq. (6) into hyperpoints (larger open circles) and bonding sites (smaller filled circles). Hyperpoints connected by an fAB-bond also have an eR-bond (dashed line) between them, while those without an fAB-bond have an fR-bond (dotted line). In this example, a molecule of component 1 has two bonding sites, and a molecule of component 2 has one bonding site. We omit the graphs having two fAB-bonds at a site because these can usually be neglected due to steric hindrance, and we do not show the label-permutation prefactor (1/2) for the z-graph from Eq. (4).

Close modal

For the model interactions defined by Eq. (1), an f-bond is separable into repulsive and attractive parts,

f(ij)(1,2)=fR(ij)(1,2)+eR(ij)(1,2)AΓiBΓj1+fAB(ij)(1,2)1,
(6)

where the subscripts denote the part of the potential. Graphically, each point in a z-graph is replaced by a hyperpoint encompassing the molecule’s sites, and the f-bonds are replaced by repulsive bonds between hyperpoints and attractive bonds between sites (Fig. 2). Each hyperpoint still carries a factor zi, and repulsive fR-bonds are edges between two hyperpoints. Attractive fAB-bonds are edges between sites A and B in different hyperpoints, which have in parallel a repulsive eR-bond between the hyperpoints.

Two hyperpoints are “directly connected” if they have a repulsive fR-bond or any attractive fAB-bond between their sites. Furthermore, two sites are “bond-connected” if there is a path of attractive bonds between them, while two sites are “constraint-connected” if they are in the same hyperpoint. Now, consider the set of z-graphs where all hyperpoints are connected by networks of attractive bonds; Wertheim called such a graph having s hyperpoints an “s-mer” analogous to the association of monomers into a polymer. The connected z-graphs in ln Ξ can all be expressed as s-mers with fR-bonds between hyperpoints in distinct s-mers, including the s = 1-mer or single hyperpoint without any attractive bonds.

Within an s-mer graph, there can be multiple networks of bond-connected sites; these networks must be constraint-connected. Consider two hyperpoints 1 and 2 in an s-mer that are not directly connected by an attractive bond. For such a hyperpoint pair, Wertheim proposed that if any site A on 1 is bond-connected to any site B on 2, two s-mer graphs with and without an fR-bond between 1 and 2 can be usefully combined into one graph with an eR-bond between 1 and 2. Hence, we will form eR-bonds between all pairs of hyperpoints in an s-mer with bond-connected sites, but pairs of hyperpoints in an s-mer that are connected only through a constraint connection will not have an eR-bond. By this procedure, all hyperpoints in a bond-connected network become irreducibly connected.43 

Define ρ(i)(1) to be the sum of z-graphs obtained by turning a field hyperpoint of component i, which is an integration variable, into a “root” hyperpoint labeled 1, which is not an integration variable. This is equivalent to the functional differentiation,

ρ(i)(1)=zi(1)δlnΞδzi(1),
(7)

so ρ(i)(1) is the average “singlet” density of component i at 1 in the grand canonical ensemble. The singlet density is related to the familiar number density ρi(r1) by integration over orientations,

ρi(r1)=dΩ1ρ(i)(1).
(8)

The various z-graphs in the sum ρ(i)(1) can be further analyzed in terms of their bonding at the root hyperpoint. A z-graph is assigned to the sum ρα(i)(1) if α is the set of bonded sites at 1 (α ⊆ Γi),44 and as a result,

ρ(i)(1)=αΓiρα(i)(1).
(9)

Now, let c(i)(1) be the sum of the subset of graphs in ρ(i)(1)/zi(1) for which 1 is not an articulation point. It can be shown that41 

lnρ(i)(1)zi(1)=c(i)(1).
(10)

Furthermore, let cα(i)(1) be the sum of the subset of graphs in ρα(i)(1)/ρ(i)(1) for which 1 is not a “constraint” articulation point, i.e., the removal of the hyperpoint and any incident repulsive bonds at 1 leaves the graph connected. ρα(i)(1) includes all graphs bonded at the sites α, which can be formed by composing at the root point various ca(i)(1) for aα, such that their union gives α, along with the graphs that are not bonded,

ρα(i)(1)=ρ(i)(1)γP(α)aγca(i)(1).
(11)

The sum is over all set partitions of α, denoted as P(α),45 and the product is over all sets in each partition.

As currently defined, cα(i) is a sum of both reducible and irreducible z-graphs, and it is computationally desirable to eliminate the reducible ones. To do so, Wertheim considered the procedure of adding at each field point of an irreducible graph either nothing or another graph with one root point. By the eR-bond replacement procedure for an s-mer, the root point of the added graph must have bonding complementary to the bonding at the field point of the irreducible graph, so the added graphs can be any ργ(i) for γ ⊆ Γi\α.46 This naturally leads to the definition of the auxiliary singlet densities,

σα(i)(1)=γαργ(i)(1),
(12)

which are the sums of all graphs whose bonded sites at root point 1 are a subset of α. For example, σA(i) contains all graphs bonded at A (and not bonded elsewhere) or having no bonded sites.47 Two special cases are σ(i)=ρ(i), the sum of graphs not bonded at any site, and σΓi(i)=ρ(i), including all possible combinations of sites. The bonded singlet densities ρα(i) can be expressed in terms of σα(i),48 

ρα(i)(1)=γα(1)|α\γ|σγ(i)(1),
(13)

with details in  Appendix A.

Each field point of an irreducible graph that has a set of bonded sites α is now reassigned a factor σΓi\α(i), which are precisely the graphs not bonded at α, and cα(i) can be expressed as a sum of only irreducible graphs. We define c to be the sum of all irreducible graphs consisting of s-mer graphs and fR-bonds between s-mers where all hyperpoints are field points and carry these factors. (Wertheim called this sum c(0), but we will use c to avoid confusion with the type superscripts.) All cα(i) can be obtained from c by turning a field point bonded at α into a root point labeled 1 and removing its factor σΓi\α(i). This procedure is given by the functional differentiation,

cα(i)(1)=δcδσΓi\α(i)(1).
(14)

cα(i) can also be expressed as functions of σα(i) directly,

cα(i)=δ|α|,1γP(α)(1)|γ|(|γ|1)!aγσa(i)σ(i)
(15)

for α, with δi,j being the Kronecker delta. [Note that this corrects a typographical error in Eq. (27) of Ref. 3.] The inverse relationship is

σα(i)=σ(i)γP(α)aγδ|a|,1+ca(i).
(16)

See  Appendix A for details required to arrive at these results.

Given σα(i) and cα(i), it can be shown using a Legendre transform of the grand potential that the Helmholtz free energy A of the associating mixture is41 

βA=i=1kd1σΓi(i)(1)lnΛi1σ(i)(1)σΓi(i)(1)+βu1(i)(1)σΓi(i)(1)+αΓiασΓi\α(i)(1)cα(i)(1)c.
(17)

The integrand can be expressed in terms of only σα(i) using Eq. (15) to replace cα(i) (see  Appendix A),

αΓiασΓi\α(i)(1)cα(i)(1)=σΓi(i)(1)+Q(i)(1),
(18)

with

Q(i)=AΓiσΓi\A(i)+σ(i)γP(Γi)γ{Γi}(1)|γ|(|γ|2)!aγσa(i)σ(i),
(19)

to give

βA=i=1kd1σΓi(i)(1)lnΛi1σ(i)(1)+βu1(i)(1)σΓi(i)(1)+Q(i)(1)c.
(20)

At equilibrium, A should be minimized with respect to variations in σα(i), subject to a constraint on σΓi(i) to give a fixed number of molecules. Taking the functional derivative of A with respect to σα(i) for α ⊂ Γi gives49 

0=δ(βA)δσα(i)(1)=δ|α|,0σΓi(i)σ(i)+Q(i)σα(i)cΓi\α(i).
(21)

That this relationship holds can be verified by taking derivatives of Eq. (19) and comparing them to Eq. (15). Functional differentiation of A with respect to the total singlet density σΓi(i) gives the chemical potential μi after some manipulation, which is the correct thermodynamic relationship.

The results so far are formally exact, but Wertheim’s theory has been most successfully applied as an approximate TPT.2,4 In TPT, the free energy is computed relative to a reference system with molecules whose pairwise interactions are given by only uR(ij) and whose free energy AR is more readily computed. The reference system has no attractive bonds and is assumed to have the same total densities σΓi(i),

βAR=i=1kd1σΓi(i)(1)lnΛi1σΓi(i)(1)σΓi(i)(1)+βu1(i)(1)σΓi(i)(1)cR,
(22)

where cR is the sum of irreducible graphs for the reference fluid without any bonding, i.e., it contains only fR-bonds and should only be a functional of the total density σΓi(i). By identifying the first term as the ideal free energy, −cR is the excess contribution to βAR. In practice, AR is usually evaluated with a known equation of state for the reference system; for example, AR for a spatially homogeneous hard-sphere mixture is well approximated by Boublík’s equation of state.50 

The association free energy, ΔA = AAR, is then

βΔA=i=1kd1σΓi(i)(1)lnσ(i)(1)σΓi(i)(1)+σΓi(i)(1)+Q(i)(1)Δc.
(23)

The last term Δc = ccR is the key to Wertheim’s theory, as it is the contributions to the sum of irreducible graphs from association. Although Δc formally contains all irreducible graphs with at least one attractive bond, Wertheim discussed a few cases of “steric hindrance” that can simplify the graphs that need to be included. Two particularly convenient ones are as follows:

  1. The formation of an fAB-bond between site A on molecule 1 and site B on molecule 2 prevents any site in third molecule 3 from bonding with either A or B. This hindrance usually occurs because of repulsion between the cores of the molecules.

  2. Two sites B and C on molecule 2 cannot bond to the same site A on molecule 1. This can also be enforced by the site geometry or by restricting the model to form exclusive bonds (e.g., covalent chemistry).

Together, these hindrances remove all graphs with multiple bonding of a single site, which greatly reduces the combinations that must be considered. Further simplification can be made by including only a limited subset of irreducible graphs, which we do next.

TPT1 usually includes the lowest-order irreducible graphs from pair association in Δc.2,4 This admits at most one attractive bond between two molecules and enforces a tree-like bonded network. To extend TPT1 to include double bonding, we also include the irreducible graphs with two attractive bonds between two molecules,

Δc12i=1kAΓij=1kBΓjd1d2×σΓi\A(i)(1)σΓj\B(j)(2)gR(ij)(1,2)fAB(ij)(1,2)+12i=1kABΓi|AB|=2j=1kCDΓj|CD|=2d1d2×σΓi\AB(i)(1)σΓj\CD(j)(2)gR(ij)(1,2)fABCD(ij)(1,2),
(24)

where gR(ij) is the pair correlation function in the reference fluid. The sums over AB and CD are taken over the two-site pairs on i and j (e.g., AB = {A, B} is a two-element subset of Γi), and

fABCD(ij)=fAC(ij)fBD(ij)+fAD(ij)fBC(ij)
(25)

accounts for the two way sites in AB, and CD can associate with each other. The first term in Eq. (24) is standard in TPT1 and accounts for all ways a single bond can form between two molecules, while the second term accounts for all ways a double bond can form. In practice, some or all these graphs may be zero or negligible for a given model due to the bonding chemistry and/or site geometry, but we do not exclude any of these at this stage.

In expressing the graph sum in this way, we have neglected any rings (cycles) in the graph sum beyond the double bond, although all such graphs are irreducible.22,51 Hence, the overall bonding network is still a tree of single and double bonds between molecules. Rings of molecules bonded at two sites can be included within the framework we develop here using the approximation of Sear and Jackson;22 we will revisit this point in Sec. IV. We also neglect the effects of steric hindrance that prevent bonding at one site when a bond forms at another; these effects might be treated by including additional graphs or with a renormalization approach.27,28,52

We now differentiate Eq. (24) with respect to σα(i), giving

ΔcA(i)(1)=j=1kBΓjd2σΓj\B(j)(2)gR(ij)(1,2)fAB(ij)(1,2)
(26)

and

ΔcAB(i)=j=1kCDΓj|CD|=2d2σΓj\CD(j)(2)gR(ij)(1,2)fABCD(ij)(1,2),
(27)

so

Δc=12i=1kAΓid1σΓi\A(i)(1)ΔcA(i)(1)+12i=1kABΓi|AB|=2d1σΓi\AB(i)(1)ΔcAB(i)(1).
(28)

The remaining Δcα(i) for α ⊂Γi are zero. Because cR is only a functional of σΓi(i), we can replace cα(i) by Δcα(i) in Eq. (18),

σΓi(i)+Q(i)=AΓiσΓi\A(i)ΔcA(i)+ABΓi|AB|=2σΓi\AB(i)ΔcAB(i).
(29)

We may further express ΔcA(i) and ΔcAB(i) in terms of σα(i) using Eq. (15),

ΔcA(i)=1+σA(i)σ(i)
(30)

and

ΔcAB(i)=σAB(i)σ(i)σA(i)σB(i)(σ(i))2.
(31)

These results can be directly substituted into Eq. (23) to define ΔA using only σα(i); however, this system of equations is still incomplete for specifying all σα(i) (and hence, ΔA).

We use Eq. (16) to find a relationship between σα(i) in terms of ΔcA(i) and ΔcAB(i), given that all other Δcα(i) are zero for α ⊂ Γi,

σα(i)=σ(i)γP(α)|a|2aγAγ|A|=11+ΔcA(i)ABγ|AB|=2ΔcAB(i).
(32)

The sum is taken over all partitions of α having only one-element or two-element subsets, and the products are taken over each type of subset; all other partitions do not contribute. For |α|≥ 3, Eq. (32) can be restated as a recursive relationship,

σα(i)=σA(i)σα\A(i)σ(i)+Bα\AσAB(i)σ(i)σA(i)σB(i)(σ(i))2σα\AB(i),
(33)

after the substitution of Eqs. (30) and (31). The first term accounts for A being a one-element subset in the first product of Eq. (32), which contributes a factor of 1+ΔcA(i); this factor then multiplies σα\A(i), which includes all partitions of α\A. The second term accounts for all two-element subsets containing A, which contribute a factor of ΔcAB(i) in the second product of Eq. (32) for all Bα\A. This must multiply σα\AB(i), which includes all remaining partitions of α\AB. For example, if α = ABC = {A, B, C} and we remove A, we obtain

σABC(i)=σA(i)σBC(i)σ(i)+σAB(i)σ(i)σA(i)σB(i)(σ(i))2σC(i)+σAC(i)σ(i)σA(i)σC(i)(σ(i))2σB(i)=σA(i)σBC(i)σ(i)+σB(i)σAC(i)σ(i)+σC(i)σAB(i)σ(i)2σA(i)σB(i)σC(i)(σ(i))2.
(34)

As expected, the same expression is obtained if site B or C is initially removed from α rather than A. The partitions of α having only one-element or two-element subsets are {{A}, {B, C}}, {{B}, {A, C}}, {{C}, {A, B}}, and {{A}, {B}, {C}}, so Eq. (32) is also equivalent after substitution and some simplification.

We now recast these equations using the notation of Chapman and co-workers11–14 and define Xα(i)(1)=σΓi\α(i)(1)/σΓi(i)(1) as the local fraction of molecules of type i not bonded at all sites in α. [Note the special case X(i)(1)=1.] The association free energy is

βΔA=i=1kd1σΓi(i)(1)lnXΓi(i)(1)+12AΓiXA(i)(1)×XΓi\A(i)(1)XΓi(i)(1)1+12ABΓi|AB|=2XAB(i)(1)×XΓi\AB(i)(1)XΓi(i)(1)XΓi\A(i)(1)XΓi\B(i)(1)(XΓi(i)(1))2,
(35)

while the “mass action” or “chemical equilibrium” conditions from Δcα(i) are

1+XΓi\A(i)(1)XΓi(i)(1)=j=1kBΓjd2σΓj(j)(2)XB(j)(2)gR(ij)(1,2)fAB(ij)(1,2)
(36)

and

XΓi\AB(i)(1)XΓi(i)(1)XΓi\A(i)(1)XΓi\B(i)(1)(XΓi(i)(1))2   =j=1kCDΓj|CD|=2d2σΓj(j)(2)XCD(j)(2)gR(ij)(1,2)fABCD(ij)(1,2).
(37)

If the mixture is spatially homogeneous, none of the singlet densities (or Xα(i)) depends on position or orientation. Carrying out the integral in Eq. (35) gives the association free-energy density Δa = ΔA/V in terms of the number densities ρi and Xα(i),

βΔa=i=1kρilnXΓi(i)+12AΓiXA(i)XΓi\A(i)XΓi(i)1+12ABΓi|AB|=2XAB(i)XΓi\AB(i)XΓi(i)XΓi\A(i)XΓi\B(i)(XΓi(i))2,
(38)

integrating Eq. (36) over 1 gives

1+XΓi\A(i)XΓi(i)=j=1kρjBΓjXB(j)ΔAB(ij),
(39)

with

ΔAB(ij)=(ΩiΩj)1dr12dΩ1dΩ2gR(ij)(r12,Ω1,Ω2)×fAB(ij)(r12,Ω1,Ω2),
(40)

and integrating Eq. (37) over 1 gives

XΓi\AB(i)XΓi(i)XΓi\A(i)XΓi\B(i)(XΓi(i))2=j=1kρjCDΓj|CD|=2XCD(j)ΔABCD(ij),
(41)

with

ΔABCD(ij)=(ΩiΩj)1dr12dΩ1dΩ2gR(ij)(r12,Ω1,Ω2)×fABCD(ij)(r12,Ω1,Ω2),
(42)

where r12 = r2r1 is the separation between molecules and Ωi is the integral over orientations of a molecule of component i. Here, we have used that gR(ij) should only depend on the relative separation between molecules in a spatially homogeneous mixture. ΔAB(ij) and ΔABCD(ij) are then “bond volumes” integrated over separations and averaged over all orientations of the two molecules.13,14

Equation (32) can also be restated as

XΓi\α(i)=XΓi(i)γP(α)|a|2aγAγ|A|=1XΓi\A(i)XΓi(i)×ABγ|AB|=2XΓi\AB(i)XΓi(i)XΓi\A(i)XΓi\B(i)(XΓi(i))2
(43)

and equivalently, as a recursive relationship,

XΓi\α(i)=XΓi\A(i)XΓi\(α\A)(i)XΓi(i)+Bα\AXΓi\AB(i)XΓi(i)XΓi\A(i)XΓi\B(i)(XΓi(i))2XΓi\(α\AB)(i).
(44)

These relationships allow XA(i) and XAB(i) to be computed in terms of XΓi(i) and all XΓi\A(i) and XΓi\AB(i), and they also self-consistently define XΓi(i). Taken together, Eqs. (39), (41), and (43) comprise a system of nonlinear equations that can be numerically solved to obtain all Xα(i), with Eq. (38) giving the free energy; in practice, this may be challenging for certain parameter values. There are, however, two cases that simplify considerably: (1) no double bonds can form and (2) a double bond can form at only one pair of sites on each molecule. We now consider these, in turn, before applying the general theory to a colloid–linker mixture in Sec. IV.

In TPT1, at most one bond can form between two molecules,2,4 which is a reasonable approximation when the site geometry on the molecules inhibits double bonding, e.g., patchy colloids with two bonding sites on opposite hemispheres, and the formation of bonded rings of molecules is unlikely, e.g., because the bonding sites are rigidly fixed. (This amounts to setting ΔABCD(ij)=0 for all pairs of sites on all molecules.) TPT1 is the foundation of SAFT,11–15 so we will show that we recover the SAFT equations originally derived by Chapman for multicomponent mixtures with multiple bonding sites14 if we neglect double bonding between molecules.

Making the simplification that ΔABCD(ij)=0, Eq. (44) becomes

XΓi\α(i)=XΓi\A(i)XΓi\(α\A)(i)XΓi(i)
(45)

for any Aα, which when evaluated for α = Γi gives

1=XΓi\A(i)XA(i)XΓi(i),
(46)

and it follows that

XΓi(i)=AΓiXA(i).
(47)

This relationship captures the statistical independence of bonding in TPT1: the fraction of molecules of component i that are not bonded at any site is the product of the fraction not bonded at each site. The free energy

βΔa=i=1kρiAΓilnXA(i)+121XA(i)
(48)

and chemical equilibrium equations

XA(i)=1+j=1kρjBΓjXB(j)ΔAB(ij)1
(49)

follow by substitution and are equivalent to Eqs. (3) and (4) of Ref. 14 after a conversion to express ΔA per molecule.

As another special case, suppose that exactly two sites A12 = {A1, A2} on a molecule of component i can form a double bond with a similar set of two sites on a molecule of component j. (We will label all these sites A1 and A2, but A1 on component i and A1 on component j need not be the same, and similarly for A2.) This restriction may physically correspond to a molecular geometry where only two of the sites are close enough to double bond with another molecule, which was studied by Marshall and Chapman for a one-component fluid28 and by Marshall for a multicomponent mixture29 or a restrictive bonding chemistry.

The immediate consequence is that only ΔA12A12(ij)0, which drastically simplifies the expressions for σα(i) because a maximum of two partitions of α contribute to Eq. (43)—the partition into all singleton sets and the partition containing one two-element subset A12—so

XA(i)=XΓi(i)XΓi\A(i),AA12XA12(i)XΓi\(A12\A)(i)XΓi(i),AA12.
(50)

The first case, when A is not one of the double-bond sites, is obtained from Eq. (44) with α = Γi and is the same as in TPT1. The second case, when A is one of the double-bond sites, is obtained from Eq. (44) with α = Γi\A1 and Γi\A2 and removing A2 and A1, respectively. Combining Eq. (50) with Eq. (44) for α = Γi and either AA12 removed gives the additional relationship

XA12(i)XΓi\A12(i)XΓi(i)XΓi\A1(i)XΓi\A2(i)(XΓi(i))2=1XA1(i)XA2(i)XA12(i).
(51)

Finally, it follows from applying Eq. (44) recursively that

XΓi(i)=XA12(i)AΓi\A12XA(i).
(52)

With these results, the free energy is

βΔa=i=1kρiAΓilnXA(i)+121XA(i)+lnXA12(i)XA1(i)XA2(i)121XA1(i)XA2(i)XA12(i).
(53)

The free energy has been written so that the sum over sites A is the same as TPT1, while the additional terms give the free energy due to double-bond formation. This expression is equivalent to Eq. (19) of Ref. 29 after some manipulations. The corresponding chemical equilibrium equations for AA12 are the same as Eq. (49), while for αA12, the three Xα(i) are given by the system of equations,

1+XA2(i)XA12(i)=j=1kρjBΓjXB(j)ΔA1B(ij),
(54)

a similar equation with the labels A1 and A2 reversed, and

1XA12(i)XA1(i)XA2(i)(XA12(i))2=j=1kρjXA12(j)ΔA12A12(ij).
(55)

Having extended Wertheim’s TPT to include general double-bond association, we now apply the theory to a model mixture of colloids and flexible linker molecules that we previously studied with molecular simulations and TPT1.16 In that work, we aimed at predicting the phase behavior of the homogeneous mixture; however, TPT also gives valuable information about the prevalence of different bonding motifs. We measured a significant fraction of linkers forming loops (a double bond with both linker ends attached to the same colloid) in our simulations. Linker loops are neglected in TPT1, so here we investigate the linker loop fraction under different conditions as a useful test of the new theory.

The linkers we studied were linear chains of M = 8 tangentially bonded hard-sphere segments with diameter dl, which we also called “polymers” in Ref. 16, while the colloids were hard spheres with diameter dc = 5 dl. The mixture composition was defined using the colloid volume fraction ηc=ρcπdc3/6, which we varied from 0.01 to 0.15, and the linker-to-colloid number ratio ρl/ρc, which we fixed at 1.5. Each linker had nl = 2 bonding sites, one at the center of each end segment. Each colloid had nc = 6 bonding sites arranged at the vertices of an octahedron at a distance dcl*=dcl+(21/61)dl3.12dl from the center of the colloid, where dcl = (dc + dl)/2 = 3 dl is the hard-sphere contact distance for a colloid and a linker segment. (Placement of the colloid bonding sites at dcl* rather than dcl was due to our original work using nearly hard potentials for the spheres.53) A linker could reversibly bond to a colloid through a short-ranged site attraction,

uAB(cl)(1,2)=εe(rAB/)2,rAB<2.50,rAB2.5
(56)

where rAB=|rB(l)(2)rA(c)(1)| is the distance between the linker and colloid sites, ε is the strength of the attraction, and = 0.2 dl sets the range of the attraction. The potential was truncated for rAB ≥ 2.5 = 0.5 dl. There was at most one bond at each site, and there were no bonds between two colloid sites or two linker sites.

We developed the theory in Secs. II and III with all sites fixed on rigid molecules, but in our colloid–linker model, the linker bonding sites fluctuate with the end-to-end vector R. This intramolecular flexibility can be incorporated into the theory with a few modifications and approximations.21,22,52 To specify the linker’s conformation, we redefined 1 = (r1, R1), where r1 is the position of one end of the chain and R1 is its end-to-end vector, effectively coarse-graining the internal linker segments. To appropriately weight the conformations, the linkers were given an intramolecular potential

βu1(l)(R)lnp(R),
(57)

where p(R) is the probability distribution of end-to-end vectors in the reference mixture without bonding. The linker singlet densities now depend on R, and this dependence may be different when both ends of a linker are bonded than when they are free. As a result, σΓl(l) may have a different R-dependence in the associating mixture than in the reference mixture, and computing AR at the same σΓl(l) may be inconvenient for equations of state that are functionals of the number density. Wertheim showed how this can be circumvented in certain cases by choosing a reference mixture having the same number density after integration over internal degrees of freedom and making a small modification to Eq. (23).52 We accordingly take the reference state to be the hard-chain mixture at the same number density and without bonding (βε = 0).

With these modifications, we specialized the graph sum given by Eq. (24) to our model. We first reduced the number of σα(c) that needed to be explicitly included. For a linker to double bond, the arc-length distance between colloid bonding sites must be less than the linker’s contour length L = (M − 1)dl = 7 dl. In our model, the arc-length distance between nearest sites on the colloid was dcl*π/24.9dl, so the linker was able to form double bonds with any of these pairs of sites. However, the arc-length distance between sites on opposite hemispheres was dcl*π9.8dl, so a linker could not form a double bond with these sites. We accordingly excluded double bonds between the three pairs of sites on opposite sides of the colloid; the remaining νc = nc(nc − 2)/2 double-bond pairs were all equivalent. We then defined the graph sum using representative sets of colloid bonding sites, A1 for a single bond and A12 = {A1, A2} for a double bond between nearest sites. The linker had nl = 2 equivalent bonding sites B12 = {B1, B2}.

We next simplified the chemical equilibrium equations for a spatially homogeneous mixture. The linker singlet densities were constant with respect to r but not R, so we redefined Xα(l)(r1)=σΓl\α(l)(r1)/σΓl(l)(r1) for the linkers using singlet densities integrated over R1. We further assumed that the distribution of end-to-end vectors was the same for linkers having neither or one end bonded and with both given by p(R), i.e., ρ(l)(1)=ρ(l)(r1)p(R1) and ρB1(l)(1)=ρB1(l)(r1)p(R1). Carrying through the integration and explicitly specializing for the difunctional linker gave

1+XΓc\A1(c)XΓc(c)=2ρlXB1(l)Δ1,
(58)
1+XB1(l)XB12(l)=ncρcXA1(c)Δ1,
(59)

and

XΓc\A12(c)XΓc(c)XΓc\A1(c)XΓc(c)2=2ρlXB12(l)Δ2,
(60)
1XB12(l)XB1(l)XB12(l)2=2νcρcXA12(c)Δ2,
(61)

with single-bond volume

Δ1Ωc1dr12dΩ1dR2p(R2)gR(cl)(r12,R2)fA1B1(cl)(r12,Ω1)
(62)

and double-bond volume

Δ2Ωc1dr12dΩ1dR2p(R2)gR(cl)(r12,R2)×fA1B1(cl)(r12,Ω1)fA2B2(cl)(r12,Ω1,R2).
(63)

The factors of 2 in Eqs. (60) and (61) account for equivalent permutations of the site labels. Additionally, Eq. (43) gives XA1(c) and XA12(c) as functions of XΓc\A12(c), XΓc\A1(c), and XΓc(c), along with an implicit equation for XΓc(c) in these variables. Together, these make a complete set of nonlinear equations in five variables that can be solved for a given mixture composition and attraction strength ε after the evaluation of Δ1 and Δ2. Complete details of how we evaluated these integrals and solved these equations are given in  Appendix B.

The pair correlation function between the colloids and linkers in the reference system gR(cl)(r12,R2) depends not only on the distance between the colloid and the linker r12 but also on the linker’s conformation R2. We considered approximating gR(cl) using a superposition of the hard-sphere correlations ghs(cl) in the mixture obtained by dissolving the linker’s internal bonds, but we suspected this approach might overestimate gR(cl) because segments of polymer chains are known to be depleted near the surface of colloids,54 whereas hard spheres are enriched.50 

We accordingly measured gR(cl) in molecular dynamics simulations of the reference hard-chain mixture using the model and simulation methods of Ref. 16. We simulated 1000 colloids and 1500 linkers in a cubic box with periodic boundary conditions using LAMMPS (22 Aug 2018).53,55,56 Starting at ηc = 0.01, we first equilibrated the mixture for 1.5 × 104τ, where τ is the unit of time in the simulations. We then sampled configurations for analysis every 10 τ over a 104τ period (1000 configurations). Finally, we linearly compressed the edge length of the simulation box over 5 × 103τ to reach ηc = 0.02. We repeated this procedure up to and including ηc = 0.15. We computed gR(cl) as a histogram of three variables: the distance from the center of the colloid to one end of the linker rB1=|r12|, the end-to-end distance R = |R2|, and the polar angle ϕ between R2 and r12, i.e., cosϕ=(r12R2)/(rB1R2). The histogram bin ranges and widths were 3dlrB14dl with width 0.25 dl, 0 dlR ≤ 7 dl with width 0.5 dl, and 0 ≤ ϕπ with width π/9 (20°). We have labeled rB1 using the end B1, but we also took the other linker end B2 as r12 to improve sampling.

We also computed the end-to-end vector distribution p(R) as a histogram of R using the same R-bins as for gR(cl) (Fig. 3). This distribution was largely independent of ηc at the compositions we simulated, so we used p(R) at ηc = 0.01 to calculate Δ1 and Δ2. [We used the p(R) measured at each composition, however, to normalize gR(cl).] We noted that the probability of having R commensurate with the distance between two colloid bonding sites (2dcl*4.4dl) was nonnegligible, meaning that it was likely that double bonds would form.

FIG. 3.

Probability density function p for the linker end-to-end vector R as a function of the end-to-end distance |R| at ηc = 0.01 (red circles) and ηc = 0.15 (blue squares). The black dotted line shows an empirical fit to the ηc = 0.01 distribution for |R| ≥ dl, which we used for numerical convenience when evaluating the bond-volume integrals ( Appendix B).

FIG. 3.

Probability density function p for the linker end-to-end vector R as a function of the end-to-end distance |R| at ηc = 0.01 (red circles) and ηc = 0.15 (blue squares). The black dotted line shows an empirical fit to the ηc = 0.01 distribution for |R| ≥ dl, which we used for numerical convenience when evaluating the bond-volume integrals ( Appendix B).

Close modal

The integrand of Δ1 is nonzero only when one end of the linker interacts with a bonding site. Anticipating the hard-sphere exclusion between the colloid and linker segment, this requires dclrB1r*, where r*=dcl*+2.5 based on the range of Eq. (56). We found that gR(cl) did not change significantly in this range of rB1, tending to increase slightly for larger rB1. In the strong association limit of large βε, the dominant contribution to the integral occurs when the bonding sites overlap at rB1dcl*. We accordingly approximated gR(cl)(rB1,R,ϕ)gR(cl)(dcl*,R,ϕ), which still depended on both R and ϕ, e.g., because the other linker end cannot penetrate the colloid. However, further inspection of the data suggested a function of only one other variable, namely, the distance rB2 from the center of the colloid to the other end of the linker,

rB2=|r12+R2|=(rB12+R2+2rB1Rcosϕ)1/2.
(64)

Conformations having rB2<dcl had gR(cl)(dcl*,rB2)=0 because this would cause the linker end to penetrate the colloid. For larger values of rB2, we typically found gR(cl)(dcl*,rB2)<1 for most conformations, indicating depletion of the linker near the surface of the colloid.54 

We noted that gR(cl)(dcl*,rB2) also had a volume fraction dependence. Conveniently, scaling gR(cl) by ghs(cl)(dcl+), which we obtained from Boublík’s equation of state,50 effectively accounted for this composition dependence and collapsed the data [Fig. 4(a)]. [We included in Fig. 4(a) only data for dlR ≤ 5 dl and 0 ≤ ϕπ/2, for which we had sufficient sampling. These conformations are expected to contribute most significantly to Δ1, as others are less probable or wholly excluded.] Empirically, the data were well fit by a linear function of rB2 for rB2>dcl, which we used to evaluate Δ1.

FIG. 4.

Simulated colloid–linker pair distribution function gR(cl) in reference mixtures having varied ηc and fixed ρl/ρc = 1.5. gR(cl) is normalized by the contact value of the pair distribution function ghs(cl)(dcl+) in the hard-sphere mixture that would be obtained by removing all bonds from the linkers.50 (a) When one linker end is fixed at rB1=dcl*, gR(cl)/[ghs(cl)(dcl+)] collapses as a linear function (dashed line) of the distance from the colloid center to the other chain end rB2. Here, we show data only for dlR ≤ 5 dl and 0 ≤ ϕπ/2, for which we have reliable sampling. (b) When both ends of the linker are fixed at colloid bonding sites so that rB1=dcl*, R=2dcl* and ϕ = 3π/4, gR(cl)/[ghs(cl)(dcl+)]2 is a constant (dashed line). The inset of (b) illustrates the definitions of rB1, rB2, R, and ϕ.

FIG. 4.

Simulated colloid–linker pair distribution function gR(cl) in reference mixtures having varied ηc and fixed ρl/ρc = 1.5. gR(cl) is normalized by the contact value of the pair distribution function ghs(cl)(dcl+) in the hard-sphere mixture that would be obtained by removing all bonds from the linkers.50 (a) When one linker end is fixed at rB1=dcl*, gR(cl)/[ghs(cl)(dcl+)] collapses as a linear function (dashed line) of the distance from the colloid center to the other chain end rB2. Here, we show data only for dlR ≤ 5 dl and 0 ≤ ϕπ/2, for which we have reliable sampling. (b) When both ends of the linker are fixed at colloid bonding sites so that rB1=dcl*, R=2dcl* and ϕ = 3π/4, gR(cl)/[ghs(cl)(dcl+)]2 is a constant (dashed line). The inset of (b) illustrates the definitions of rB1, rB2, R, and ϕ.

Close modal

The integrand of Δ2 is nonzero only when both ends of the linker interact with bonding sites. Using similar reasoning as for Δ1, we approximated gR(cl) by its value when all the bonding sites overlap, gR(cl)(rB1,R,ϕ)gR(cl)(dcl*,2dcl*,3π/4). The value of gR(cl) again depended on ηc, but we found that it was essentially a constant across volume fractions when scaled by [ghs(cl)(dcl+)]2 [Fig. 4(b)]. Qualitatively, the additional factor of ghs(cl)(dcl+) here accounts for the presence of the other chain end near the surface.21,22,51 The value of gR(cl) is significantly less than the superposition of the hard-sphere correlations, presumably due to the large entropic penalty to confine the linker near the surface of the colloid.

With all inputs determined, we proceeded to compute the fraction of colloids and linkers in different bonding states at various attraction strengths ε and colloid volume fractions ηc. We first determined the fractions of linkers having no, one, or both ends bonded using TPT; these fractions are XB12(l), 2(XB1(l)XB12(l)), and 12XB1(l)+XB12(l), respectively, using Eq. (13). Figure 5 shows a representative result for ηc = 0.10 as a function of ε, which was a composition where the mixture remained a single (homogeneous) phase in our previous simulations.16 The fraction of unbonded linkers decreased monotonically as the attraction strength ε increased, the fraction of linkers bonded at both ends concomitantly increased monotonically, and the fraction of linkers with only one end bonded had a maximum near βε ≈ 12.5. The TPT calculation (blue dashed line) agreed nearly quantitatively with the simulation data (black circles).16 However, first-order TPT, which completely neglects double bonds between colloids and linkers, yielded comparable results that we omitted here for clarity.

FIG. 5.

Fractions of linkers with (a) no bonded ends, (b) one bonded end, and (c) two bonded ends at ηc = 0.10 and ρl/ρc = 1.5 as a function of attraction strength ε. The lines are TPT predictions accounting for only loops (blue dashed lines) and for both loops and two-colloid rings (red solid lines), and the black circles are simulation data.16.

FIG. 5.

Fractions of linkers with (a) no bonded ends, (b) one bonded end, and (c) two bonded ends at ηc = 0.10 and ρl/ρc = 1.5 as a function of attraction strength ε. The lines are TPT predictions accounting for only loops (blue dashed lines) and for both loops and two-colloid rings (red solid lines), and the black circles are simulation data.16.

Close modal

The usefulness of the newly developed theory is its ability to predict the fraction of linkers in loops χ, which is absent from TPT1. This fraction can be computed from the double-bond graphs in ρB12(l), which here is

χ=d1ρ(l)(1)ΔcB12(l)(1)d1ρB12(l)(1)=1(XB1(l))2XB12(l).
(65)

In TPT1, χ = 0 because ΔcB12(l)=0 (Sec. III A), but with double bonding included, χ should depend on ε and on the mixture composition. We computed χ as a function of ε at volume fractions ηc = 0.01 and ηc = 0.10 (red dotted lines in Fig. 6). At both compositions, more loops formed as ε increased, and there were significantly more loops at the dilute composition ηc = 0.01 than when ηc = 0.10. The latter is more obvious considering χ as a function of ηc in the strong association limit, which we take here as βε = 20 (Fig. 7). Both of these trends in the TPT predictions are qualitatively consistent with our simulations (red circles).

FIG. 6.

Fraction of linkers χ in loops (red circles/lines) and two-colloid rings (blue squares/lines) as a function of attraction strength ε at ρl/ρc = 1.5 and (a) ηc = 0.01 and (b) ηc = 0.10. Simulation data from Ref. 16 (symbols) are compared to two TPT calculations: one including only loops (χ = χloop as red dotted lines), and one including both loops (χloop as red solid lines) and two-colloid rings (χring as blue dashed lines).

FIG. 6.

Fraction of linkers χ in loops (red circles/lines) and two-colloid rings (blue squares/lines) as a function of attraction strength ε at ρl/ρc = 1.5 and (a) ηc = 0.01 and (b) ηc = 0.10. Simulation data from Ref. 16 (symbols) are compared to two TPT calculations: one including only loops (χ = χloop as red dotted lines), and one including both loops (χloop as red solid lines) and two-colloid rings (χring as blue dashed lines).

Close modal
FIG. 7.

Fraction of linkers χ in loops (red circles/lines) and two-colloid rings (blue squares/line) at βε = 20 and ρl/ρc = 1.5 as a function of colloid volume fraction ηc. Simulation data from Ref. 16 (symbols) are compared to two TPT calculations: one including only loops (χ = χloop as red dotted line), and one including both loops (χloop as red solid line) and two-colloid rings (χring as blue dashed line).

FIG. 7.

Fraction of linkers χ in loops (red circles/lines) and two-colloid rings (blue squares/line) at βε = 20 and ρl/ρc = 1.5 as a function of colloid volume fraction ηc. Simulation data from Ref. 16 (symbols) are compared to two TPT calculations: one including only loops (χ = χloop as red dotted line), and one including both loops (χloop as red solid line) and two-colloid rings (χring as blue dashed line).

Close modal

By manipulating the chemical equilibrium equations for the linkers, it can be shown in the strong association limit that

χ1+Kρcnc22νcΔ12Δ21,
(66)

where K=(XA1(c))2/XA12(c) approaches a limiting value as βε becomes large. Here, we assume that Δ2 ≫ Δ1 and Δ2Δ12, which is motivated by Eqs. (62) and (63) and supported by our calculations. For the range of compositions we studied, K depended weakly on the mixture composition in our numerical calculations but was approximately 1. Consistent with our calculations and simulations (Fig. 7), smaller colloid number density ρc favors loop formation (larger χ). χ also depends on the number of colloid bonding sites nc relative to the number of potential double bonding pairs νc and, significantly, on the ratio of bond volumes Δ12/Δ2. This suggests a potential strategy for limiting loop formation. The number of double bonding pairs νc depends on the linker length, and Δ2 is also roughly proportional to p(rA1A2), the probability of the linker having an end-to-end vector commensurate with the distance between colloid bonding sites ( Appendix B). Reducing the compatibility between the linker and the typical distance between colloid bonding sites could decrease both νc and p(rA1A2), and as a result χ. This might be achieved by modifying the linker’s length or flexibility and presents an opportunity for engineering colloidal self-assembly.

Despite qualitatively capturing the ε and ηc dependences of χ, our TPT calculations consistently overestimated χ compared to the simulations. We note that the simulations at smaller colloid volume fractions (ηc ≲ 0.06) phase separated,16 but our TPT calculations assume that the mixture is homogeneous. Some differences might then be expected between TPT and simulations in this regime; however, similar differences were also obtained at larger volume fractions where the simulated structures were homogeneous. Hence, we suspected there might be additional bonding motifs impacting the thermodynamics that we neglected in our approximation of the irreducible graph sum Δc. In particular, our simulations suggested that “rings,” where two linkers redundantly bridged two colloids [Fig. 1(b)], were also prevalent, but only double-bond graphs [Fig. 8(a)] representing loops were included in Eq. (24).

FIG. 8.

(a) Double-bond (“loop”) graph between one colloid (larger red hyperpoint) and one linker (smaller blue hyperpoint). (b) Ring graph with two colloids and two linkers. Both graphs require simultaneously bonding at two sites on the colloids and linkers. The bonding sites are labeled with general indices, and for clarity, only the fAB-bonds (dashed lines) are drawn.

FIG. 8.

(a) Double-bond (“loop”) graph between one colloid (larger red hyperpoint) and one linker (smaller blue hyperpoint). (b) Ring graph with two colloids and two linkers. Both graphs require simultaneously bonding at two sites on the colloids and linkers. The bonding sites are labeled with general indices, and for clarity, only the fAB-bonds (dashed lines) are drawn.

Close modal

Rings of two colloids linked by two linkers can be incorporated into our TPT using the approximation of Sear and Jackson,22 which adds irreducible graphs [Fig. 8(b)] to Δc,

14ABΓc|AB|=2CDΓc|CD|=2d1d2d3d4σΓc\AB(c)(1)   ×σΓc\CD(c)(2)σ(l)(3)σ(l)(4)gR(4)(1,2,3,4)   ×16fAE(cl)(1,3)fCF(cl)(2,3)fDG(cl)(2,4)fBH(cl)(1,4),
(67)

where gR(4) is the four-body correlation function in the reference mixture approximated by the superposition,

gR(4)(1,2,3,4)gR(cl)(1,3)gR(cl)(2,3)gR(cl)(2,4)×gR(cl)(1,4)eR(cc)(1,2).
(68)

The sums run over all pairs of bonding sites on the colloids labeled 1 and 2, {E, F} and {G, H} are the bonding sites on the linkers labeled 3 and 4, respectively, and the factor of 16 accounts for permutations of bonding within these sets of sites. Because this ring graph involves only singlet densities not bonded at two sites, it will contribute only to ΔcA12(c) and ΔcB12(l) in a similar manner to the double-bond graphs. A similar approximation was also made by Haghmoradi et al. for mixtures of ring-forming colloids with two bonding sites in an article that appeared after our work was completed.57 

We carried out the functional differentiation of Δc including these ring graphs, and we integrated the chemical equilibrium equations for a spatially homogeneous fluid. Two-colloid rings are only likely to form at colloid bonding sites that are nearest neighbors because the linkers would need to stretch significantly (Rdc) to bond at sites on opposite hemispheres. Hence, the equivalent σα(c) remains the same as before, and we need only add terms for the ring graph to Eqs. (60) and (61). The result is

XΓc\A12(c)XΓc(c)XΓc\A1(c)XΓc(c)2=2ρlXB12(l)Δ2+16νc2ρcXA12(c)ρlXB12(l)2Δ4,
(69)
1XB12(l)XB1(l)XB12(l)2=2νcρcXA12(c)Δ2+16νc22ρcXA12(c)2ρlXB12(l)Δ4,
(70)

with the ring-bond volume

Δ4=Ωc2dr12dr13dr14dΩ1dΩ2dR3dR4×p(R3)p(R4)gR(4)(r12,r13,r14,R3,R4)×fAE(cl)(r13,Ω1)fCF(cl)(r12,r13,Ω2,R3)×fDG(cl)(r12,r14,Ω2)fBH(cl)(r14,Ω1,R4).
(71)

Details of how we approximated Δ4 are given in  Appendix B.

As in Sec. IV D, we proceeded to solve the new chemical equilibrium equations. The fractions of linkers with zero, one, or two bonded ends predicted by the TPT with both loops and rings (red line in Fig. 5) were very similar to those predicted by the TPT with only loops, and both agreed well with the simulations. We then determined the fractions of linkers in either loops or rings. With rings included in the TPT, χ given by Eq. (65) is now the total fraction of linkers in either motifs. We accordingly separated the loop-graph contributions to ΔcB12(l) and χ,

χloop=2νcρcXA12(c)XB12(l)Δ2,
(72)

and by subtraction, the ring-graph contribution was χring = χχloop.

The TPT-predicted fraction of linkers in rings (blue dashed line in Figs. 6 and 7) is in very good agreement with the simulations (blue squares). The predicted fraction of loops also decreased in the TPT with loops and rings (red line) compared to the TPT with only loops. This improved the agreement between the TPT and simulations, although the fraction of linkers in loops was still slightly overpredicted in the TPT. We believe that this overprediction might be due to the neglect of other competing bonding motifs such as larger rings or of other irreducible graphs that include effects from steric hindrance.52 Nonetheless, the agreement between the simulations and the TPT with loops and two-colloid rings is good, and the TPT with loops and rings is a significant improvement over TPT1, which predicted that χ = 0 for both bonding motifs.

We have extended Wertheim’s thermodynamic perturbation theory (TPT) for fluids with strong directional attractions to include double bonds between an arbitrary number of pairs of molecular bonding sites in a multicomponent mixture. This extension, which relaxes restrictions in prior TPTs on the number of potential double-bond site pairs, was required to model the assembly of colloidal particles (“colloids”) and difunctional flexible chain molecules (“linkers”). We showed that the fraction of linkers that “looped” to make a double bond with both ends attached to the same colloid and/or that formed “rings” of two linkers bridging between the same two colloids could be reliably predicted using TPT across a range of compositions and attraction strengths. A large linker loop fraction may inhibit assembly of percolated, self-supporting networks of colloids because looped linkers do not form new bridges between colloids. Our work suggests that the loop fraction can be reduced by making the end-to-end distance of the linker incompatible with the distance between colloid bonding sites, which might be achieved by modifying the linker’s molecular weight or flexibility.58 It would be interesting to use the developed theory to compute not only loop fractions but also phase boundaries in order to demonstrate how loops modify the phase behavior of the mixture, including conditions amenable to gelation or equilibrium gels.7,59

M.P.H. thanks Ryan Jadrich for introducing him to Wertheim’s elegant theory. This research was primarily supported by the National Science Foundation through the Center for Dynamics and Control of Materials: an NSF MRSEC under Cooperative Agreement No. DMR-1720595, with additional support from an Arnold O. Beckman Postdoctoral Fellowship (Z.M.S.) and the Welch Foundation (Grant Nos. F-1696 and F-1848). We acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources.

The data that support the findings of this study are available from the authors upon reasonable request.

Useful relationships between ρα(i), σα(i), and cα(i) can be derived through a formalism of “site operators” (or dual numbers) that permits for the manipulation of equations such as Eq. (11).3 For completeness, we will explicitly derive some of Wertheim’s key results using this formalism. Let εA(i)(1) be an operator of a site A on a molecule of component i at a root point labeled 1. All site operators commute and are defined to satisfy

εA(i)2=0.
(A1)

For some α ⊆ Γi, define

εα(i)=AαεA(i),
(A2)

and as a result,

εα(i)εγ(i)=εαγ(i),αγ=0,otherwise.
(A3)

Wertheim considered numbers of the form

ẍ(i)=x(i)+αΓiαxα(i)εα(i)
(A4)

that follow the arithmetic and algebra for real numbers with the exception that division by ẍ(i) having x(i)=0 is not permitted. Analytical functions of ẍ(i) can be computed by series expansion.

With this algebra, ρ̈(i) is defined by Eq. (A4),

ρ̈(i)=ρ(i)+αΓiαρα(i)εα(i),
(A5)

and we will now connect σ̈(i) to ρ̈(i). We start with this same definition for σ̈(i) and use Eq. (12) to replace σα(i),

σ̈(i)=ρ(i)+αΓiαγαργ(i)εα(i)=ρ(i)+αΓiαρ(i)+ρα(i)+γαγργ(i)εα(i)=ρ(i)+αΓiαρ(i)+ρα(i)+γΓiγργ(i)εγ(i)εα(i)=ρ(i)+αΓiαρα(i)εα(i)1+αΓiαεα(i)=ρ̈(i)AΓi1+εA(i).
(A6)

We made use of the identity

αΓiαεα(i)γΓiγxγ(i)εγ(i)=αΓiαεα(i)γαγxγ(i)
(A7)

in the third line and

AΓi1±εA(i)=1+αΓiα(±1)|α|εα(i)
(A8)

in the last line.

It is now straightforward to invert this relationship,

ρ̈(i)=σ̈(i)AΓi1εA(i),
(A9)

using the series expansion (1+εA(i))1=1εA(i). To find expressions for ρα(i) in terms of σα(i), we expand σ̈(i),

ρ̈(i)=σ(i)+αΓiασα(i)εα(i)1+αΓiα(1)|α|εα(i),
(A10)

and collect coefficients of εα(i). This can only include terms in subsets γα, giving Eq. (13).

To find expressions for cα(i) in terms of σα(i), we start from the definition of c̈(i) and consider its exponential

expc̈(i)c(i)=expαΓiαcα(i)εα(i)=αΓiα1+cα(i)εα(i),
(A11)

where in the second equality, we have separated the exponential and used the series expansion exp(εα(i))=1+εα(i). The product can be re-expressed as a sum over partitions of α because only products of terms in disjoint sets will survive the expansion,

expc̈(i)c(i)=1+αΓiαγP(α)aγca(i)εα(i)=ρ̈(i)ρ(i)=σ̈(i)σ(i)AΓi1εA(i),
(A12)

using Eq. (11) to arrive at the second equality and Eq. (A9) with the replacement ρ(i)=σ(i) to obtain the last equality.

Using Eqs. (A11) and (A12) and then inverting, we can immediately show

σ̈(i)σ(i)=αΓiα1+cα(i)εα(i)AΓi1+εA(i)=αΓiα1+δ|α|,1+cα(i)εα(i).
(A13)

Expanding this product and collecting terms as in Eq. (A12) gives Eq. (16).

Alternatively, we can invert Eq. (A12) directly,

c̈(i)c(i)=lnσ̈(i)σ(i)+AΓiln1εA(i)=lnσ̈(i)σ(i)AΓiεA(i),
(A14)

where the last equality used the series expansion ln(1εA(i))=εA(i). In order to determine expressions for cα(i), we now series expand the logarithm,

lnσ̈(i)σ(i)=ln1+αΓiασα(i)εα(i)σ(i).
(A15)

Recall that in general,

ln(1+x)=n=1(1)n1xnn,
(A16)

so the expansion can be performed in multiple variables by collecting all terms involving disjoint sets whose union is α, i.e., these sets must be partitions of α,

lnσ̈(i)σ(i)=αΓiαγP(α)(1)|γ|1|γ|!|γ|aγσa(i)σ(i)εα(i).
(A17)

The factor of |γ|! accounts for all permutations of γ that occur in the expansion. This gives Eq. (15) for cα(i).

We will now show one additional relationship that is useful for computing the thermodynamics. We start by explicitly expanding the sum

αΓiσΓi\α(i)cα(i)=σΓi(i)c(i)+αΓiασΓi\α(i)cα(i)+σ(i)cΓi(i).
(A18)

The second term can be re-expressed as a sum over partitions of Γi, excluding the improper partition, by substituting Eq. (15), and regrouping terms,

αΓiασΓi\α(i)cα(i)=AΓiσΓi\A(i)+αΓiασΓi\α(i)×γP(α)(1)|γ|1(|γ|1)!aγσa(i)σ(i)=AΓiσΓi\A(i)+σ(i)×γP(Γi)γ{Γi}(1)|γ|2|γ|(|γ|2)!aγσa(i)σ(i).
(A19)

Note in the second equality that the sum excludes the improper partition of Γi and that there is a factor of |γ| that accounts for all permutations removing one of the terms in the partition. By direct substitution of Eq. (15), the third term can be similarly expressed,

σ(i)cΓi(i)=σΓi(i)+σ(i)γP(Γi)γ{Γi}(1)|γ|1(|γ|1)!aγσa(i)σ(i).
(A20)

Combining the two sums leads to

αΓiσΓi\α(i)cα(i)=σΓi(i)1+c(i)+Q(i),
(A21)

where Q(i) is given by Eq. (19). Note that Q(i) does not depend on σΓ(i).

To evaluate Δ1, we fixed the colloid at the origin and specified its orientation by a director through A1 having three intrinsic Euler angles (θc, ϕc, ψc)—0 ≤ θc ≤ 2π is a right-handed rotation around the z axis, 0 ≤ ϕcπ is a right-handed rotation about the rotated y axis, and 0 ≤ ψc ≤ 2π is a right-handed rotation about the rotated z axis—so Ωc = 8π2 and θc and ϕc have the same meaning as the azimuthal and polar angles in spherical coordinates. We defined the position rB1 of the bonded linker site B1 in spherical coordinates (rB1,θB1,ϕB1) using the colloid as the origin and rA1 as an axis, i.e., in a frame that rotates with the colloid. We will denote the vector between two sites A and B by rAB = rBrA and the distance by rAB = |rAB|. The end-to-end vector R=rB1B2 was also defined in spherical coordinates (rB1B2,θB1B2,ϕB1B2) using B1 as the origin and rB1 as the axis. These coordinates are illustrated in Fig. 9(a).

FIG. 9.

Schematic of linker coordinates (blue) used to evaluate (a) Δ1 and (b) Δ2 in the frame that rotates with the colloid (red). Arrows designate integration variables, solid lines with bars designate dependent variables, and dotted lines designate the axis of the spherical coordinate system. Refer to the text for the meaning of each variable.

FIG. 9.

Schematic of linker coordinates (blue) used to evaluate (a) Δ1 and (b) Δ2 in the frame that rotates with the colloid (red). Arrows designate integration variables, solid lines with bars designate dependent variables, and dotted lines designate the axis of the spherical coordinate system. Refer to the text for the meaning of each variable.

Close modal

The distance between A1 and B1 was

rA1B12=(dcl*)2+rB122dcl*rB1cosϕB1,
(B1)

and with the potential given by Eq. (56), fA1B1 was nonzero rA1B1<2.5, or equivalently

cosϕB1>x(rB1)=rB12+(dcl*)2(2.5)22dcl*rB1.
(B2)

The end-to-end distance was limited by the segment diameter and the contour length, dlrB1B2L. Conformations where the linker wraps around the colloid are either improbable or forbidden by the hard-sphere exclusions, so we limited 0ϕB1B2π/2, effectively approximating the colloid as a flat surface. With these considerations and after integration over all colloid orientations, θB1, and θB1B2, the integral for Δ1 [Eq. (62)] was

Δ1(2π)2dclr*drB1rB12x(rB1)1dcosϕB1dlLdrB1B2rB1B22×01dcosϕB1B2p(rB1B2)gR(cl)(dcl*,rB2)fA1B1(cl)(rA1B1).
(B3)

We used a linear fit [Fig. 4(a)] for the pair correlation function

gR(cl)(dcl*,rB2)=[ghs(cl)(dcl+)](0.137rB20.0960),
(B4)

with rB2 given by Eq. (64). For numerical convenience, we also fit p(R) for |R| ≥ dl at ηc = 0.01 (Fig. 3) with a piecewise function,

p(r)=pmax+a(rrmax)2,rrmaxpmaxexpb(rrmax)c,r>rmax,
(B5)

where rmax = 2.76, pmax = 3.00 × 10−3, a = −5.24 × 10−4, b = 0.276, and c = 2.45 (all units implicit and consistent with r having units of dl).

To evaluate Δ2, we modified the scheme we used for Δ1 so that B2 was now bonded with a site A2 on the colloid that was a nearest neighbor of A1. The position rB2 of B2 was redefined in spherical coordinates (rB2,θB2,ϕB2) with the colloid as the origin and rA2 as the axis [Fig. 9(b)]. The distance between A1 and B1 was unchanged, and the distance between A2 and B2 was given by an analogous formula. The distance between B1 and B2 was

rB1B22=rB12+rB222rB1rB2sinϕB1sinθB1cosϕB2sinϕB2sinθB2cosϕB1+sinϕB1cosθB1sinϕB2cosθB2.
(B6)

It was then not possible to trivially integrate over θB1 and θB2, and we found that this increased dimensionality added a significant computational cost for evaluating Δ2. However, the distance between the colloid sites rA1A2 is much larger than the range of Eq. (56), so to a good approximation rB1B2rA1A2. As a result,

Δ2p(rA1A2)gR(cl)(dcl*,2dcl*,3π/4)vA1B12,
(B7)

with

gR(cl)(dcl*,2dcl*,3π/4)0.0485[ghs(cl)(dcl+)]2
(B8)

from Fig. 4(b) and

vA1B1=2πdclr*drB1rB12x(rB1)1dcosϕB1fA1B1(cl)(rA1B1).
(B9)

To evaluate Δ4, which is an even higher dimensional integral than Δ2, we adopted a similar approximation for determining p and gR(4) using the positions of the colloid bonding sites rather than the linker ends. We replaced each pairwise gR(cl) in the superposition approximation for gR(4) using Eq. (B4). Labeling the additional colloid bonding sites as C12 = {C1, C2} and linker bonding sites as D12 = {D1, D2}, we obtained

Δ44πvA1B14dc3dcl*+Ldr12r122p(rA1C1)p(rA2C2)×gR(cl)(dcl*,rB1)gR(cl)(dcl*,rB2)gR(cl)(dcl*,rD1)×gR(cl)(dcl*,rD2)fA1B1(cl)(rA1B1)fC1B2(cl)(rC1B2)×fC2D1(cl)(rC2D1)fA2D2(cl)(rA2D2)Ω1,Ω2.
(B10)

We evaluated the integral over the range of colloid separations dcr123dcl*+L, for which the integrand is nonzero based on eR(cc) and p(R), using the trapezoidal rule with 21 uniformly spaced points. The angle brackets denote an unweighted average over the orientations of the two colloids, which we evaluated at each separation r12 using Monte Carlo sampling. We placed one colloid with sites A12 at the origin (0, 0, 0) and the other with sites C12 along the z axis at (0, 0, r12). The bonding sites on both colloids were initially oriented along the +z and +y axes, and we generated 4 × 106 configurations by drawing and applying uniformly random axis–angle rotations to each colloid.60 In averaging, we rejected any configurations in which (1) the distance rA1C1 or rA2C2 was greater than the contour length L or less than the segment diameter dl, (2) either of the site z-coordinates zA1 or zA2 was less than zero, or (3) either of the site z-coordinates zC1 or zC2 was greater than r12. The first condition rejects linker conformations excluded based on p(R), while the other two conditions reject ring configurations that are unlikely or forbidden because the linkers would need to wrap around or penetrate the colloids.

The chemical equilibrium equations can be written as five nonlinear algebraic equations in five variables: XΓc\A12(c), XΓc\A1(c), XΓc(c), XB1(l), and XB12(l). In addition to Eqs. (58) and (59) and Eqs. (69) and (70), the equations from Eq. (43) were

1=11XΓc\A1(c)624XΓc\A12(c)XΓc\A1(c)4XΓc(c)+6XΓc\A12(c)2XΓc\A1(c)2XΓc(c)2+8XΓc\A12(c)3XΓc(c)3XΓc(c)5,
(B11)
XA1(c)=3XΓc\A1(c)512XΓc\A12(c)XΓc\A1(c)3XΓc(c)+10XΓc\A12(c)2XΓc\A1(c)XΓc(c)2XΓc(c)4,
(B12)
XA12(c)=2XΓc\A1(c)4+XΓc\A12(c)XΓc\A1(c)2XΓc(c)+2XΓc\A12(c)2XΓc(c)2XΓc(c)3.
(B13)

The expressions for XA1(c) and XA12(c) were directly substituted before solving.

To obtain numerical solutions, we first evaluated the bond volume integrals for 0 ≤ βε ≤ 20 in increments of 0.5 using the multivariable quadrature method in SciPy (version 1.3.1) with numba (version 0.50.1).61–63 We then iteratively solved the equations for a given composition ηc and ρl/ρc as a function of ε, starting from βε = 0 where all Xα(i)=1 and increasing in steps of 0.01. The bond volumes at intermediate values of ε were interpolated; for βε < 5, we used a cubic spline interpolation of the bond volumes, while for larger ε, we used a cubic spline interpolation of the natural logarithms of the bond volumes. The equations were solved using the SciPy implementation of the Levenberg–Marquadt algorithm with an explicitly specified Jacobian matrix, a relative error tolerance of 10−12 in the sum of squares, a relative error tolerance of 10−12 in the solution, and a maximum of 5000 iterations. We used the solution from the previous value of ε as an initial guess, and we checked that 0Xα(i)1 for both components and all α (including dependent values of Xα(c)) to ensure a physically meaningful converged solution.

1.
M. S.
Wertheim
, “
Fluids with highly directional attractive forces. I. Statistical thermodynamics
,”
J. Stat. Phys.
35
,
19
34
(
1984
).
2.
M. S.
Wertheim
, “
Fluids with highly directional attractive forces. II. Thermodynamic perturbation theory and integral equations
,”
J. Stat. Phys.
35
,
35
47
(
1984
).
3.
M. S.
Wertheim
, “
Fluids with highly directional attractive forces. III. Multiple attraction sites
,”
J. Stat. Phys.
42
,
459
476
(
1986
).
4.
M. S.
Wertheim
, “
Fluids with highly directional attractive forces. IV. Equilibrium polymerization
,”
J. Stat. Phys.
42
,
477
492
(
1986
).
5.
L.
Rovigatti
,
F.
Bomboi
, and
F.
Sciortino
, “
Accurate phase diagram of tetravalent DNA nanostars
,”
J. Chem. Phys.
140
,
154903
(
2014
).
6.
E.
Locatelli
,
P. H.
Handle
,
C. N.
Likos
,
F.
Sciortino
, and
L.
Rovigatti
, “
Condensation and demixing in solutions of DNA anostars and their mixtures
,”
ACS Nano
11
,
2094
2102
(
2017
).
7.
E.
Bianchi
,
J.
Largo
,
P.
Tartaglia
,
E.
Zaccarelli
, and
F.
Sciortino
, “
Phase diagram of patchy colloids: Towards empty liquids
,”
Phys. Rev. Lett.
97
,
168301
(
2006
).
8.
E.
Bianchi
,
P.
Tartaglia
,
E.
La Nave
, and
F.
Sciortino
, “
Fully solvable equilibrium self-assembly process: Fine-tuning the clusters size and the connectivity in patchy particle systems
,”
J. Phys. Chem. B
111
,
11765
11769
(
2007
).
9.
J.
Russo
,
P.
Tartaglia
, and
F.
Sciortino
, “
Reversible gels of patchy particles: Role of the valence
,”
J. Chem. Phys.
131
,
014504
(
2009
).
10.
J.
Russo
,
J. M.
Tavares
,
P. I. C.
Teixeira
,
M. M.
Telo da Gama
, and
F.
Sciortino
, “
Reentrant phase diagram of network fluids
,”
Phys. Rev. Lett.
106
,
085703
(
2011
).
11.
W. G.
Chapman
,
K. E.
Gubbins
,
C. G.
Joslin
, and
C. G.
Gray
, “
Theory and simulation of associating liquid mixtures
,”
Fluid Phase Equilib.
29
,
337
346
(
1986
).
12.
C. G.
Joslin
,
C. G.
Gray
,
W. G.
Chapman
, and
K. E.
Gubbins
, “
Theory and simulation of associating liquid mixtures. II
,”
Mol. Phys.
62
,
843
860
(
1987
).
13.
G.
Jackson
,
W. G.
Chapman
, and
K. E.
Gubbins
, “
Phase equilibria of associating fluids: Spherical molecules with multiple bonding sites
,”
Mol. Phys.
65
,
1
31
(
1988
).
14.
W. G.
Chapman
,
G.
Jackson
, and
K. E.
Gubbins
, “
Phase equilibria of associating fluids: Chain molecules with multiple bonding sites
,”
Mol. Phys.
65
,
1057
1079
(
1988
).
15.
E. A.
Müller
and
K. E.
Gubbins
, “
Molecular-based equations of state for associating fluids: A review of SAFT and related approaches
,”
Ind. Eng. Chem. Res.
40
,
2193
2211
(
2001
).
16.
M. P.
Howard
,
R. B.
Jadrich
,
B. A.
Lindquist
,
F.
Khabaz
,
R. T.
Bonnecaze
,
D. J.
Milliron
, and
T. M.
Truskett
, “
Structure and phase behavior of polymer-linked colloidal gels
,”
J. Chem. Phys.
151
,
124901
(
2019
).
17.
B. A.
Lindquist
,
R. B.
Jadrich
,
D. J.
Milliron
, and
T. M.
Truskett
, “
On the formation of equilibrium gels via a macroscopic bond limitation
,”
J. Chem. Phys.
145
,
074906
(
2016
).
18.
C. A.
Saez Cabezas
,
G. K.
Ong
,
R. B.
Jadrich
,
B. A.
Lindquist
,
A.
Agrawal
,
T. M.
Truskett
, and
D. J.
Milliron
, “
Gelation of plasmonic metal oxide nanocrystals by polymer-induced depletion attractions
,”
Proc. Natl. Acad. Sci. U. S. A.
115
,
8925
8930
(
2018
).
19.
M. N.
Dominguez
,
M. P.
Howard
,
J. M.
Maier
,
S. A.
Valenzuela
,
Z. M.
Sherman
,
J. F.
Reuther
,
L. C.
Reimnitz
,
J.
Kang
,
S. H.
Cho
,
S. L.
Gibbs
,
A. K.
Menta
,
D. L.
Zhuang
,
A.
van der Stok
,
S. J.
Kline
,
E. V.
Anslyn
,
T. M.
Truskett
, and
D. J.
Milliron
, “
Assembly of linked nanocrystal colloids by reversible covalent bonds
,”
Chem. Mater.
32
,
10235
10245
(
2020
).
20.
A.
Stukowski
, “
Visualization and analysis of atomistic simulation data with OVITO–the open visualization tool
,”
Modell. Simul. Mater. Sci. Eng.
18
,
015012
(
2010
).
21.
R. P.
Sear
and
G.
Jackson
, “
Thermodynamic perturbation theory for association into douby bonded dimers
,”
Mol. Phys.
82
,
1033
1048
(
1994
).
22.
R. P.
Sear
and
G.
Jackson
, “
Thermodynamic perturbation theory for association into chains and rings
,”
Phys. Rev. E
50
,
386
394
(
1994
).
23.
A.
Galindo
,
S. J.
Burton
,
G.
Jackson
,
D. P.
Visco
, Jr.
, and
D. A.
Kofke
, “
Improved models for the phase behaviour of hydrogen fluoride: Chain and ring aggregates in the SAFT approach and the AEOS model
,”
Mol. Phys.
100
,
2241
2259
(
2002
).
24.
A. S.
Avlund
,
G. M.
Kontogeorgis
, and
W. G.
Chapman
, “
Intramolecular association within the SAFT framework
,”
Mol. Phys.
109
,
1759
1769
(
2011
).
25.
J. M.
Tavares
,
L.
Rovigatti
, and
F.
Sciortino
, “
Quantitative description of the self-assembly of patchy particles into chains and rings
,”
J. Chem. Phys.
137
,
044901
(
2012
).
26.
L.
Rovigatti
,
J. M.
Tavares
, and
F.
Sciortino
, “
Self-assembly in chains, rings, and branches: A single component system with two critical points
,”
Phys. Rev. Lett.
111
,
168302
(
2013
).
27.
B. D.
Marshall
and
W. G.
Chapman
, “
Thermodynamic perturbation theory for associating fluids with small bond angles: Effects of steric hindrance, ring formation, and double bonding
,”
Phys. Rev. E
87
,
052307
(
2013
).
28.
B. D.
Marshall
and
W. G.
Chapman
, “
Molecular theory for the phase equilibria and cluster distribution of associating fluids with small bond angles
,”
J. Chem. Phys.
139
,
054902
(
2013
).
29.
B. D.
Marshall
, “
A general mixture equation of state for double bonding carboxylic acids with ≥2 association sites
,”
J. Chem. Phys.
148
,
174103
(
2018
).
30.
K.
Hoppe
,
E.
Geidel
,
H.
Weller
, and
A.
Eychmüller
, “
Covalently bound CdTe nanocrystals
,”
Phys. Chem. Chem. Phys.
4
,
1704
1706
(
2002
).
31.
W.
Maneeprakorn
,
M. A.
Malik
, and
P.
O’Brien
, “
Developing chemical strategies for the assembly of nanoparticles into mesoscopic objects
,”
J. Am. Chem. Soc.
132
,
1780
1781
(
2010
).
32.
R. J.
Macfarlane
,
B.
Lee
,
M. R.
Jones
,
N.
Harris
,
G. C.
Schatz
, and
C. A.
Mirkin
, “
Nanoparticle superlattice engineering with DNA
,”
Science
334
,
204
208
(
2011
).
33.
S.
Borsley
and
E. R.
Kay
, “
Dynamic covalent assembly and disassembly of nanoparticle aggregates
,”
Chem. Commun.
52
,
9117
9120
(
2016
).
34.
Y.
Wang
,
P. J.
Santos
,
J. M.
Kubiak
,
X.
Guo
,
M. S.
Lee
, and
R. J.
Macfarlane
, “
Multistimuli responsive nanocomposite tectons for pathway dependent self-assembly and acceleration of covalent bond formation
,”
J. Am. Chem. Soc.
141
,
13234
13243
(
2019
).
35.
N.
Marro
,
F.
della Sala
, and
E. R.
Kay
, “
Programmable dynamic covalent nanoparticle building blocks with complementary reactivity
,”
Chem. Sci.
11
,
372
383
(
2020
).
36.
C. A.
Mirkin
,
R. L.
Letsinger
,
R. C.
Mucic
, and
J. J.
Storhoff
, “
A DNA-based method for rationally assembling nanoparticles into macroscopic materials
,”
Nature
382
,
607
609
(
1996
).
37.
A. P.
Alivisatos
,
K. P.
Johnsson
,
X.
Peng
,
T. E.
Wilson
,
C. J.
Loweth
,
M. P.
Bruchez
, and
P. G.
Schultz
, “
“Organization of ’nanocrystal molecules’ using DNA
,”
Nature
382
,
609
611
(
1996
).
38.
H.
Xiong
,
D.
van der Lelie
, and
O.
Gang
, “
Phase behavior of nanoparticles assembled by DNA linkers,
Phys. Rev. Lett.
102
,
015504
(
2009
).
39.
D.
Zanchet
,
C. M.
Micheel
,
W. J.
Parak
,
D.
Gerion
, and
A. P.
Alivisatos
, “
Electrophoretic isolation of discrete Au nanocrystal/DNA conjugates
,”
Nano Lett.
1
,
32
35
(
2001
).
40.
B. D.
Marshall
and
W. G.
Chapman
, “
Thermodynamic perturbation theory for associating molecules
,” in
Advances in Chemical Physics
, edited by
S. A.
Rice
and
A. R.
Dinner
(
John Wiley & Sons, Inc.
,
2016
), Vol. 160, pp.
1
47
.
41.
W.
Zmpitas
and
J.
Gross
, “
Detailed pedagogical review and analysis of Wertheim’s thermodynamic perturbation theory
,”
Fluid Phase Equilib.
428
,
121
152
(
2016
).
42.
H. C.
Andersen
, “
Cluster methods in equilibrium statistical mechanics
,” in
Statistical Mechanics
, Modern Theoretical Chemistry Vol. 5, edited by
B. J.
Berne
(
Springer
,
1977
), pp.
1
45
.
43.

An articulation point of a connected graph is a point that, if removed, will disconnect the graph into two or more graphs.42 An irreducible graph is free of articulation points. For example, a pair of points or any closed cycle of points is irreducible, but three points bonded collinearly are not irreducible because the middle point is an articulation point.

44.

AB denotes that A is a subset of B, including the improper subset A = B.

45.

A partition of set A is a grouping of the elements of A into one or more non-empty sets using every element exactly once. For example, if A = {a, b, c}, then {{a}, {b}, {c}}, {{a}, {b, c}}, {{a, b}, {c}}, {{a, c}, {b}}, and {{a, b, c}} are all partitions of A. P(A) denotes the set of all possible partitions of A. The last partition, into only a single subset {A}, is called an improper partition.

46.

AB = {aA|aB} denotes the set difference, i.e., all elements that are in A but not in B.

47.

Forgiving some abuse of notation, the label of a single site A should be replaced by a set {A} when it represents a set of bonded sites.

48.

|A| denotes the number of elements in set A.

49.

AB denotes that A is a proper subset of B, i.e., there is at least one element of B that is not in A so AB.

50.
T.
Boublík
, “
Hard-sphere equation of state
,”
J. Chem. Phys.
53
,
471
472
(
1970
).
51.
R. P.
Sear
and
G.
Jackson
, “
The ring integral in a thermodynamic perturbation theory for association
,”
Mol. Phys.
87
,
517
521
(
1996
).
52.
M. S.
Wertheim
, “
Thermodynamic perturbation theory of polymerization
,”
J. Chem. Phys.
87
,
7323
7331
(
1987
).
53.
J. D.
Weeks
,
D.
Chandler
, and
H. C.
Andersen
, “
Role of repulsive forces in determining equilibrium structure of simple liquids
,”
J. Chem. Phys.
54
,
5237
5247
(
1971
).
54.
M.
Fuchs
and
K. S.
Schweizer
, “
Structure of colloid–polymer suspensions
,”
J. Phys.: Condens. Matter
14
,
R239
R269
(
2002
).
55.
S.
Plimpton
, “
Fast parallel algorithms for short-range molecular dynamics
,”
J. Comput. Phys.
117
,
1
19
(
1995
).
56.
G. S.
Grest
and
K.
Kremer
, “
Molecular dynamics simulation for polymers in the presence of a heat bath
,”
Phys. Rev. A
33
,
3628
3631
(
1986
).
57.
A.
Haghmoradi
,
B. D.
Marshall
, and
W. G.
Chapman
, “
Beyond Wertheim’s multi-density theory: Steric hindrance and associated rings in a two-density formalism for binary mixtures of molecules with two associating sites
,”
J. Chem. Eng. Data
,
65
,
5743
5752
(
2020
).
58.
M. P.
Howard
,
Z. M.
Sherman
,
A. N.
Sreenivasan
,
S. A.
Valenzuela
,
E. V.
Anslyn
,
D. J.
Milliron
, and
T. M.
Truskett
, “
Effects of linker flexibility on phase behavior and structure of linked colloidal gels
,” arXiv:2011.12512 (
2020
).
59.
E.
Zaccarelli
, “
Colloidal gels: Equilibrium and nonequilibrium routes
,”
J. Phys.: Condens. Matter
19
,
323101
(
2007
).
60.
R. E.
Miles
, “
On random rotations in R3
,”
Biometrika
52
,
636
639
(
1965
).
61.
C. R.
Harris
,
K. J.
Millman
,
S. J.
van der Walt
,
R.
Gommers
,
P.
Virtanen
,
D.
Cournapeau
,
E.
Wieser
,
J.
Taylor
,
S.
Berg
,
N. J.
Smith
,
R.
Kern
,
M.
Picus
,
S.
Hoyer
,
M. H.
van Kerkwijk
,
M.
Brett
,
A.
Haldane
,
J. F.
del Río
,
M.
Wiebe
,
P.
Peterson
,
P.
Gérard-Marchant
,
K.
Sheppard
,
T.
Reddy
,
W.
Weckesser
,
H.
Abbasi
,
C.
Gohlke
, and
T. E.
Oliphant
, “
Array programming with NumPy
,”
Nature
585
,
357
362
(
2020
).
62.
P.
Virtanen
,
R.
Gommers
,
T. E.
Oliphant
,
M.
Haberland
,
T.
Reddy
,
D.
Cournapeau
,
E.
Burovski
,
P.
Peterson
,
W.
Weckesser
,
J.
Bright
,
S. J.
van der Walt
,
M.
Brett
,
J.
Wilson
,
K. J.
Millman
,
N.
Mayorov
,
A. R. J.
Nelson
,
E.
Jones
,
R.
Kern
,
E.
Larson
,
C. J.
Carey
,
İ.
Polat
,
Y.
Feng
,
E. W.
Moore
,
J.
VanderPlas
,
D.
Laxalde
,
J.
Perktold
,
R.
Cimrman
,
I.
Henriksen
,
E. A.
Quintero
,
C. R.
Harris
,
A. M.
Archibald
,
A. H.
Ribeiro
,
F.
Pedregosa
,
P.
van Mulbregt
and
SciPy 1.0 Contributors
, “
SciPy 1.0: Fundamental algorithms for scientific computing in Python
,”
Nat. Methods
17
,
261
272
(
2020
).
63.
S. K.
Lam
,
A.
Pitrou
, and
S.
Seibert
, “
Numba: A LLVM-based Python JIT compiler
,” in
Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, LLVM ’15
,
2015
.