剛?cè)岣叻肿芋w系自組裝對(duì)稱膠束的自洽場(chǎng)模擬

羅永順，李思銳，趙方鑫，王曉淵

貴州大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，貴州貴陽(yáng) 550025

高分子也稱聚合物(polymer)，是由若干小分子通過(guò)共價(jià)鍵重復(fù)鏈接而成，并具有較大相對(duì)分子質(zhì)量的長(zhǎng)鏈分子，它具有復(fù)雜的空間鏈構(gòu)像分布，這一特殊結(jié)構(gòu)賦予高分子材料質(zhì)輕、透明、電絕緣及耐腐蝕等諸多優(yōu)良性能．嵌段共聚物是由兩種或兩種以上的嵌段(其基元是化學(xué)性質(zhì)相同的單體)通過(guò)化學(xué)鍵聚合而成的一類特殊高分子．兩嵌段共聚物是最簡(jiǎn)單的情形．若高分子長(zhǎng)鏈僅由一種單體聚合而成，則稱為均聚物．共聚物在不同條件下會(huì)自組裝形成豐富的微觀結(jié)構(gòu)，這些納米尺度的微觀結(jié)構(gòu)因具有良好的可調(diào)控性，使其在制備各種納米材料等工業(yè)應(yīng)用中發(fā)揮重要的作用[1]．

近年來(lái)，為更好的控制和應(yīng)用高分子共聚物的誘人性質(zhì)，除實(shí)驗(yàn)手段之外，研究?jī)A向于通過(guò)理論或數(shù)值模擬的方法，分析高分子的自組裝等豐富相行為．目前，絕大多數(shù)研究基于柔性(coil)體系的Gauss鏈模型展開(kāi)，并已取得很大成功．但實(shí)際的高分子常帶有一定的剛性(rod)，如共軛高分子、液晶高分子及生物DNA分子等．此時(shí)，柔性分子鏈模型失效，必須引入帶剛性的分子鏈模型(wormlike鏈模型)．最簡(jiǎn)單的帶有剛性分子的共聚物是剛?cè)?rod-coil)兩嵌段共聚物．由于剛性鏈段的各向異性與取向關(guān)聯(lián)性，rod-coil共聚物展現(xiàn)出更豐富復(fù)雜的液晶相等自組裝結(jié)構(gòu)[2]．

自洽平均場(chǎng)理論(self-consistent field theory，SCFT)是研究嵌段共聚物自組裝行為最成功、最廣泛的理論．其核心思想是對(duì)高分子鏈進(jìn)行粗?；幚恚雎哉鎸?shí)高分子基元的復(fù)雜化學(xué)結(jié)構(gòu)，將一條高分子鏈看成空間中運(yùn)動(dòng)“粒子”走過(guò)的路徑．換言之，SCFT引入一個(gè)外場(chǎng)，將一個(gè)多體作用體系變成等效外場(chǎng)下的單體體系，將基于粒子的體系變成基于場(chǎng)的體系．基于wormlike鏈模型的SCFT，求解高分子體系平衡態(tài)的問(wèn)題就轉(zhuǎn)化成求解相應(yīng)自洽場(chǎng)方程組的數(shù)學(xué)問(wèn)題，而這類數(shù)學(xué)問(wèn)題是一類復(fù)雜的非線性、非局部，且具有多解和多參數(shù)的優(yōu)化問(wèn)題．

