基于分子晶體序參數(shù)與K-means聚類的TNT晶型轉(zhuǎn)化有限溫度弦研究

常玲玲,任福德,劉英哲,葛忠學(xué),王曉磊,邱麗莉,孟子暉,王艷紅,曹端林

(1.中北大學(xué) 化學(xué)與化工學(xué)院,山西 太原 030051; 2.西安近代化學(xué)研究所,陜西 西安 710065; 3.北京理工大學(xué) 化學(xué)與化工學(xué)院,北京 100081)

引 言

轉(zhuǎn)晶是物質(zhì)的基本相變過程,成核是轉(zhuǎn)晶的關(guān)鍵步驟,稀有事件方法是揭示成核本質(zhì)的重要途徑。目前,對于分子晶體轉(zhuǎn)晶的研究[1],稀有事件方法仍具挑戰(zhàn)性。一方面,分子晶體的分子結(jié)構(gòu)復(fù)雜、不同晶型的結(jié)構(gòu)相似,構(gòu)建能夠描述分子晶體成核的多維空間集合變量 (collective variable, CV)[2]極其困難;另一方面,利用增強采樣進行最低自由能路徑 (minimum free-energy path,MFEP)[3]的刻畫一直是稀有事件挑戰(zhàn)性的課題。由于分子晶體成核受尺寸效應(yīng)的影響,基于平均采樣的稀有事件方法致使自由能難以收斂。時至今日,稀有事件方法大多僅應(yīng)用于蛋白質(zhì)和核酸構(gòu)象的探索[4],而對于具有明顯稀有事件特征的分子晶體成核過程鮮有報道。

有限溫度弦(finite temperature string, FTS) 是一種基于路徑采樣的稀有事件方法[5],可通過對弦路徑的重新參數(shù)化,借助集合變量弦方法(string method in collective variables,SMCV)實現(xiàn)MFEP的快速演化?；趯Ψ植己瘮?shù)參數(shù)化模型的分子晶體序參數(shù) (order parameters, OPs) 作為CV已成功用于晶型的確定[6]。K-means聚類是一種將包含n維向量的數(shù)據(jù)集聚類迭代為k個子聚類的增強采樣算法,通過對分子晶體OPs的K-means聚類,利用FTS方法有望獲得MFEP快速收斂,實現(xiàn)分子晶體結(jié)構(gòu)演化與分子運動的統(tǒng)一,從而在分子水平上揭示其晶型轉(zhuǎn)化的本質(zhì)。

TNT(2,4,6-Trinitrotoluene)是應(yīng)用最廣泛的熔鑄炸藥載體,在分子水平上對其相變過程的剖析有助于人們對熔鑄過程本質(zhì)的認(rèn)識。本研究構(gòu)建基于TNT分子晶體的對分布函數(shù)OPs,結(jié)合基于Euclidean距離和權(quán)重分布的K-means聚類方法構(gòu)建平滑F(xiàn)TS,通過SMCV方法獲得TNT分子晶體轉(zhuǎn)晶的MFEP,揭示轉(zhuǎn)晶機制,提出了一種在分子水平上探索炸藥分子晶體轉(zhuǎn)晶的有效途徑,以期為高能鈍感材料轉(zhuǎn)晶工藝的優(yōu)化提供理論依據(jù)。

1 理論與方法

1.1 序參數(shù)構(gòu)建

TNT有正交(orthorhombic)和單斜(monoclinic)兩種晶型,分別記為TNT(O)(CCDC:227800)和TNT(M)(CCDC:227799)[7]。TNT(O)的晶型參數(shù)為:a=14.910(2)?,b=6.034(2)?,c=19.680(4)?,β=90.0°,Pca21;TNT(M)的晶型參數(shù)為:a=14.9113(1)?,b=6.0340(1)?,c=20.8815(3)?,β=110.365(1)°,P21/a。

根據(jù)實驗晶體結(jié)構(gòu)[8]構(gòu)建6×6×6 (216個分子)的TNT(O)和TNT(M)超晶胞,如圖1所示,圖中紅色、藍(lán)色、灰色和白色小球分別表示O、N、C和H原子。

圖1 6×6×6 (216個分子)的TNT(O)和TNT(M)的重要晶面Fig.1 Important lattice planes of TNT(O) and TNT(M) with the crystal size of 6×6×6 (216 molecules)

