V01.36 2015年11月 高等学校化学学报 CHEMICAL JOURNAL OF CHINESE UNIVERSITIES No.11 2179~2188 耗散粒子动力学与可极化ABEEM力场结合 模拟嵌段共聚物自组装形态 刘林林,杨忠志 (辽宁师范大学化学化工学院,大连116029) 摘要本文采用介观尺度上的耗散粒子动力学与基于ABEEM可极化力场的全原子分子动力学相结合的方 式,从理论上探究了嵌段共聚物结构、对称性、分子组成及温度等对嵌段共聚物最终自组装结构形态的影 响.模拟结果表明,这些因素均会对嵌段共聚物在选择性溶剂中的最终自组装结构产生一定影响.这一理论 研究为实现可控操纵自组装结构,提供了有意义的参考. 关键词耗散粒子动力学;ABEEM力场;嵌段共聚物;囊泡;胶束;自组装结构 O641 文献标志码A 中图分类号嵌段共聚物是由2种或多种具有不同化学结构的单体序列以共价键相连接而形成的聚合物.嵌段 共聚物有多种分类方法:按照结构可将其分为线型共聚物、星型共聚物、梳型共聚物、H型共聚物和 环形共聚物等…;按照组成嵌段的数目可分为两嵌段共聚物、三嵌段共聚物以及多嵌段共聚物等[2]. 由于不同嵌段具有不同的物理和化学性质,因此嵌段共聚物展现出了区别于均聚物的特殊优异的性 能.嵌段共聚物的一个显著性质是其能够在体相或者溶液中自组装成多种有序结构,如球状胶束、棒 状胶束、层状结构、管状结构和囊泡结构等 .这些胶束和囊泡等结构可用于制造微反应器及药物载 体等,因此具有广泛的用途和商业前景,受到越来越多的重视和研究 j.另外,两亲性嵌段共聚物自 组装成的囊泡在结构上与细胞膜极为相似,两者均具有双分子层的膜状结构,结构中间的空腔可以包 含溶剂分子和药物分子,因此这种囊泡结构具有很强的仿生意义 . 近年来,嵌段共聚物在选择性溶剂中自组装的潜在应用前景广阔,已成为实验和理论研究的热 点.研究表明,嵌段共聚物的分子量、嵌段的顺序、分子结构_6 J、亲疏溶剂嵌段的相对长度、溶剂选择 性、溶液pH、离子强度和温度等因素都会对自组装行为产生影响 .因此对这些影响因素进行分析, 从而对自组装结构进行定向,具有至关重要的意义. 目前,对于软物质体系,如聚合物溶液、熔体、胶体溶液及多种生物体系,采用透射电子显微镜、 原子力显微镜和扫描电子显微镜等实验技术手段对其进行研究.但这并不能充分提供微观结构信息及 动态性质.因此,需要理论模拟来探究其中的规律和机理,揭示介观结构与宏观性质之间的关系.考 虑到支配这些软物质的性质所在的尺度,并不完全是微观细节所决定的,而是发生在介观的长度尺度 和时间尺度上,而这远在原子尺度模拟方法所能处理的范围之外.针对这种情况,亟需从一种更为宏 大的视角来看待这一问题.因此在实际情况下,必须设计方法来对基础体系进行尽可能的简化,只保 留控制所感兴趣的关键性质.粗粒化动力学模拟是研究大空间尺度和长时间尺度下生物体系的一种实 际方法.粗粒化模拟采用将原子尺度的结构和相互作用进行简化来提高计算效率,实现对原子尺度的 分子动力学模拟所无法探究的过程进行研究. 最近已有多种粗粒化方法被开发,并用于对多种过程进行理论模拟.其中,耗散粒子动力学在聚 合物的粗粒化模拟中得到了较为广泛的使用,并且与实验结果取得了很好的一致性.尽管耗散粒子动 力学 中所使用的粗粒化模型在对一些现象的理解上已经变得很强大,但其预测力仍主要是定性的, 收稿日期:2015-07—17.网络出版日期:2015.10—10. 基金项目:国家自然科学基金(批准号:21133005,21073083)资助. 联系人简介:杨忠志,男,博士,教授,博士生导师,主要从事理论与计算化学研究.E-mail:zzyang@lnnu.edu.cn 高等学校化学学报 Vo1.36 而非定量的.事实上,使用无维度的单位来衡量长度、时间和能量尺度,当与实际的体系进行关联时, 这种物理量上的投影通常是事后进行的,导致模型参数和实际的分子参数之间缺乏直接的关联.典型 的处理方式是使用一些可观测量作为标准设定模拟参数 J.同时,基于要保留的现象来进行模拟参数 的发展还需要相当大的努力 .研究¨ 指出,在微观尺度与介观尺度之间寻找对应关系是一项极具 挑战性的任务.然而,一旦这种对应关系被发现,那么就可以从分子动力学模拟中推导出耗散粒子动 力学中所使用的相互作用参数,而不再是以宏观可测物理量作为参考进行参数的设定.1997年,Groot 和Warren_l 在这方面做出了重要的贡献,建立了耗散粒子动力学中保守排斥的简单泛函与 Flory.Huggins参数之间的关系.同时在聚合物理论方面,也指出Flory—Huggins参数与聚合物的溶解参 数及混合能之间也存在关联.而聚合物的溶解度参数和混合能是可以从原子尺度的模拟或者实验中获 得的.目前通过实验获得的溶解度参数仍很有限 .因此,将介观尺度的耗散粒子动力学与全原子分 子动力学模拟相结合尤为必要.Groot和Warren 为耗散粒子动力学中的保守排斥参数可以从原子尺 度模拟推导获得,提供了合理的理论基础.但存在的问题是:一种粗粒化模型如何实现事先从全原子 模型进行参数化,以便进行实际体系的模拟,然后用于定量预测;另外这种方式能否很好地重现实验 上所测得的溶解度参数及自组装形态.在之前的研究 中,我们进行了一些尝试,并在聚合物的组装 方面取得了一定的阶段性成果¨ .但是全原子分子动力学模拟均基于非极化力场,而非极化力场中非 键相互作用的准确描述是以牺牲非键参数的可转移性为代价的.因此,将耗散粒子动力学与基于可极 化力场的全原子分子动力学相结合非常必要.由于内聚能及溶解度参数敏感地依赖于所使用的力 场lL】 ,这样全原子分子动力学模拟就必须采用一种高质量的力场. 聚苯乙烯一聚异戊二烯(Ps.PI)嵌段共聚物存在多种多样的结构 ,其在聚异戊二烯的选择性溶剂 正庚烷中可以自组装成多种结构,并且已经得到实验证实.当前的理论研究主要局限于其在体相中的 自组装,而对其在选择性溶剂中的自组装则很少有研究.溶剂的引入将会导致自组装形态的复杂化, 同时也将获得具有新奇性质的新型纳米结构.因此对嵌段共聚物在选择性溶剂中自组装的深入研究非 常必要.对于PS—P1,其最终的自组装形态受多种因素影响,如嵌段共聚物的结构¨ 、对称性 、分 子组成 。。及温度 等,并已得到实验证实.然而,对选择性溶剂中该类嵌段共聚物自组装机制及调 控因素的理论探究仍然欠缺. 我们采用将介观尺度的耗散粒子动力学与基于ABEEM可极化力场 的全原子分子动力学模拟 相结合的方法,对PS—PI嵌段共聚物在正庚烷中的自组装形态及其因素开展了一系列模拟研究. ABEEM可极化力场 。 呈现出极好的性能和良好的参数可转移性,其对体系中非键相互作用有很好的 描述,尤其是对静电相互作用的描述具有独特的优势 . 1理论模型与方法 1.1基础理论 耗散粒子动力学(DPD)是Hoogerbrugge和Koelman 引在1992年引入的一种介观尺度上的模拟方 法,主要用于复杂流体的模拟研究 .它实现了在更大的长度和时间尺度下对体系进行模拟研究.耗 散粒子动力学中的珠子所表示的是一组原子或者一定体积的流体,并按照化学特性和环境对其进行分 类.这种珠子在原子尺度上是很宏大的,但是在宏观尺度上依然是很微小的,这样就实现了模拟尺度 上几个数量级的增加.其中,所有珠子的位置和速度随时间的演变遵从牛顿运动方程.在一定截断半 径内,作用于粒子上的作用力可以分为5个部分: F =∑(F +F D+F R)+F + (1) 式中: 表示保守力; 表示耗散力,是移动的珠子之间的黏性曳力,起耗散能量的作用; 表示随 机力,代表随机排斥,起补偿能量的作用;F 和F:共同作用充当珠子的热浴;以上所有的力都是短程 的,即在一定的截断半径内起作用; 表示弹簧力; 表示键角力.保守力 是软排斥力,作用于连 接2个珠子中心的连线上: 刘林林等:耗散粒子动力学与可极化ABEEM力场结合模拟嵌段共聚物自组装形态 C—f口 (1一rSrc (ro.<r ) , 一1.0 (r ≥≥r ) ‘ 式中: = 一 ,r = l,e =rij/r ,而,. 和0是作用于粒子i和 上的位置矢量;r 是截断半径, 是体系中唯一的长度尺度,因此使用截断半径作为长度单位, :1.0.软排斥力的使用,允许在更大 的长度和时问尺度下进行模拟.n 是2个相互作用珠子之间的最大排斥相互作用参数,依赖于混合物 的特性,在确定聚合物溶液的平衡纳米结构中起重要作用. 对速度依赖的耗散力对应于黏性曳力,作用于粒子的相对速度,降低了粒子之间的速度差,依赖 于粒子的位置和相对速度.而随机力则弥补了耗散力所导致的动能损失.耗散力和随机力共同充当了 体系的热浴,并且使得模拟中的总动量得以守恒.它们的具体形式如下: FD =一TO) ( )(', f・e )eif (3) F =or(o (rij) At-1/2e (4) 式中: =l, —l,,,l, 和l, 为粒子的速度矢量; 为噪声强度,可控制每个时间步内由随机力输人到体 系中的能量; 为摩擦系数;O9。(r )和∞ (r )为依赖于距离的权重函数,分别描述耗散力和随机力随 粒子间距离的增加而衰减的情况; 为随机涨落变量,呈高斯几率分布,平均值为0,平方差为1.在 每一时间步下,对每对相互作用的珠子而言,这一数值是选择的. ( f(t))=0,( f( ) (t ))=(6 +6 z )6(t—t ) (5) 为了满足涨落耗散定理,以及使得体系演变到对应于吉布斯系综的平衡状态,需要满足下式: 山。( )=[ ( )] (6) 按照Groot和Warrenl1 的设定方法,在符合涨落耗散定理的条件下,为了满足正则系综的条件,对于 水溶液中的软链而言,在一定的截断半径内要满足 (, f):[ (,ff)z]:f 一 / c) L0 <rc) ( ≥r ) (7) 另有, =2yk T,式中: 为绝对温度;k 为玻尔兹曼常数;7和or通常被设定为4.5和3.0_1 ,保 守相互作用势k T被选定为1.0. 弹簧力将相邻的珠子相连构成聚合物链,其形式如下: F =∑Cr 有略微小的距离.而键角力则常用于包含刚性嵌段的共聚物研究中,以表征刚性嵌段的刚性. (8) 式中:c为高斯链中2个邻近珠子之间的弹簧常数,选定为4.0,这种设定会使得成键粒子比非键粒子 在耗散粒子动力学中,珠子表示真实原子或聚合物单体的集合,导致原子水平细节的缺失,不同 类型珠子之间的相互作用参数口 成为控制介观体系行为的主要参数,对重现嵌段共聚物自组装形态 起着至关重要的作用.耗散粒子动力学理论与聚合物的Flory—Huggins理论具有一定联系.耗散粒子动 力学中的相互作用参数0 可以投影到Flory—Huggins参数 具体地,为了获得 ,首先要确定口 和 这2个相同类型珠子之间的排斥相互作用参数可通过下式获得: 口 =75kBT/p (9) 式中:P为珠密度,通常选择为3;而保守相互作用势 仍设定为1.实际上,珠密度是一个可以自由 选择的参数,但若选择一种较大的密度,则要花费较长的计算时间.按照Groot和Warren L1 的研究, 对于液体而言P=3和n=25k T是一种较为合适的参数设置. 在P=3的情况下,不同类型珠子之间的排斥参数0 是按照如下的线性关系进行确定的: 8 =口 +3.27xi (1O) 式中:X j是从溶解度参数计算出来的,其依赖于珠子的化学特性,实验上可以通过混合能法 。“和溶饵 度参数法 ¨获得,但是实验上可获得的数据是有限的.本文从基于ABEEM可极化力场 的全原子分 子动力学模拟来获得这一参数. 力场选择的决定性因素是所研究体系性质的体现,其参数化是一项重要和艰巨的工作.当一种新 高等学校化学学报 Vol_36 的力场或者是原始的力场之外引入新参数用于计算新性质时,需要进行参数的拟合,以便获得良好的 一致性.通常,一种力场对几个参数是非常敏感的,而这些通常是用于计算非键相互作用和二面角相 互作用的参数.实际上,非键相互作用描述的是否准确是衡量一种力场准确与否的标准.对极化作用 进行明确处理,代表了生物分子力场中对非键相互作用处理上的重要进展 .在这方面,ABEEM可 极化力场 对体系中的非键相互作用给出了很好的描述,尤其是静电相互作用的处理上,更凸显其特 色和优势.事实上,分子间相互作用不仅来自于范德华相互作用和静电相互作用,同时键长和键角的 改变也有一定贡献.相对于范德华相互作用和静电相互作用在非键相互作用中所体现的重要性,这种 贡献并不是主要的.对于ABEEM可极化力场而言,分子问非键相互作用的计算是针对所有位点,包 括原子、键和孤对电子加和获得的 . 一个体系能量中的非键部分,通常采用成对的分子内和分子间相互作用的库仑和Lennard—Jones贡 描述的是非键的原子之间的相互作用,具有下面的形式: 献加和计算出来. E d =∑4 s ( 12/r 12一 6 /r6 ) i<J (11) 对于Lennard—Jones系数,采用如下的几何组合规则: f=( ) 。, f=(gii ̄ ̄ ) .对所有的非键 原子对i和J.进行加和.如果i和 在同一分子,如成键1~2原子对和成键角的1~3原子对之间, = 0;而对于组成二面角的1~4原子对之间,. =0.5.对于分子内和分子间相互作用, =1.0.与其它 力场相比,ABEEM可极化力场的显著特色是在静电相互作用的处理上 .静电相互作用通常表示为 /5…1 (J kijq ̄qj/rij l l2 式中:r 是位点i和 之间的距离;q 和qj是2个位点的电荷,通过ABEEM计算获得 ;k 是考虑了交 换、渗透以及屏蔽效应的参数; 的数值通常设定为0.57,但对氢键相互作用区域需要进行特殊处理. 在模拟中,ABEEM被用于计算原子、 键、仃键和孤对电子区域的电荷及体系的总能量 .如果键 长、键角、二面角或分子之间的相对位置发生改变时,将重新进行电荷和总能量的计算.由此可以获 得与分子环境变化相一致的部分电荷 ]. 对静电相互作用能的计算通常使用固定的以原子为中心的单极确定.这种表示的不变性,阻止了 对外界静电场变化的响应,因此缺失了对电子极化的准确描述.极化导致额外的分子间的吸引作用 力.对于缺失极化效应所导致的对分子间相互作用描述上的误差,可以通过对色散相互作用以及原子 的部分电荷的高估来补偿.这种补偿导致非可极化力场中从气相二聚体高水平的量子化学计算所获得 的非键参数,并不能简单地转移到凝聚相体系.因此,明确处理极化,将对从头算所获得的气相参数 转移到大尺度的生物化学体系是尤为重要的.因此,基于ABEEM可极化力场 的全原子分子动力学 模拟可以用来计算内聚能以及溶解度参数. 耗散粒子动力学模拟的关键是获得珠子之间的相互作用参数.而这种相互作用参数与珠子之间的 混溶性是紧密相关的.溶解度参数6具有温度依赖性,广泛用于对溶解行为的粗略和近似描述.溶解 度参数是由内聚能计算获得的,而内聚能表示的是孤立链体系的非键能与周期性边界条件下的无定型 元胞中计算的非键能二者的差值 ,即 Ee0h=[∑ Ei ̄o ! (i)一既】/n 式中: iSO ( )为真空中第i条孤立链的非键能;E n 为周期性边界条件下的 条链的非键能. (13) 溶解度参数6是内聚能密度的平方根.内聚能密度是将单位体积的分子从近邻处移动到无限远处 所需要的能量,它衡量的是分子间力的强度: 6=. (14) 式中: 是平衡时盒子的体积.Flory.Huggins参数 可以通过下面的方程从溶解度参数计算出来: . := ,iJB1 (6 一6 ) (15) 式中: 是聚合物片段的体积,即耗散粒子动力学中珠子的大小. 2188 2O 21 22 23 24 高等学校化学学报 V0l_36 Zhulina E.B.,Adam M.,LaRue I鹞 .,Shei∞ ko S.S.,Rubinstein M.,Macr”omol弘 ecules,2005,38(12),5330-∞ 躬* -5351盯 LaRue I.,Adam M.,Pitsikalis M.,Hadjichristidis N.,Rubinstein M.,Sheiko S.S.,Macromolecules,2006,39(1),309_314 Wang F.F.,Zhao D.X.,Yang Z.Z.,Chem. .,2009,360(1),141一l49 Wang F.F.,Zhao D.X.,Gong L.D.,Theor.Chem.Ace.,2009,124(1/2),139—15O Zhao D.X.,uu C.,Wang F.F.,Yu C.Y.,Gong L.D.,Liu S.B.,Yang Z.Z., Chem.Theory Comput.,2010,6(3),795— 8o4 Yang Z.Z.,Wang C.S., Phys.Chem.A,1997,101(35),6315--6321 Wang C.S.,Li S.M.,Yang Z.Z., Mo1.Struct:Theochem.,1998,430,191—199 Yang Z.Z.,Wang C.S., Theor.Comput.Chem.,2003,2(2),273_299 Hoogerbrugge P.J.,Koelman J.M.V.A.,Europhys.Lett.,1992,19(3),155—16O Li X.,Tang Y.H.,Liang H.,Karniadakis G.E.,Chem.Commun.,2014,50(61),8306---8308 Zhao Y.,You L.Y.,Lu Z.Y.,Sun C.C.,Polymer,2009,50(22),5333--5340 “Y.,Xu G.,Luan Y.,Yuan S.,Zhang Z.,Colloid Su A,2005,385,257_258 Yang Z.Z.,Wang J.J.,Zhao D.X., Compm.Chem.,2014,35(23),1690--1706 Chinese Universities,2014,35(12),2645--2653(刘翠,张千 Maekerell A.D.,J.Comput.Chem.,2004,25(13),1584—16O4 Liu C.,ZhangQ.H.,GongL.G.,LuL.N.,Yang Z.Z.,Chem. 慧,宫利东,卢丽男,杨忠志.高等学校化学学报,2014,35(12),2645--2653) Guan Q.M.,Yang Z.z., Theor.Comput.Chem.,2007,6(4),731.746 Ponder J.W.,Tinker,Sotware ToolfsforMolecularDesign,Venion 4.2,2004 Phillips C.L.,Anderson J.A.,Glotzer S.C.,J.Chem.Phys.,2011,230(19),7191 DuH.,Zhu J.,JiangW., Phys.Chem.B,2007,111(8),1938--1945 Femyhough C.M.。Chalari I.,Pispas S.,Hadjiehristidis N.,Eur.POtym.J.,2004,40(1),73 9 Mizerovskii L.N.,Smimova K.P.,Russ.Chem.曰+,2009,58(8),1547--1561 Guo H.,Qiux.,Zhou J., Chem.P .,2013,139(8),084907 Soto-Figueroa C.,Rodrfguez-Hidalgo M.D.R.,Martinez—Magaddn J.M.,Polymer,2005,46(18),7485--7493 Luo Z.,Jing aJ.,Polymer,2010,51(1),291—299 Takahashi T.。R&D Review ofToyota CRDL,2003,粥(4),1__9 Uu L.L.,Yang Z.Z.,Zhao D.X.,Gong L.D.,Lju C.,RSCAdv.,2014,4(94),52083--52087 Wells S.L.,The& ofAmphiphilic Block Copolymers in Selectie Solventvs,University of Noah Carolina,Chapel Hill,2000 Hart M.,Hang M.,Sim E., Chem.Phys.,2011,134(20),204901 Block Copolymers Simulated v/a Self-assembly Morphological Control of the Combination of Dissipative Particle Dynamics and Force Field* ABEEM Polarizable LIU Linlin,YANG Zhongzhi (School ofChemistry and Chemical Engineering,Liaoning Normal University,Dalian 116029,China) Abstract Mesoscopic dissipative particle dynamics(DPD)simulations and all—atom molecular dynamics simulations were combined to investigate the influence of the architecture,symmetry,moleculr composiation of block copolymers and temperature on final assembled structures of block copolymers.From the simulation resuhs,these factors have an influence on the final assembled structures.This theoretical investigation pro- vides a reference for experiments and practice to manipulate self-assembled structures. Keywords Dissipative particle dynamics;ABEEM force field;Block copolymer;Vesicle;Micelle; Self-assembly structure control (Ed.:V,Z) t Supposed by the National Natural Science Foundation of China(Nos.21 133005,21073083)