Rod-coil兩嵌段聚合物自組裝已取得豐富的研究成果．實(shí)驗(yàn)上已發(fā)現(xiàn)多種有趣的微觀結(jié)構(gòu)，如wavy相、zigzags相、arrowheads相、continuous cubic相、六邊堆積的柱狀相、穿孔層狀相及球狀相等[3-6]．在理論上，SEMENOV等[7-8]首次引入微觀模型研究rod-coil共聚物的向列相、層狀相-A和層狀相-C等一維結(jié)構(gòu)．WILLIAMS等[9]將SEMENOV等的模型發(fā)展到二維hockey-puck-like micellar結(jié)構(gòu)，但這些理論都比較粗略，難以分析復(fù)雜結(jié)構(gòu)．在數(shù)值模擬方面，早期的典型代表是實(shí)空間和復(fù)空間方法，它們各有利弊．后來(lái)的重要改進(jìn)是擬譜方法[10]．LIANG等[11]給出rod-coil體系自洽場(chǎng)方程的多種數(shù)值格式，并建立相應(yīng)的高效算法．JIANG等[12]分析多個(gè)不同化學(xué)成分高分子體系的自洽場(chǎng)模型，給出能量上升和下降的方向以及相應(yīng)的自洽場(chǎng)高效算法．GAO等[13]利用有限體積法，獲得四方堆積的柱狀相．此外，分子模擬、動(dòng)力學(xué)方法和蒙特卡洛模擬等方法也應(yīng)用于rod-coil高分子體系的研究[14-16]．2014年，LI等[17]研究rod-wormlike兩嵌段形成的液晶相結(jié)構(gòu)，理論預(yù)測(cè)了改變?nèi)嵝枣湺蔚膭傂詫?duì)液晶相變的影響．2017年，CAI等[18]數(shù)值研究了rod-coil自組裝雙層膜的豐富液晶相行為．2018年，MAO等[19]利用場(chǎng)論蒙特卡洛(field-theoretic Monte Carlo)方法研究半剛性高分子的密度漲落效應(yīng)與相變．

本研究基于wormlike鏈模型，給出巨正則系綜下rod-coil/coil混合體系的自洽場(chǎng)理論．將相應(yīng)自洽場(chǎng)算法運(yùn)用于rod-coil/coil體系下臨界膠束結(jié)構(gòu)的數(shù)值模擬，得到了與coil-coil體系下臨界膠束相似的結(jié)構(gòu)和性質(zhì). 即會(huì)形成不同的膠束結(jié)構(gòu)，且當(dāng)嵌段A的體積分?jǐn)?shù)小于0.54時(shí)，平面膠束的臨界化學(xué)勢(shì)最低；大于0.54且小于0.67時(shí)，柱狀膠束的最低；大于0.67時(shí)，球狀膠束最低. 隨著嵌段A體積分?jǐn)?shù)的增加，膠束結(jié)構(gòu)的穩(wěn)定形態(tài)將從平面轉(zhuǎn)化為柱狀再到球狀．

1 rod-coil/coil混合體系的統(tǒng)計(jì)力學(xué)

首先考慮高分子鏈條數(shù)不變，體積為V的正則系綜．假定體系含nc條rod-coil(AB)兩嵌段共聚物和nh條coil(A)均聚物，每條rod-coil高分子由NA個(gè)A單體和NB個(gè)B單體連接組成，聚合度為N=NA+NB， 每條coil高分子由Nh個(gè)A單體組成，記每個(gè)單體的數(shù)密度為ρ0=(nhNh+ncN)/V， 即體系中每個(gè)單體占有相同體積1/ρ0． 另一條AB分子鏈的構(gòu)型表示為Rkαc(s)，k∈(A, B)，s∈[0,Nk]，αc=1,2,…,nc． 類似地， A分子鏈的構(gòu)型可以表示為Rαh(s)，s∈[0,Nh]，αh=1,2,…,nh． rod-coil/coil混合體系的簡(jiǎn)單表示如圖1.

圖1 Coil均聚物和rod-coil兩嵌段共聚物的示意圖Fig.1 (Color online) Schematic diagram of coil homopolymer and rod-coil diblock copolymer

(1)

(2)

(3)

(4)

(5)

假設(shè)2種不同粒子之間存在相互排斥的Flory-Huggins作用力，而另外相鄰剛性鏈之間有平行排列的趨勢(shì)，采用Maier-Saupe模型描述此種作用力，則體系具有4部分能量：coil鏈段的拉伸能H0、 rod鏈段的彎曲能H1、 兩類粒子之間的排斥能HF及rod鏈之間平行排列的能量HS， 其相應(yīng)的表達(dá)式分別為

(6)

(7)

(8)

(9)

按照統(tǒng)計(jì)力學(xué)原理，首先必須寫(xiě)出多鏈高分子體系的配分函數(shù)．根據(jù)上述定義的序參量以及引入的4部分能量，對(duì)于不可壓體系，可以寫(xiě)出基于分子構(gòu)型Rkαc和Rαh的正則系綜下的配分函數(shù)，為將其表達(dá)為基于場(chǎng)的形式，則由路徑積分的高斯公式(Hubbard-Stratonovich變換)，將粒子形式的配分函數(shù)轉(zhuǎn)化為配分函數(shù)，即