圖2 基于對分布函數(shù)的序參數(shù)構(gòu)建Fig.2 Construction of the order parameter based on pair distribution functions

(1)

(2)

(3)

基于鍵序參數(shù)的鍵取向和相對取向分別由下式給出 (在后面討論中,分別簡稱“鍵取向序參數(shù)”和“相對取向序參數(shù)”):

(4)

(5)

理論研究表明[11],局部序參數(shù)是描述成核反應(yīng)坐標(biāo)重要的集合變量。因此,在本研究中采用了兩類局部序參數(shù)作為集合變量。一類是針對特定分子i的序參數(shù),定義為體系中i的所有α-峰序參數(shù)之和,見式(6);另一類是針對體系中單元C的序參數(shù),定義為單元C序參數(shù)的平均,見式(7):

(6)

(7)

式中:*分別表示公式(1)～(5)中的d、b、r、db和dr。

1.2 弦方法

弦方法是一種通過評估反應(yīng)路徑上集合變量演化的局部平均力勢及其張量進行最小能量路徑或MFEP的一種技術(shù)[12]。給定一組用于描述反應(yīng)的N個集合變量,θ(x)=(θ1(x),…,θN(x)),x∈Rn為笛卡爾坐標(biāo),以z=z1,z2,…,zN表示的自由能定義為:

(8)

使t時刻弦z(α,t)(α∈[0,1]為弦的演化參數(shù))[3]演化收斂到MFEP的方程為:

(9)

給定自由能位于za和zb的兩個最小值,根據(jù)邊界條件z(0,t)=za和z(1,t)=zb求解方程(9)。當(dāng)t→∞,式(9)的解收斂于連接za和zb的MFEP。同時,在每個點z處的切向量平行于M(z)F(z), 即:

(10)

(11)

(12)

(13)

(14)

(15)

同理,根據(jù)式(11)可求得:

(16)

1.3 FTS的實現(xiàn)

1.3.1 初始軌跡的構(gòu)建

1.3.2 K-means聚類與初始弦的平滑

眾所周知,String方法在多維情況下極易陷入局部極優(yōu),這是其理論所決定的(在確定初始string的情況下,無法克服與初始string正交空間的能壘)。為盡可能防止其陷入局部極優(yōu),進行初始弦的平滑,采用K-means聚類算法進行采樣,包括以下4個步驟:

步驟1:獨立空間限制性采樣。在每個給定代表點(副本)周圍的獨立空間中根據(jù)如下諧波函數(shù)采樣:

(17)

(18)

式中:kd、kb、kr分別對應(yīng)距離、鍵取向和相對取向spring系數(shù),其值分別為20910kJ/mol、4182kJ/mol和4182kJ/mol。

通過這種方式,每個模擬都位于由每個代表點副本確定的獨立空間中,在每個獨立空間進行獨立模擬(3×10ns)。

步驟2:K-means聚類。給出兩個假設(shè):(1)在每個獨立空間中存在具有最大權(quán)重的聚類中心;(2)不同空間中最大權(quán)重聚類中心對應(yīng)的序參數(shù)不同。顯然,這樣的假設(shè)得到的聚類中心可被視為局部重要采樣點,可顯示弦代表點的演化方向。

對于3×10ns的模擬軌跡進行K-means聚類。將每個獨立空間的軌跡分為15個簇,借助維度加權(quán)Euclidean距離測量樣本點之間的距離,在求得所有樣本密度權(quán)重之后,在每個獨立的采樣空間中選擇權(quán)重最高的點作為弦的代表點。權(quán)重ωi由下式給出:

(19)

式中:ρi、Si和ai分別表示樣本點i的密度、不同數(shù)據(jù)集之間的距離以及同一聚類中不同樣本點之間的距離,分別通過公式(20)～(22)計算得到:

(20)

(21)

(22)

顯然,基于密度的方法在理論上更為自然,這就是本研究將樣本點i的密度引入聚類算法的原因。然后,連接最大權(quán)重聚類中心形成臨時弦。

步驟3:平滑函數(shù)搜索連續(xù)弦。為了獲得集合坐標(biāo)空間中連續(xù)弦,定義平滑函數(shù):

