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.