exp(-H(wk,φk,ξ,M))

(10)

其中，H是有效哈密頓量(也稱自由能)．

在巨正則系綜下，體系的分子數(shù)目是變化的，但分子的化學(xué)勢(shì)不變．記均聚物的化學(xué)勢(shì)為μh， AB兩嵌段共聚物的化學(xué)勢(shì)為μc， 則巨正則系綜下的配分函數(shù)為

(11)

將式(10)代入式(11)可得巨正則系綜下的有效哈密頓量為

H(wk,φk,ξ,M)=

eμcVQc-eμhVQh

(12)

其中，Qc和Qh為rod-coil分子鏈和均聚物分子鏈在場(chǎng)wk、φk、M和ξ下的單鏈配分函數(shù)，單鏈配分函數(shù)可通過(guò)對(duì)高分子鏈傳播子q和相應(yīng)的鏈逆向傳播子q+積分得到．且傳播子的計(jì)算可通過(guò)求解擴(kuò)散方程得到[12]，即

wA(r)qA(r,s),s∈[0,NA]

(13)

(14)

(15)

(16)

至此，對(duì)于任意給定場(chǎng)wk、φk、M和ξ， 均可通過(guò)求解擴(kuò)散方程得到傳播子，進(jìn)一步得到單鏈配分函數(shù)Qc、Qh以及有效哈密頓量．

2 平均場(chǎng)近似與自洽場(chǎng)方程

在自洽場(chǎng)理論中，要精確計(jì)算配分函數(shù)或自由能十分困難．所幸式(11)的配分函數(shù)是關(guān)于指數(shù)型泛函的路徑積分，這樣可以利用指數(shù)型積分的Laplace近似或復(fù)變量情形的鞍點(diǎn)近似來(lái)逼近配分函數(shù)．FREDRICKSON等[20]首次將這種近似思想用于討論兩嵌段共聚物體系的自洽場(chǎng)方程，并建立高效的求解方法．

因此，利用鞍點(diǎn)近似思想，可將尋求體系的平衡態(tài)問(wèn)題轉(zhuǎn)化為求自由能泛函的鞍點(diǎn)問(wèn)題，即在鞍點(diǎn)處自由能泛函關(guān)于各個(gè)場(chǎng)的一階變分為0，

(17)

在rod-coil/coil混合體系的自洽場(chǎng)模型中，參數(shù)較多，必須對(duì)模型進(jìn)行量綱為一化．量綱為一化傳播子滿足的修正擴(kuò)散方程為

(18)

(19)

(20)

相應(yīng)的初值條件為

其中，f為體系中高分子鏈的體積分?jǐn)?shù)，即fA、fB及fh分別表示體系中兩嵌段共聚物中的柔性、剛性高分子鏈的體積分?jǐn)?shù)及體系中均聚物中高分子鏈的體積分?jǐn)?shù). 另外， 選取化學(xué)勢(shì)μ0=(μh-ln(ρ0/N))/fh作為標(biāo)準(zhǔn)，并作變換使得體系的有效哈密頓量保持不變，則變換后的自由能形式為

zclnQc-lnQh

(21)

其中，zc=eμc指化學(xué)勢(shì)的活性度(即化學(xué)勢(shì))，而序參量為

(22)

(23)

(24)

于是，相應(yīng)的自洽場(chǎng)方程組為

φA(r)+φB(r)=1

(25)

(26)

(27)

(28)

M(r)=ηNS(r)

(29)

綜上，基于wormlike鏈模型已給出巨正則系綜下rod-coil/coil混合體系的自洽場(chǎng)模型．

3 算法描述與數(shù)值算例

3.1 自洽場(chǎng)算法分析

數(shù)值求解剛?cè)狍w系下，自洽場(chǎng)方程的關(guān)鍵在于初值的選取、求解傳播子方程以及鞍點(diǎn)問(wèn)題的迭代法．對(duì)于半剛性鏈求解已有一些高效算法[11]，但大多都針對(duì)周期邊界條件．為數(shù)值研究臨界膠束的需要，本研究采用反射邊界條件，并用平面直角坐標(biāo)、柱坐標(biāo)和球坐標(biāo)將空間的計(jì)算維度降到一維．此時(shí)，柔性鏈的傳播子方程對(duì)應(yīng)的算子記為

