Self-Assembly of Patchy Particles into Polymer Chains: A Parameter-Free Comparison between Wertheim Theory and Monte Carlo Simulation

0.58 MB·English

by  F. Sciortino

Self-Assembly of Patchy Particles into Polymer Chains: A Parameter-Free
Comparison between Wertheim Theory and Monte Carlo Simulation
Francesco Sciortino,1 Emanuela Bianchi,2 Jack F. Douglas,3 and Piero Tartaglia2
1 Dipartimento di Fisica and INFM-CNR-SOFT,
Universit`a di Roma La Sapienza, Piazzale A. Moro 2, 00185 Roma, Italy
2 Dipartimento di Fisica and INFM-CNR-SMC, Universit`a di Roma La Sapienza, Piazzale A. Moro 2, 00185 Roma, Italy
7 3 Polymers Division, National Institutes of Standards and Technology, Gaithersburg, Maryland 20899, USA
0
Wenumerically studya simple fluid composed of particles havinga hard-corerepulsion, comple-
0
mented bytwo short-ranged attractive (sticky) spots at theparticle poles, which provides a simple
2
modelforequilibriumpolymerizationoflinearchains. Thesimplicityofthemodelallowsforaclose

n comparison, with no fitting parameters, between simulations and theoretical predictions based on
a the Wertheim perturbation theory, a unique framework for the analytic prediction of the proper-
J
tiesofself-assembling particle systemsin termsofmolecular parameter andliquid statecorrelation