(23)

通過蒙特卡洛方法計算臨時弦的平滑度。首先,計算當(dāng)前臨時弦的平滑分?jǐn)?shù)。然后,在同一空間中選擇與第一個聚類中心的權(quán)重差較小的另一個聚類中心,重新計算平滑分?jǐn)?shù),根據(jù)Metropolis準(zhǔn)則決定是否接受新的聚類中心。最后,連接新聚類中心得到平滑分?jǐn)?shù)低的初始弦。

1.3.3 MFEP的刻畫

(24)

1.4 MD模擬

在300K和1標(biāo)準(zhǔn)大氣壓下,依據(jù)實驗晶體結(jié)構(gòu),對216個分子TNT(O)和TNT(M)進行1.0ns MD模擬。借助particle mesh Ewald方法對周期性靜電相互作用進行建模,引進振蕩頻率為25ps-1的Langevin熱浴控溫。平衡后,在恒定體積下退火至零溫度,進行能量最小化弛豫。然后,升溫至300K,在恒定體積下進行2.5ns MD模擬,時間步長0.5fs,根據(jù)軌跡借助最大似然估計方法計算式(1)～(7)的參數(shù),確定初始序參數(shù)。

對于FTS演化的MD模擬,系綜仍為NPT(300K和1標(biāo)準(zhǔn)大氣壓)。每一步初始構(gòu)型均選擇與前一步最接近新目標(biāo)序參數(shù)的構(gòu)型。為保持模擬的穩(wěn)定性,在前一段時間內(nèi)序參數(shù)約束的中心從上一步的相應(yīng)值移動到新值;在后一段時間內(nèi)平均每步的約束力,且初始速度均從所處溫度的Maxwell-Boltzmann分布隨機分配。CHARMM22力場[14]用于苯的熔化,計算結(jié)果與實驗值吻合[9]。所有模擬借助自編程序和NAMD軟件包[15]的CHARMM22力場進行,利用PLUMED軟件包進行了序參數(shù)和自由能計算[16]。模擬流程見圖3。

圖3 MD模擬流程圖Fig.3 MD simulation process

2 結(jié)果與討論

2.1 TNT(M)和TNT(O)對分布函數(shù)峰值

TNT(M)和TNT(O)的對分布函數(shù)峰值在355K下獲得 (實驗轉(zhuǎn)晶溫度為351K),截至值為7.8?。平均峰位置和集中參數(shù)見表1(括號內(nèi)為355K溫度下計算得到的參考結(jié)構(gòu)對分布函數(shù)峰值)。與參考結(jié)構(gòu)(括號內(nèi))相比,TNT(M)和TNT(O)的對分布函數(shù)峰值變化不大;對于兩種晶型,對分布函數(shù)的峰值非常接近。兩種晶型結(jié)構(gòu)的高度相似性,不僅使構(gòu)建能夠精確描述其轉(zhuǎn)晶的多維空間集合變量極其困難,轉(zhuǎn)晶序參數(shù)發(fā)生“維數(shù)爆炸”,而且在弦的演化過程中,可能使自由能難以收斂。另外,較小的標(biāo)準(zhǔn)偏差和較大的集中參數(shù)顯示了兩種晶型平均質(zhì)心距離、鍵取向角和相對取向角分布的定域性特征。

表1 355K溫度下TNT(O)和TNT(M)平均峰位置及其集中參數(shù)Table 1 Average peak locations and concentration parameters for the TNT(O) and TNT(M) crystals at 355K

2.2 K-means聚類的收斂性

對于21個獨立空間得到的10ns×3軌跡,分別進行了鍵取向和相對取向序參數(shù)的平均采樣。結(jié)果發(fā)現(xiàn),對于鍵取向序參數(shù),ln(Convergence)值較大(-4.0～0.0),接近0.0的ln(Convergence)值可能是序參數(shù)“維數(shù)爆炸”引起的;對于相對取向序參數(shù),迭代次數(shù)達(dá)到150次仍未見收斂,如圖4所示(I、II和III分別表示經(jīng)K-means聚類采樣的鍵取向和相對取向序參數(shù)以及基于平均采樣的相對取向序參數(shù)收斂),這表明基于平均的采樣具有明顯的分散性和不確定性。在這種情況下,來自一個采樣空間中的模擬可能覆蓋沿過渡路徑的相鄰區(qū)域,使連接相鄰點的初始弦不平滑。眾所周知,對于非光滑的弦,總有一些獨立的采樣空間不垂直于當(dāng)前弦,這使得每個代表點周圍的受限樣本相互干擾并交織在一起,導(dǎo)致弦的演化方向混亂。