(30)

其中，n=0,1,2分別對(duì)應(yīng)直角坐標(biāo)系、柱坐標(biāo)及球坐標(biāo)系．對(duì)于柔性鏈的傳播子方程中，算子直接使用緊致差分格式半離散，再采用Runge-Kutta格式．對(duì)于剛性鏈的傳播子對(duì)應(yīng)的算子作如下分裂

LB=L1+L2+LCS

(31)

L1=a?x-wB(x)+

(32)

(33)

其中，LCS為不同坐標(biāo)系下引入的空間曲率算子．

分裂后的算子使用不同算法：在空間維度的擴(kuò)散方程L1， 為保證能量守恒以及算法的靈活性，采用緊致離散格式將傳播子方程半離散，得到常微分方程組再采用Runge-Kutta格式；而指向方向的擴(kuò)散方程L2和LCS則采用基于球諧函數(shù)的擬譜方法[13]．

3.2 傳播子方程的數(shù)值格式

本節(jié)介紹柔性傳播子方程及剛性傳播子方程的空間維度和指向維度的數(shù)值格式．柔性傳播子方程可統(tǒng)一記為

(34)

對(duì)于a(x)為常數(shù)的情形，構(gòu)造柔性傳播子方程的三點(diǎn)四階緊致格式為

(35)

代入式(34)可得四階半離散格式為

(36)

o(Δx3)

(37)

o(Δx3)

(38)

同樣可得求解剛性傳播子方程(空間維度)的四階半離散格式

(39)

這里的邊界條件采用如下近似

o(Δx3)

(40)

o(Δx3)

(41)

以上便是柔性和剛性傳播子方程(空間維度)的半離散格式，即線性常系數(shù)常微分方程組，它們可以統(tǒng)一寫(xiě)為

Ut=(AB-1C-W)U

(42)

對(duì)于式(42)的求解使用隱式三階Runge-Kutta格式．以下構(gòu)造剛性傳播子方程(指向維度)的數(shù)值格式，其傳播子方程統(tǒng)一記為

(43)

(44)

然后利用球諧函數(shù)性質(zhì)將傳播子方程(43)的第1項(xiàng)化為常微分方程組，且這個(gè)方程組的解可以記為

(45)

這樣，將式(45)代入式(44)就可得第1項(xiàng)的解q(u,s)， 這里使用SPHEREPACK軟件包來(lái)加速．同樣地，式(43)中的第2項(xiàng)也利用球諧函數(shù)展開(kāi)，得到相應(yīng)的線性常微分方程組為

Q′=ATQ

(46)

(47)

并對(duì)n做截?cái)鄟?lái)近似求解．至此，對(duì)于求解柔性和剛性傳播子方程的數(shù)值格式都已給出．

3.3 序參量計(jì)算和場(chǎng)迭代格式

序參量計(jì)算主要涉及的是數(shù)值積分，因此，對(duì)序參量的空間r和時(shí)間t積分時(shí)，采用四階復(fù)合數(shù)值積分公式

(48)

在對(duì)指向u采用Gauss積分公式和復(fù)合梯形公式．相應(yīng)地，在求解自洽場(chǎng)方程組時(shí)，需要對(duì)場(chǎng)進(jìn)行迭代，為加速解的收斂速度和減少解的誤差，先使用Picard迭代，當(dāng)場(chǎng)誤差達(dá)到一定范圍時(shí)，再啟用Anderson迭代．具體迭代格式為

ξ(n)+ψ(n)Gε(r-r1)]

(49)

ξ(n)+ψ(n)Gε(r-r1)]

(50)

M(n+1)=(1+αM)M(n)+αMηNS(n)

(51)

其中，αi(i∈A,B,M)為更新系數(shù)，通常取為0.01～0.10．

因此，序參量和自洽場(chǎng)模型數(shù)值求解的數(shù)值格式也已給出，整個(gè)自洽場(chǎng)模型數(shù)值求解步驟如下．

步驟1給出場(chǎng)函數(shù)wA、wB及M的初始猜測(cè)值；

步驟3通過(guò)迭代法更新場(chǎng)函數(shù)wA、wB及M；