2 functions. Thistheoryhasnotbeensubjectedtostringenttestsagainstsimulationdataforordering
2 across thepolymerization transition. Wenumerically determinemanyof thethermodynamicprop-
erties governing this basic form of self-assembly (energy per particle, order parameter or average
] fraction ofparticles intheassociated state, averagechainlength,chainlength distribution,average
t
f end-to-end distance of the chains, and the static structure factor) and find that predictions of the
o
Wertheim theory accord remarkably well with thesimulation results.
s
.
at PACSnumbers: 81.16.Dn,61.20.Ja,61.20.Qg,82.70.Gg: Version: Jan19
m
- I. INTRODUCTION tion.
d
n As a starting point for this type of fundamental inves-
o Recently, there has been great interest in exploiting tigation of self-assembly, we consider the problem of the
c self-assemblytocreatefunctionalnanostructuresinman- equilibrium polymerization of linear polymer chains [17,
[
ufacturing, and this challenge has stimulated a great 18, 19], which is arguably the simplist variety of ther-

1 deal of experimental and theoretical activity [1, 2, 3, mally reversible particle assembly into extended objects
v 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. Self-assembly has (polymersin the caseofourmolecularmodel). Toinves-
1 been considered for over 50 years to be central to un- tigatethisproblemanalytically,weexploittheWertheim
3
derstanding structure formation in living systems and thermodynamicperturbationtheory(W-TPT),whichof-
5
1 modeling and measurements of naturally occurring self- fersa systematic molecularly-basedframeworkfor calcu-
0 assembling systems has long been pursued in the biolog- lating the thermodynamic properties of self-assembling
7 ical sciences [13]. Even the term self-assembly derives systems, although this theory has rarely before been ap-
0 from an appreciation of the capacity of viruses to spon- pliedtothispurpose[20]. TheWertheimtheoryhasbeen
t/ taneously reconstitute themselves from their molecular previously considered to better understand properties of
a components [15, 16], much as in the familiar example of associating fluids and a similar sticky sphere model to
m
micelle formation by block copolymers, lipids, and other the one we consider below has been considered to deter-
- surfactantmoleculesexhibiting amphiphilic interactions. mine how particle association affects critical properties
d
n Thediversityandmorphologicalandfunctionalcomplex- associated with fluid phase separation (critical temper-
o ityofviruses,andthevastnumberofbiologicalprocesses ature and composition, binodals, critical compressibility
c thatformbyasimilarprocessinliving systems,pointto factor, etc. [21, 22, 23], to determine the effect of asso-
v: thepotentialofthistypeorganizationalprocessforman- ciation on nucleation [24] and model antigen-antibody
i ufacturing new materials. While the potential of self- bonding [25]. In contrast, we are concerned here with
X assembly as a manufacturing strategy is clear, our un- thethermodynamictransitionthataccompaniestheself-
r derstandingofhowthis processactuallyworksis stillin- assemblyoftheparticlesintoorganizedstructuresdueto
a
completeandmanyofthebasicprinciplesgoverningthis their anisotropic interactions.
type of organization are unclear. An evolutionary (trial The Wertheim theory is certainly not a unique theory
and error) approach to this problem is not very efficient of the thermodynamics of self-assembling particle sys-
formanufacturing. Thereisevidentlyaneedfordevelop- tems. Models of equilibrium polymerization of linear,
ing a first principles understanding of this phenomenon branched and compact structures have all been intro-
where no free parameters are involved in the theoreti- ducedbasedontheconceptofanassociationequilibrium
cal description to elucidate the fundamental principles being established between the assembling particles and
governing self-assembly and the observables required to comparison of this type of theory to simulation has led
characterize the interactions governing thermodynamic to remarkably good agreement [13, 26, 27, 28, 29, 30,
self-assembly transitions, at least for simple model sys- 31, 32, 33]. Up to the presenttime, however,it has been
tems that can be subjected to high resolution investiga- necessaryto adjust the entropyof associationparameter
2
units of the potential depth (i.e. Boltzmann constant
k = 1). Geometric considerations for a three touching
B
spheres configuration show that the choice of well-width
δ =0.5( 5 2√3 1) 0.119guaranteesthateachsite
− − ≈
isengagedatmostinonebond. Hence,eachparticlecan
p
be form only up to two bonds and, correspondingly, the
lowest energy per particle is u .
0
−
The choice of a simple square-well interaction model
FIG. 1: Pictorial representation of the model studied. Par- todescribetheassociationprocessbetweendifferentpar-
ticles are modeled as hard-core spheres (grey large sphere of ticles is particularly convenient from a theoretical point
diameter σ), decorated by two sites located on the surface of view. It allows for a clear definition of bonding and a
along a diameter(centersof thesmall gold spheresof diame- clear separation of the bond free energy in an energetic
ter2δ). Thevolumeofthegold sphereoutsidethegreylarge andentropiccontributions,beingunambiguouslyrelated
sphereisthebondingvolume. Whenthesiteofanotherparti-
to the depth of the well and to the bonding volume, re-
clesislocatedwithinthebondingvolume,thepairinteraction
spectively.
energy is taken equal to−u0.
III. WERTHEIM THEORY
in this class of theories, so that the modeling is not re-
allyfullypredictive(seediscussioninSect. IIIwherethis
quantityisexplicitlydeterminedfromWertheimtheory). The first-order Wertheim thermodynamic perturba-
The unique aspect of W-TPT is that all the interaction tion theory [34, 35, 37] provides an expression for the
parametersofthistheorycanbedirectlycalculatedfrom free energy of associating liquids. The Helmholtz free
knowledge of the intermolecular potential and standard energy is written as a sum of the HS reference free en-
liquid state correlation functions, so that the theory is ergy AHS plus a bond contribution Abond, which derives
fully predictive. It has also been shown that W-TPT is by a summation over certain classes of relevant graphs
formally equivalent to association-equilibrium models of in the Mayer expansion [37]. In the sum, closed loops
self-assembly [20], so that the Wertheim theory also of- graphs are neglected. The fundamental assumption of
fers the prospect of being able to improve the predictive W-TPT is that the conditions of steric incompatibilities
characteroftheseothertheoriesifthetheoryitselfcanbe are satisfied: (i) no site can be engaged in more than
validatedasbeingreliable. TheWertheimtheoryitselfis one bond; (ii) no pair of particles can be double bonded.
based on a formal perturbation theory [34, 35, 36] and These steric incompatibilities are satisfied in the present
therearenaturallyquestionsabouttheaccuracythatcan modelthankstothelocationofthetwositesandthecho-
be expected from this theory. The present paper consid- sen δ value. In the formulation of Ref. [38], for particles
ersastringenttestofWertheimtheoryasamodelofthe with two identical bonding sites,
thermodynamics of self-assembly by comparing precise
βA
numericalMC data for the thermodynamic propertiesof bond =2lnX X +1 (2)
N −
our model associating fluid to the analytic predictions
of the Wertheim theory where there are no free parame- Here β = 1/k T and X is the fraction of sites that are
B
ters in the comparison. Notably, many of the properties not bonded. X is calculated from the mass-action equa-
that we consider have never been considered before in tion
Wertheim theory.
1
X = (3)
1+2ρX∆
II. TWO PATCHY SITES PARTICLE MODEL
where ρ = N/V is the particle number density and ∆ is
defined by
Wefocusonasystemofhard-sphere(HS)particles(of
diameterσ,theunitoflength)whosesurfaceisdecorated
∆=4π g (r ) f(12) r2 dr (4)
by M =2 identical sites oppositely located (see Fig. 1). HS 12 h iω1,ω2 12 12
The interaction V(1,2) between particles 1 and 2 is Z
Here g (12) is the reference HS fluid pair corre-
HS
V(1,2)=VHS(r12)+ VW(ri1j2) (1) leaxtpio(nVfun(rc1t1io)/nk, Tth)e 1M,aanyedr ff(-1f2u)nction[3i9s]rfe(p1r2e)sen=ts
iX=1,2jX=1,2 − W 12 B − h iω1,ω2
anangle-averageoverallorientationsofparticles1and2
where V is the hard-sphere potential, V (x) is a at fixed relative distance r . Since all bonding sites are
HS W 12
square-wellinteraction (of depth u for x δ, 0 other- identical (same depth and width of the square-well in-
0
wise) and r12 and rij are respec−tively the≤vectors join- teraction), ∆ refers to a single site-site interaction. The
12
ing the particle-particle centers and the site-site (on dif- number of attractive sites on each particle is encoded in
ferent particles) locations. Temperature is measured in the factor 2 in front of ∆ in Eq. 3. In the W-TPT, the
3
resulting free energy is insensitive to the location of the The last expression shows that, within Wertheim the-
attractive sites, i.e. to the bonding geometry of the par- ory, bonding can be seen as a chemical reactionbetween
ticle. Note thatthe angle averagedMayerf-functionco- two unreacted sites forming a bonded pair. In this lan-
incides with the bonding interaction contribution to the guage the quantity K 2∆ plays the role of equilib-
b
≡
virial coefficient. At low T (i.e. βu>>1) the hard-core riumconstant(inunitofinverseconcentration). Writing
contributiontothevirialbecomesnegligibleascompared ρK = exp β(∆U T∆S ) (introducing the energy
b b b
{− − }
tothebondingcomponent. Inthislimit,itisalsopossible andentropychangeinthebondprocess),itispossibleto
toassumeexp(βu ) 1 exp(βu ),sothattheaveraged provide precise expressions for ∆U and ∆S within the
0 0 b b
− ≈
Mayerf-functioncanbeapproximatedwiththe virialas Wertheim theory. Specifically, when eβu0 >> 1 (a very
wellaswiththeintegraloftheBoltzmannfactoroverthe minorapproximationsinceaggregationrequiresT <<u
0
bond volume[40]. to be effective) it is possible to identify
For a site-site square-wellinteraction,the Mayerfunc-
tion can be calculated as [39] ∆Ub = u0 (12)
−
σ+δ
hf(12)iω1,ω2 =[exp(βu0)−1]S(r) (5) ∆Sb =ln 8πρ gHS(r)S(r)r2dr (13)
" Zσ #
where
(note that ∆S is in dimensionless entropy units, since
b
(δ+σ r)2(2δ σ+r) Boltzmann’s constant is taken to equal 1.). If g
S(r)= − − (6) HS ≈
6σ2r 1, ∆S = ln(2NV /V). Hence the change in energy is
b b
givenby the bond energy, while the change in entropy is
is the fraction of solid angle available to bonding when
essentiallyprovidedbythelogarithmoftheratiobetween
two particles are located at relative center-to-center dis-
the bonding volume and the volume per site (V/2N).
tance r. Thus the evaluation of ∆ requires only an ex-
pression for g (r ) in the range where bonding occurs
HS 12
(σ < r < σ +δ). We have used the linear approxima- IV. CLUSTER SIZE DISTRIBUTIONS AND
tion [41] ASSOCIATION PROPERTIES
1 0.5φ 9φ(1+φ) r σ
g (r)= − − (7) It is interesting to discuss the prediction of the
HS (1 φ)3 − 2 (1 φ)3 σ
− − (cid:20) (cid:21) Wertheimtheoryintermofclustersofphysicallybonded
particles [43, 44]. In the case of square-well interactions
(where φ= πσ3ρ) which provides the correctCarnahan-
6 (differently from a continuous potential) we can define a
Starling [42] value at contact. This gives
bond between two particles unambiguously. Evidently,
when there is a bond the interaction energy equals - u .
V (eβu0 1) 0
b
∆= − (8) To make the discussion more transparent, we can de-
(1 φ)3 ×
− fine pb 1 X as the probability that an arbitrary site
≡ −
5 3+8δ+3δ2 3 12δ+5δ2 is bonded. It is thus easy to convince oneself that the
1 φ φ2 number density of monomers is ρ = ρ(1 p )2 = ρX2,
" − 2(cid:0) (15+4δ) (cid:1) − 2(cid:0)(15+4δ)(cid:1) # 1 − b
sincebothsitesmustbeunbonded[20]. Similarlyachain
of l particles has a number density ρ equal to
where we have defined the spherically averaged bonding l
vthoelusmpeecVifibc≡va4luπeoσσf+δδsStu(rd)ire2ddhre=re,πVδ4=(150.+0040δ3)3/23208.5Fσo3r. ρl =ρ(1−pb)2pb(l−1) =ρX2(1−X)(l−1) (14)
R b
At low φ, g (r) tends to the ideal gas limit value
HS since one site of the first and one of the last particle
g (r) 1. In this limit
HS in the chain must be unbonded and l 1 bonds link the
≈
−
l particles.
∆=V (eβu0 1). (9)
b
−
Once the cluster size distribution of chains is known,
Eq. (3) can be easily solved, providing the T and ρ
it is possible to calculate the average chain length L as
dependence of X as
theratiobetweenthetotalnumberdensityandthenum-
2 ber density of chains in the system, (i.e., as the ratio
X = , (10) between the first l1 and the zero l0 moments of the
1+√1+8ρ∆ h i h i
ρ distribution)
l
which has a low T limit X √2ρ∆. ∞
≈ lρ 1
In the more transparent chemical equilibrium form, L l=1 l = (15)
Eq. (10) can be written as ≡ P lρl X
1 X where we have substitutPed, by summing the geometric
X−2 =2ρ∆=ρKb. (11) seriesoverallchainlengths, ∞l=1lρl =ρ and ∞l=1ρl =
P P
4
ρX. Thus, using Eq. (10), extrema in the T ρ plane. The location of the max-
−
imum in the specific heat has been utilized to estimate
1+√1+8ρ∆ the polymerization temperature [17, 18, 19, 27, 49, 50]
L= (16)
2 Within the theory, it is also possible to evaluate the
extent of polymerization Φ, defined as the fraction of
At low T, L √2ρ∆ and hence L grows in density as particles connected in chains (chain length larger than
≈
√ρ(ifthedensitydependenceof∆canbeneglected[45]) one), i.e.
and in T as L exp(βu /2).
0
∼ ∞
The potential energy of the system coincides with the lρ ρ
npuarmtibceler Eof/Nboinsds (times −u0). Hence, the energy per Φ= Pl∞l==21lρll =1− ρ1 =1−X2, (20)
where we have uPsed Eq. 14. Φ plays the role of order
∞
E/N = u0 l=1∞(l−1)ρl = u0(1 X)= u0pb (17) parameterforthepolymerizationtransition. Thedensity
− lρ − − − and temperature dependence of Φ are shown in Fig. 3.
P l=1 l
Thecross-overfromthe monomericstateathighT (Φ
ThesameresultPisofcourseobtainedcalculatingE/N as 0) to the polymeric thermodynamic state at low T (Φ≈
≈
∂(βAbond/N) where (βA /N) is given by Eq. (2). The 1) takes place in a progressively smaller T-window on
∂β bond
decreasing ρ
energyapproachesitsgroundstatevalue(E /N = u )
gs 0
−
as
1
E Egs −1 ρ=0.1
−N =u0X ≈u0(2ρ∆) 2 (18) ρ=10−2
0.8
−3
ρ=10
i.e. with an Arrhenius law with activation energy u /2
0
in the low T limit, a signature of independent bonding 0.6 ρ=10−4
sites [46, 47, 48]. −5
Φ ρ=10
From Eq. (17) it is possible to calculate the (constant
0.4 ρ=10−6
volume) specific heat C as
V
−7
ρ=10
X2eβu0ρ 0.2
C =C (19)
V √1+8ρ∆T2 (a)
0
σ+δ 0 0.05 0.1 0.15 0.2
where C = 8πu2 g (r)S(r)r2dr. In the low T
0 HS T
Zσ
and ρ regime, the specific heat becomes CV X/2T2. 1
≈
At each ρ, the specific heat shows a maximum at a fi- T=0.050
T=0.06
0.8 T=0.07
ρ=0.1 T=0.08
80 ρ=10−2 80 0.6 TT==00..0190
ρ=10−3
x 60
ρ=10−4 ma Φ
60 ρ=10−5 Cv40 0.4
B
k ρ=10−6
/ 20
Cv40 ρ=10−7 0.2
0 (b)
-7 -6 -5 -4 -3 -2 -1
10 10 10 10 10 10 10
ρ
20 0
-8 -6 -4 -2
10 10 10 10
ρ
0.05 0.1 0.15
T
FIG. 3: Extent of polymerization Φ vs. temperature T (a)
and density ρ (b). Symbols are GC simulation data.
FIG.2: W-TPTpredictionsfortheT-dependenceofthespe-
cific heat CV for different valuesof densities ρ (Eq. 19). The
inset shows the value of the specific heat at the maximum Indeed, for each ρ, a transition temperature T can
CVmax. be defined as the inflection point of Φ as functiΦon of
T [17, 18, 19, 49, 50]. The locus T (ρ) provides an es-
Φ
nite T (see Fig. 2), which define a lines of specific heat timate of the polymerization line in the phase diagram.
5
The location of the inflection point of the energy E and providea wayto locate the transitiontemperature when
ofΦ aredifferent, sincedE/dT dX/dT (Eq.17), while experimental (or numerical) data are noisy [32].
∼
dΦ/dT XdX/dT (Eq. 20). In other models of equi-
librium∼po−lymerization,incorporatingthermalactivation 3
or chemical initiation TΦ and TCmax coincide [49, 51]. TΦ(ρ)
Anotherestimate ofthe transitVionlineofthis rounded T (ρ)
2.5 Cvmax
thermodynamic transition can be defined as the locus in
the T ρ plane at which Φ = 0.5, i.e. half of the parti-
−
cles are in chain form (the analog of the critical micelle L 2
concentration [52]), i.e. ρ (T ,ρ) = ρ/2. The corre-
1 1/2
sponding temperature T is then given by the solution
1/2
1.5
of the equation X =√0.5, or equivalently
ρ∆=1 √0.5 (21)
− 1 -4 -3 -2 -1
10 10 10 10
In the present model, the Φ = 0.5 locus is a line at con- ρ
stantp ,correspondingtoaconstantvalueoftheproduct
b
ρ [exp(βu0) 1]. FIG.5: Valueoftheaveragechainlengthalongthepolymer-
−
To provide a global view of the polymerization transi- ization linesTΦ(ρ)andTCmax(ρ),accordingtotheWertheim
V
tion, we show in Fig. 4 the Wertheim theory predictions theory.
forthespecificheatmaximumandpolymerizationtransi-
tion lines. Curves becomes progressively more and more A final relevant consideration concerns the expression
similar on cooling. As clearly shown by the simple ex- for the bonding free energy. Under the assumption ofan
pression for T (Eq. 21) the quantity ρ∆ is constant, ideal state of chains, the system free energy A can be
1/2
which implies that lnρ 1/T. This behavior is also ap- written as
∼
proximativelyfound for the loci defined by the inflection ∞
point of Φ and E, as shown in Fig. 4. βA/V = ρ [lnρ 1+(l 1)lnK ] (22)
l l b
− −
l=1
-1 X
10
which accounts for the translational entropy of the clus-
10-2 T max ters as well as for the bond free energy lnKb which
C
10-3 T v is assumed to be linear in the number of bonds. At
1/2 high T only monomers are present, ρl = ρδl1 and hence
10-4 TΦ βA/V =ρ[lnρ 1]. The bonding part of the free energy
0.12 −
βA canthusbeevaluatedasdifferenceofβAandthe
ρ -5 bond
10
0.1 corresponding high T limit, i.e.
-6 T
10 ∞
0.08
-7 βA /V = ρ [lnρ 1+(l 1)lnK ] ρ[lnρ 1]
10 bond l l b
− − − −
-8 0.06 0.10 0.20 Xl=1
10 ρ (23)
InsertingtheWertheimclustersizedistribution(Eq.14),
10 15 20 25
summing over all cluster size and using the definition
1/T
K 2∆ (Eq. 11), one exactly recovers the Wertheim
b
≡
bonding free energy (Eq. 2). The Wertheim theory
FIG. 4: Specific heat maxima and polymerization transition
lines,aspredictedbytheWertheimtheory,intheT−ρplane. can thus be considered as a mean-field theory of chain
Theinsetshowsthelinearscale. Observethesimilarityofthe associations, with the hard-sphere free energy as its
data in Fig. 2-4 to that of Fig. 4 - 6 in Ref. [29] correspond- high-T limit. Differently from other mean-field ap-
ing to equilibrium polymerization in the Stockmayer fluid, a proaches[30, 45, 49, 54], the theory providesalso a well-
Lennard-Jones particle with a superimposed point dipole. defined prescription for calculating the equilibrium con-
stant, which partially account for the structure of the
Itisinterestingtoevaluatethe valueofaveragedegree reference hard-sphere fluid via the gHS(r) contribution
of polymerization L along the polymerization transition in Kb. It is this special feature that makes the theory
line. This information helps estimating the polymeriza- fully predictive.
tion transition itself, but it is also relevant for the re-
cently proposed analogiesbetween polymerizing systems
V. MONTE CARLO SIMULATION
and glass-forming liquids [53]. As shown in Fig. 5 the
transitiontakesplaceforL 2,consistentwithprevious
≈
findings(comparewiththe insetofFig. 7 inRef.[29]for Wehaveperformednumericalsimulationsofthemodel
the Stockmayer fluid). The fact that, at T , L 2 can inthegrand-canonical(GC)ensemble[55]forseveralval-
φ
≈
6
uesofT andofactivitiestoevaluatethestructuralprop-
erties of the system as a function of density andtemper-
ature. We have performed two types of grand canonical
simulations: particle and chain insertion/removal.
Thefirstmethod(particle)isaclassicalGCMCsimu-
lationwheremonomersareindividuallyaddedtoorelim-
inatedfromthe systemwith insertionandremovalprob-
abilities given by
zV
P = min(1, e−β∆E) (24)
insertion
N +1
N eβ∆E
P = min(1, ), (25)
removal
V z
where N is the number of monomers, ∆E is the change
in the system energy upon insertion (or removal) and z
is the chosen activity. We have simulated for about two
million MC steps, where a MC step has been defined as
50000 attempts to move a particle and 100 attempts to
insertordeleteaparticle. Amoveisdefinedasadisplace-
ment in each direction of a random quantity distributed
uniformlybetween 0.05σ andarotationaroundaran-
±
domaxisofrandomangledistributeduniformlybetween
0.1radiant. The box size has been fixedto 50 σ. With
±
this type of simulations we have studied three tempera-
tures(T=0.08,0.09and0.1)andseveraldensitiesrang-
ing from 0.001 up to 0.2 (corresponding to number N of
FIG. 6: Snapshots of a fraction of the simulated system at
particles ranging from N 1000 to N 20000. ρ≈0.0035 and T =0.10 and T =0.055. At this density,the
≈ ≈
The second method is a grand-canonical simulations polymerizationtransitionislocatedatTΦ≈0.08. Thelength
where we had and remove chains of particles (see for ex- of the shown box edge is about 30σ.
ampleRef.[31]). Theinsertionandremovalprobabilities
for a chain of length l are:
study densities ranging from 10−6 up to 0.2 for four low
temperatures(T=0.05,0.055,0.6and0.7)wherechains
P = min(1,zle−(l−1)βfbV e−β∆E) (26) of length up to 400 monomers are observed. We have
insertion
Ni+1 alsostudiedthesamethreetemperatures(T=0.08,0.09
N eβ∆E and 0.1) examined with the particle insertion/removal
l
P = min(1, ) (27)
removal V zle−(l−1)βfb method to compare the two MC approaches.
As a result of the long simulations performed (about
where Nl is the number of chains of size l, ∆E is the threemonthsofCPUtime foreachstatepoint)resulting
change in energy and βfb is the ln of the integral of the average quantities calculated from the MC data (L, Φ,
Boltzmann factor over the bond volume, i.e., E/N) areaffected by less than3 % relativeerror. Chain
length distributions ρ (whose signal covers up to six or-
l
βf =βu ln[2V ] (28)
b 0 b der of magnitude) are affected by an error proportional
−
to log(ρ )whichprogressivelyincreasesondecreasingρ ,
The activity of a chain of l particles is thus written as | l l
reaching 70% at the smallest reported ρ values.
zl =zle−(l−1)βfb, orequivalently intermofchemicalpo- l
tentialµ =lµ (l 1)f ,whereµisthemonomerchem-
l b
− −
ical potential. Using the chain insertion algorithm, we
VI. SIMULATION RESULTS
have followed the system for about 105 MC steps, where
a MC step has now been defined as 50000 attempts to
move a particle (as in the previous type) and 100 at- First,webeginbyshowinginFig.6representativepar-
tempts to insert or delete a chain of randomly selected ticle configurations both above and below the polymer-
length. The geometry of the chain to be inserted is also ization transition temperature, T . Evidently, the parti-
φ
randomly selected. The box size was varied from 50 to cles are dispersed as a gas of monomers and as a gas of
400 σ, according to density and temperature, to guar- chains above and below this characteristic temperature.
antee that the longest chain in the system was always The low-T configurations is composed by semi-flexible
shorter than the box size. At the lowest T, N 50000. chains, with no rings. Indeed, the short interaction
≈
Using chain moves we have been able to equilibrate and rangeδ introducesa significantstiffness inthe chainand
7
-2
10
4
10
-3 T=0.08
10
-4
10
>
2<Re102 ρl10-5
2 1.1
<R > ~ l 10-6
e
2 2
<R > ~ l -7
e 10
0 -8
10 10
0 1 2 3
10 10 10 10
0 10 20 30 40 50
chain length l
Chain length l
FIG. 7: Mean squared end-to-end distance < Re2 > for iso- FIG. 8: Comparison of the chain length distributions ρl
lated chains of length l built according to the model studied
at T = 0.08 obtained independently with the monomer
in thisarticle.
(full symbols) and chain (lines connecting open symbols)
grand-canonicalsimulationsforseveralvaluesofthemonomer
activity z. From left to right, the activity values are:
a persistence length extending over several monomers. 10−3,1.5 × 10−3,2.4 × 10−3,3 × 10−3,5 × 10−3,6 × 10−3.
To quantify the linearity of the chains for the present Similar agreement is also found at the othertemperatures.
model we show in Fig. 7 the chain end-to-end squared
distance<R2 >forisolatedchainsofchainlengthupto
e
l O(103). Single chains are generated by progressively choose the ideal gas as reference state, i.e. we approxi-
∼
adding monomers to a pre-existing chain in a bonding mate the reference radial distribution function with one.
configurations, after checking the possible overlap with In the more realistic approximation, we use the small r
allpre-existingmonomers. Sincethebondinteractionisa expansionofthe hard-sphereradialdistributionfunction
well,allpointsinthebond-volumehavethesamea-priori (see Eq. 7). Comparison between simulation data and
probability. AsshowninFig.7,theend-to-endchaindis- theoreticalpredictions(notethere areno fitting parame-
tance scalesas a power-law(<R2 > l2ν) both atsmall ters)is reportedin Fig.9 for two different temperatures.
e ∼
andlargelvalues,withacrossingbetweenthetwobehav- At low densities (sampled at low T) the approximation
iorsaroundl O(10). At smalll,the chainis persistent g 1 is already sufficient to properly describe ρ . At
≃ HS ≈ l
informandthusisrod-like. Forlargerl,thebest-fitwith higher densities (sampled at higher T), the full theory is
a power law suggests an apparent exponent 2ν 1.1, requested to satisfactory predict the chain length distri-
≈
that is expected to evolve — for very long chains — to- butions.
wardtheself-avoidingvalue2ν 1.18[56,57]. Werecall Fig. 10(a) compares the Wertheim theory predictions
≈
that2ν =6/5intheFlorymeanfieldprediction[58]and for the averagechainlength with the correspondingsim-
2ν =1 in the simple random walk model. ulationresults. In the entire investigatedρ and T range,
Toprovideevidencethatthe chainGC MCsimulation the Wertheim theory provides an accurate description
provides the correct sampling of the configurations (and of the equilibrium polymerization process. The limiting
hencethattheactivityofthechainoflengthliscorrectly growth law in √ρ is clearly visible at the lowest tem-
assigned)wecompareinFig.8thechainlengthdensities peratures. At the highest temperatures, it is possible to
ρ calculatedwiththetwomethodsatT =0.08. Thedis- access the region of larger densities (ρ 0.1) where the
l
∼
tributions calculated with the two different methods are presence of other chains can not be neglected any longer
identical. Similar agreement is also found at T = 0.09 and ∆ becomes ρ dependent. In this limit, the growth
and T = 0.1. This strengthens the possibility of using low L √ρ is no longer obeyed [45]. As a further check
∼
the chain MC method, which does not significantly suf- on Wertheim theory, we collapse all the L data to the
fer from the slow equilibration process associated to the universal functional form predicted by Wertheim theory
increase of the Boltzmann factor exp(βu ) on cooling. using the scaling variable 2∆ρ (Fig. 10(b))
0
All the following data are based on chain GC-MC simu- As an ulterior confirmation of the predictive capabili-
lations. tiesofWertheimtheory,wereportacomparisonbetween
We nextcomparethe chainlengthdistributions calcu- simulations and theory for the density dependent of the
latedusingthechainMCmethodwiththepredictionsof extent of polymerization Φ (Fig. 3) and for the energy
the Wertheim theory. We compare the simulation data perparticle(Fig.11). Bothfiguresclearlyshowsaexcel-
with two different approximation: in the first one we lent agreement between the simulated and the predicted
8
-4
10
T=0.050
-5 T=0.06 (a)
10 T=0.055
T=0.06
-6 101 T=0.07
10
T=0.08
T=0.09
ρl10-7 L T=0.10
-8
10
-9
10
100 (a)
0 50 100 150 -7 -6 -5 -4 -3 -2 -1
10 10 10 10 10 10 10
Chain length l ρ
-2
10 T=0.09 (b)
T=0.050
-3
10 T=0.055
T=0.06
10-4 101 T=0.07
T=0.08
ρl -5 T=0.09
10 L T=0.10
Wertheim
-6
10
-7
10
(b)
-8
10
0 5 10 15 20 25 30 0
10
Chain length l 10-2 10-1 100 101 102 103
2 ρ ∆
FIG.9: Chainlengthdensitiesρlforseveralactivityvaluesat
two different temperatures, T =0.06 and T =0.09. Symbols FIG.10: Averagechainlengthasafunctionofthedensityfor
are simulation data. Lines are Wertheim theory predictions: allstudiedtemperatures. (a): LinesaretheWertheimtheory
dashed lines assume gHS = 1 (Eq. 9 for ∆), while full lines predictions: dashedlinesassumegHS =1(Eq.9for∆),while
are based on the full radial dependence of gHS (Eq. 8 for fulllinesarebasedonthefullradialdependenceofgHS(Eq.8
∆). At the lowest T (and average densities), the ideal gas for ∆). (b): Scaled representation of L vs. 2ρ∆. Symbols
approximation is sufficient to model theMonte Carlo data. aresimulationdata. ThelineisthefunctionL(x)= 1+√1+2x
2
(SeeEq. 16)
ρ dependence of Φ and E at all T investigated.
As a test of the validity of the approach of the poly-
merizingsystemasanidealgasofequilibriumchains,we
1
compare the monomer activity z and the monomer den- S(q) e−iq~·(~ri−~rj) , (29)
≡hN i
sity ρ inFig.12. Indeedthe activityofa clusterofsizel
1 i,j
X
coincides with ρ , which is consistentwith ideal gasscal-
l
ing. In particular, the activity of the single particle (the whereri isthecoordinateofparticlei,thesumrunsover
inputintheMCgrand-canonicalsimulation)canbecom- allN particlesinthesystem,andtheaverage ... isover
h i
paredwiththeresultingdensityofmonomersρ . Datain equilibrium configurations. In the ideal gas limit, cor-
1
Fig.12showsthattheideal-gaslawiswellobeyedatlow relations between different chains can be neglected and
T and ρ, confirming that, in the investigated range, the S(q) can be formally written as,
systemcanbevisualizedasanidealgasmixtureofchains ∞
ρ lS(q)
of different lengths, distributed according to Eq. 14. S(q)= l=1∞l l (30)
ρ l
The existence of a large T-ρ window where an ideal P l=1 l
mixture of chains provides a satisfactory representation
whereS (q)isthestructurePfactor(formfactor)ofachain
of the system suggests that, in this window, correlations l
of length l:
between different chains can be neglected. In this limit,
the structure of the system should be provided by the l
1
structure of a single chain, weighted by the appropriate S (q)= e−iq~·(~ri−~rj) (31)
l
lh i
chain length distribution. Specifically, we have
i,j=1
X
9
0 where < lm > lmρ denotes the m moment of the
≡ l l
cluster size distribution ρ . Fig. 13 shows a comparison
l
P
-0.2 between the S(q) calculated in the simulation and the
theoreticalS(q)evaluatedaccordingtoEq.30and32ata
lowT,wheretheidealgasapproximationisvalid. Devia-
-0.4
tionsareonlyobservedatthehighestdensity,suggesting
N
E/ that the ideal gas of chains is a good representation of
-0.6 T=0.050 thestructureofthesystem,inagreementwiththeequiv-
T=0.055 alence between activity and monomer density shown in
T=0.06
T=0.07 Fig. 12. This observation is particularly relevant, since
-0.8 T=0.08 it suggests that (in the appropriate T-ρ window) also a
T=0.09
T=0.10 description of the dynamics of the model based on the
-1 assumption of independent chains can be attempted.
10-7 10-6 10-5 10-4 10-3 10-2 10-1 100
ρ
T=0.055
30
FIG.11: Energyperparticleinunitofu0 asafunctionofthe z=0.000010
densityforallstudiedtemperatures. Symbolsaresimulations z=0.000012
datawhilelinesaretheWertheimtheorypredictions: dashed z=0.000014
linesassumegHS =1(Eq.9for∆),whilefulllinesarebased z=0.000016
on thefull radial dependenceof gHS (Eq.8 for ∆) 20 z=0.000018
)
q
(
-1 S
10
10
-2
10
-3
10
T=0.050
z 0
T=0.055 0.01 0.1 1
10-4 T=0.06 qσ
T=0.07
-5 T=0.08
10 T=0.09 FIG. 13: Comparison of the structure factor calculated from
T=0.10 the simulated configurations and from Eq. 30-32 at T=0.055
-6 for different values of theactivity.
10
-6 -5 -4 -3 -2 -1 0
10 10 10 10 10 10 10
ρ
1
VII. CONCLUSIONS
FIG. 12: Relation between the activity z and the monomer
density ρ1 for all investigated temperatures. Line indicates
theideal gas law z=ρ1. In our pursuit of a fully predictive molecularly-based
theory ofself-assembly interms ofmolecular parameters
andliquidstatecorrelationfunctions,wehaveconsidered
Since the persistence length of the chains is 10 20
a direct comparison of a liquid in which fluid particles
≈ −
particles (see Fig. 7), one can assume, as a first approx-
have sticky spots on their polar regions to the predic-
imation, that in the investigated T-ρ region chains are
tions of the Wertheim theory for the relevant properties
linear. When this is the case, averaging over all possible
governing the self-assembly thermodynamics. In the in-
orientation of the chain gives:
vestigated region of temperatures and densities (which
l−1 spans the polymerization transition region), the predic-
1 sin(jqσ)
S (q)=1+ 2(l j) . (32) tions of Wertheim theory describe simulation data re-
l
l − jqσ markably well (without the use of any free parameters
j=1
X
in these comparisons). This success means that we can
The small q expansion of Sl(q) is be confident in pursuing more complicated types of self-
assembly based on the foundation of Wertheim theory.
l(l2 1)
S (q) l − (qσ)2 (33) For example, it is possible to extend the Wertheim the-
l
≈ − 36 oryanalysisin adirect wayto describe the self-assembly
Correspondingly, S(q) behaves at small q as of branched chains having multifunctional rather than
dipolar symmetry interactions as in the present paper
<l2 > q2σ2 <l4 > <l2 > and to compare these results in a parameter free fash-
S(q) 1 − , (34)
≈ <l > − 36 <l2 > ion to corresponding MC simulations for the these mul-
(cid:20) (cid:18) (cid:19)(cid:21)
10
tifunctional interaction particles [59]. Recent work has solutionisadifficultproblemthathasbeenaddressedby
shown that the Wertheim theory describes the critical many authors previously (Ref. [69] and refs. therein). It
properties of these mutifunctional interaction particles would clearly be interesting to extend the present work
rather well [60]. According to this recent study, liquid to determine how well the Wertheim theory could pre-
phases of vanishing density can be generated once small dict entropies of association for self-assembly processes
fraction of polyfunctional particles are added to chain- thatoccurinasolventratherthaninthegasphase. The
forming models like the one studied here. With the new most interesting solvent in this connection, water, is a
generation of non-spherical sticky colloids, it should be particular challenge since water itself can be considered
possibleto realize”emptyliquids” [60]andobserveequi- an associating fluid, so that we are confronted with the
librium gelation [46, 61], i.e. approach dynamical arrest problemofhowthe waterassociationcouples tothe par-
under equilibrium conditions. ticle self-assembly. The problem of understanding the
The Werhtheim theory has also been applied success- common tendency of particle self-assembly in aqueous
fullytodescriptionofmolecularassociatedliquids[62,63, solutions to occur upon heating requires particular in-
64] and to the thermodynamics of hard sphere polymer vestigation. In the future, we look forward to exploring
chains with short range attractive interactions [65, 66]. thesemorecomplex mixturesofassociatingfluids,which
Thus, the theory could be adapted to describing mutu- are so prevalent in real biological systems and in a ma-
allyassociatingpolymersandtheformationofthermally terials processing context.
reversible gels in these fluids upon cooling. As a final comment, we note that numerical work on
In summary, the Wertheim theory provides a promis- this class of simple models (playing with the particle in-
ingframeworkfortreatingthethermodynamicsofawide teractionsymmetries)canhelpunderstandingmorecom-
range of self-assembling systems. The development of plicated ordered structures (as for example sheet-like,
thistheoryanditsvalidationbysimulationandmeasure- nanotube and closed nanoshell structures), as recently
mentshouldprovidevaluabletoolsinthepracticaldevel- found when particles have multipole interaction poten-
opmentofself-assemblyasapracticalmeansofsynthetic tials [8, 9, 11, 13]. We also note that the results dis-
manufacturing. This theory also offers the prospect of cussed here apply to the growing field of functionalized
improving the existing equilibrium association theories colloidalparticles,colloidalparticleswithspecificallyde-
thatarelargelybasedonalatticefluidmodelframework. signed shapes and interaction sites [1, 2, 3, 4, 5, 70].
This couldallowprogresstobe mademorerapidly,since
this type of computation often offers computational ad-
vantagesand because many problems suchas chemically
initiatedchainbranching[67]andthermallyactivatedas-
VIII. ACKNOWLEDGMENTS
semblyprocesseshavealreadybeenconsideredbylattice
approaches [50, 68].
The problem of estimating the entropy of association We acknowledge support from MIUR-Prin and
in real self-assembling molecular and particle systems in MCRTN-CT-2003-504712.
[1] V.N.Manoharan,M.T.Elsesser,andD.J.Pine,Science teins, colloids and patchy models (2007), URL
301, 483 (2003). http://www.citebase.org/abstract?id=oai:arXiv.org:cond-mat/070
[2] C. A. Mirkin and et al., Nature (London) 382, 607 [12] F.W.StarrandF.Sciortino, J.Phys.: Condens.Matter
(1996). 18, L347 (2006), cond-mat/0512260.
[3] Y.-S. Cho, G.-R.Yi, J.-M. Lim, S.-H. Kim, V. N. [13] V.WorkumandJ. F.Douglas, Phys.Rev.E 73, 031502
Manoharan, D. J. Pine, and S.-M. Yang, J. Am. Chem. (2006).
Soc. 127, 15968 (2005). [14] S.I.Stupp,S.Son,H.C.Lin,andL.S.Li,Science259,
[4] G. Yi, V. N. Manoharan, E. Michel, M. T. Elsesser, 59 (1993).
S.Yang, and D. J. Pine, Adv.Mater. 16, 1204 (2004). [15] H. Fraenkel-Conrat and R. C. Williams, Proc. Natl.
[5] Y.-S.Cho and et al., Chem. Mater. 17, 5006 (2005). Acad. Sci. U.S.A 41, 690 (1955).
[6] C. Mirkin, R. Letsinger, R. Mucic, and J. Storhoff., Na- [16] P. Buttler, J. Gen. Virol. 65, 253 (1984).
ture382, 607 (1996). [17] S. C. Greer, J. Phys. Chem. B 102, 5413 (1988).
[7] F. W. Starr, J. F. Douglas, and S. C. Glotzer, J. Chem. [18] S. C. Greer, Adv.Chem. Phys. 94, 261 (1996).
Phys.119, 1777 (2003). [19] S. C. Greer, Ann.Rev.Phys. Chem. 53, 173 (2002).
[8] S.C. Glotzer, Science 306, 419 (2004). [20] I.G.Economou andM.D.Donohue,AIChEJ.37,1875
[9] S. C. Glotzer, M. J. Solomon, and N. A. Kotov, AIChE (1991).
Journal 50, 2978 (2004). [21] G. Jackson, W. G. Chapman, and K. Gubbins, Mol.
[10] Z. Zhangand S.C. Glotzer, Nanoletters 4, 1407 (2004). Phys. 65, 1 (1988).
[11] J. P. K. Doye, A. A. Louis, I.-C. Lin, L. R. Allen, [22] N. A. Busch, M. S. Wertheim, and M. L. Yarmush, J.
E. G. Noya, A. W. Wilber, H. C. Kok, and R. Lyus, Chem. Phys.104, 3962 (1996).
Controlling crystallization and its absence: Pro- [23] E. A.Mu¨ller and K. E. Gubbins,Ind.Engr. Chem. Res.