圖4 弦演化過程中集合變量的收斂 Fig.4 Convergence of the collective variables during the string evolution

對于FTS方法,沿路徑的每個代表點可被視為晶型轉(zhuǎn)換過程中的中間狀態(tài)。平滑的弦意味著這些中間狀態(tài)可以從一端到另一端連續(xù)變化,而使弦平滑的目的是消除不相關(guān)的信息。為了得到平滑的初始弦,本文采用維度加權(quán)Euclidean距離結(jié)合密度權(quán)重的K-means聚類算法分別對鍵取向和相對取向序參數(shù)進行采樣,發(fā)現(xiàn)經(jīng)K-means聚類增強采樣后,兩種序參數(shù)的弦分別經(jīng)14次和11次迭代后收斂。應(yīng)用聚類算法不僅可以找到權(quán)重最大的聚類中心,解決相鄰副本中高度相似序參數(shù)的區(qū)分,而且可以加快弦的收斂速度。一般地,維度加權(quán)Euclidean距離結(jié)合密度權(quán)重的K-means聚類算法具有較強的異常值干擾能力,在一定程度上,可有效避免局部最優(yōu),快速收斂到全局最優(yōu)。近年來,K-means聚類算法已廣泛應(yīng)用于DNA、RNA和蛋白質(zhì)序列的篩選和預(yù)測[17]。

2.3 最小自由能路徑

圖5 PMF隨弧長的變化曲線 Fig.5 Curves of PMF with arc length

由圖5可知,TNT(O)與過渡態(tài)之間的PMF差值分別為102.5和99.5kJ/mol。TNT(O)與TNT(M)之間PMF的差值分別為34.3和30.1kJ/mol。由于PMF對應(yīng)于27維(3×3×3)OP空間,而不是一維反應(yīng)坐標(biāo),PMF比真實自由能壘高。另外,PMF值與體系的大小以及Spring系數(shù)有關(guān)。

圖6 FTS路徑上局域序參數(shù)的變化 Fig.6 Changes in local order parameters on the FTS path

從圖6可以看出,沿著TNT(O)→(IA)→(IIA)→(IB)→(IIB)→(IC)→(IIC)→TNT(M)路徑,分子的類型從“完全的藍(lán)色”逐漸變?yōu)椤氨凰{(lán)色包裹的零星綠色”分子簇的生成(即“界面誘導(dǎo)”)、“零星綠色”分子簇的增大(即“局部引發(fā)、多核非同步生成與增大”)到“藍(lán)色被綠色吞噬”,最后形成“完全的綠色”,即“由TNT(O)向TNT(M)晶型的完全轉(zhuǎn)化”。這表明,TNT(O)與TNT(M)晶型轉(zhuǎn)變是界面誘導(dǎo)、局部引發(fā)、多核非同步發(fā)展的過程,這與Trout等對分子晶體晶型轉(zhuǎn)變的研究結(jié)果吻合[19-20]。

圖7 沿FTS路徑的PMF隨全局序參數(shù)變化的3D圖 Fig.7 3D diagram of PMF along FTS path with global order parameters

3 結(jié) 論

(1)基于Euclidean距離和密度權(quán)重的K-means聚類算法可以用于TNT炸藥分子晶體轉(zhuǎn)晶序參數(shù)的增強采樣,可以改進稀有事件常規(guī)有限溫度弦方法,使自由能快速收斂。

(2)基于分子晶體序參數(shù)與K-means聚類有限溫度弦方法可用于TNT炸藥晶型轉(zhuǎn)化自由能的計算,因而可以揭示炸藥分子晶體晶型轉(zhuǎn)化稀有事件的本質(zhì)。

(3)TNT(O)與TNT(M)晶型轉(zhuǎn)變是界面誘導(dǎo)、局部引發(fā)、多核非同步發(fā)展的過程。