步驟4若更新的場(chǎng)函數(shù)wA、wB及M是自洽的，且收斂誤差Etotal<ε， 則停止，得出數(shù)值結(jié)果；否則，返回步驟2．Etotal定義為自洽方程組(25)～(29)的l2范數(shù)，

(52)

場(chǎng)函數(shù)的初值選取過(guò)程為：① 計(jì)算出無(wú)序相對(duì)應(yīng)的序參量，猜測(cè)膠束結(jié)構(gòu)的厚度、剛性分子鏈的指向參數(shù)和場(chǎng)函數(shù)；② 設(shè)置計(jì)算區(qū)域、分片常數(shù)函數(shù)以及其他場(chǎng)函數(shù)；③ 利用移動(dòng)平均方法對(duì)場(chǎng)函數(shù)進(jìn)行多次光滑化，即得初值場(chǎng)函數(shù)．

3.4 數(shù)值算例

將上述自洽場(chǎng)算法用于計(jì)算臨界膠束結(jié)構(gòu)．在嵌段共聚物和均聚物體系中，隨著嵌段共聚物濃度的不斷增加，達(dá)到膠束的臨界濃度后，將自組裝形成膠束結(jié)構(gòu)[21]．至于形成何種膠束結(jié)構(gòu)，則取決于分子參數(shù)和相互作用力參數(shù)．為確定膠束的結(jié)構(gòu)，即臨界形態(tài)，需要考查不同膠束結(jié)構(gòu)的自由能隨濃度的變化．當(dāng)濃度增加時(shí)，自由能遞減，自由能最小時(shí)即為穩(wěn)定的膠束，相應(yīng)濃度稱為臨界膠束濃度．在巨正則系綜下，嵌段共聚物的濃度由其化學(xué)勢(shì)確定，

μc=lnzc=lnφbulk-ln(1-φbulk)+

(53)

其中，φbulk表示本體相(即自洽場(chǎng)方程組常數(shù)解對(duì)應(yīng)的相)中rod-coil共聚物的濃度．因此，巨正則系綜下的臨界膠束濃度即為臨界化學(xué)勢(shì)．

下面給出數(shù)值結(jié)果．選取空間的離散點(diǎn)為350～400個(gè)，鏈長(zhǎng)的離散點(diǎn)為300～800個(gè)，指向的離散點(diǎn)為32×36個(gè)，球諧函數(shù)的基函數(shù)個(gè)數(shù)為42～162個(gè)．結(jié)果表明，自由能、場(chǎng)函數(shù)及單鏈單鏈配分函數(shù)的誤差小于1×10-4，如圖2．

圖2 自由能、場(chǎng)函數(shù)和單鏈配分函數(shù)的誤差Fig.2 (Color online) Iterative errors of free energy, field functions and single chain partition functions

圖3 球狀、柱狀和平面膠束的臨界化學(xué)勢(shì)Fig.3 Critical chemical potentials of spherical, columnar and planar micelles

圖4 平面膠束在徑向方向的序參量Fig.4 Order parameters of planar micelles

圖5 柱狀膠束在徑向方向的序參量Fig.5 Order parameters of cylindrical micelles

圖6 球狀膠束在徑向方向的序參量Fig.6 Order parameters of spherical micelles

結(jié) 語(yǔ)

本研究基于描述帶剛性的wormlike鏈高分子模型，給出rod-coil/coil混合高分子體系在巨正則系綜下的自洽場(chǎng)模型及相應(yīng)算法．將自洽場(chǎng)算法運(yùn)用于rod-coil/coil體系的自組裝膠束，數(shù)值模擬結(jié)果得到3種膠束結(jié)構(gòu)的形成與轉(zhuǎn)變特性，與coil-coil/coil體系下形成自組裝臨界膠束的結(jié)構(gòu)與性質(zhì)類似．球狀膠束結(jié)構(gòu)穩(wěn)定性最高，平面膠束穩(wěn)定性最低，柱狀膠束位于兩者之間，且隨著嵌段的體積分?jǐn)?shù)增加，平面膠束會(huì)轉(zhuǎn)化為柱狀膠束再到球狀膠束.

致 謝：衷心感謝新加坡國(guó)立大學(xué)蔡永強(qiáng)博士的悉心指導(dǎo)！