原子间相互作用对无序大分子...和动力学影响的中子散射研究_韩泽华.pdf

原子间相互作用对无序大分子结构和动力学影响的中子散射研究韩泽华1,2,左太森1,2,马长利1,2,李雨晴1,2,3,程贺1,2*1.中国科学院高能物理研究所,中国散裂中子源,东莞 5238002.散裂中子源科学中心,东莞 5238003.中国科学院大学,北京 100049*通讯作者,E-mail:收稿日期:2022-12-30;接受日期:2023-02-07;网络版发表日期:2023-03-31国家自然科学基金(编号:22273046,U1932161)资助项目摘要不同单体在原子尺度的相互作用使得不同大分子体系在微观、介观的多尺度复杂结构及动力学行为有明显差异,从而进一步影响了体系的宏观性质.联用中子散射和计算机模拟,利用相同分子量和分子量分布的氘代大分子与氢化大分子具有相同分子结构、不同中子衬度,以及中子散射和计算机模拟宽广的时间、空间观察尺度,我们可以得到无序大分子体系最可几全原子结构,进而分析其从原子到纳米的多尺度空间结构,与从皮秒到微秒的多模式动力学行为的成因.近年来,我们使用该方法,从小分子到大分子的稀溶液、溶胀体系,从大分子熔体到玻璃态,成功分析了原子间相互作用对不同空间尺度结构和跨时间运动模式动力学行为的影响.本文介绍了这些典型的例子,希望将该方法推广到更广阔的研究领域,把大分子原子结构的多样性与多尺度的复杂结构和动力学有机地联系起来.关键词无序大分子中子全散射,准弹性中子散射,计算机模拟,全原子结构,多尺度结构和动力学1引言因为高分子材料的性能由组成其的单体成分和加工过程中形成的多尺度结构共同决定,所以高分子领域的物理学家和化学家分别采用独特的方法来提高材料的性能.物理学家的解决方式是建立具有一般性的理论来描述体系中结构、动力学行为和宏观性质的关系,并将其应用到未知体系中预测性能和表现.比如,Flory1的平均场理论就建立在粗粒化的格子模型上,通过引入不可压缩假定,描述了高分子的基本相分离过程;de Gennes的标度理论同样建立在粗粒化的串滴(blob)模型上,描绘了温度blob和浓度blob的标度变化等2.而高分子化学家的解决方案则是合成由不同单体组成的新型高分子,利用不同单体间复杂的相互作用组装成各种各样的多尺度结构.在两种方式相互联系的诸多桥梁中,经验参数和经验公式发挥了重要的作用.例如,依据溶度参数预测高分子在不同溶剂中的溶解能力,从而选择不同溶引用格式:Han Z,Zuo T,Ma C,Li Y,Cheng H.The influence of interatomic interaction on the structure and dynamics of disordered macromolecules:a neutronscattering study.Sci Sin Chim,2023,53:678692,doi:10.1360/SSC-2022-0257 2023 中国科学杂志社中国科学:化学2023 年第 53 卷第 4 期:678 692SCIENTIA SINICA C聚合物结构与性能专刊专题论述剂使大分子自组装3,4;依据侧基的空间位阻来推测高分子膜的气体选择能力5;依据Vogel-Fulcher-Tam-mann(VFT)方程来计算高分子玻璃化转变过程的动力学特征弛豫时间等6.无法定量这些经验参数与经验公式的主要原因是研究手段的匮乏:传统的研究手段无法在高分子材料本体中实时观测其全原子结构及动力学,只能通过微观、介观、宏观尺度上的一些结构与性能或者光谱的变化推测多尺度结构成因的原子图像7,8,由于这种间接观察的方式不直观,所以研究者很难从单独个体特性出发,总结出系统性的结论.中子散射与计算机模拟相结合,可以有效地解决这个问题.一方面,二者都可以覆盖从埃到微米、从皮秒到微秒的多尺度范围;另一方面,二者又可以互补不足,互相支撑.中子的优点主要表现在:强穿透性,可以在复杂的样品环境下实现对高分子体系结构(变化)的实时观测;氘代技术,可以选择性观察复杂体系中不同组分的结构和动力学等9,10.中子散射的主要缺点是不直观:作为一种倒空间的系综平均,中子散射无法给出体系中微观过程的直观物理图像;而计算机模拟最大的优势就是直观,通过模拟结果可以直接确定体系实空间的结构与运动:通过使用全原子或者粗粒化力场对实验体系进行多尺度模拟,直接与中子散射实验结果对比,分析多尺度结构和动力学的成因,进而拓展计算,从而得到许多无法直接从中子散射实验中得到的细节,解释经验参数背后的物理.迄今为止,该研究方法已经取得众多成果,如英国散裂中子源(ISIS)的Soper等1116使用无序小分子中子全散射结合蒙特卡罗模拟的方法分析了多种小分子溶液的结构;美国麻省理工的Chen等17使用准弹性中子散射与计算机模拟结合的方式研究分析了水在220 K附近时的动力学转变过程;西班牙材料物理中心的Colmenero等18使用小角中子散射、准弹性中子散射结合分子动力学模拟的方式系统研究了线性共聚高分子中的特征结构、动力学行为与单体间关系;加拿大不列颠哥伦比亚大学的Cullis等19利用小角中子散射与计算机模拟结合的方式得到了包裹siRNA的脂质体结构;瑞典隆德大学的Wacklin等20利用分子动力学模拟分析中子反射的实验结果表征了生物脂质体膜的特征结构;埃克森美孚石油的Lpez-Barrn等21使用小角中子散射、广角X射线衍射及分子动力学模拟结合的方式分析了聚乙烯基联苯(poly(vinylbiphenyl),PVBP)中的局部特征结构并得到了实空间图像;上海交通大学的洪亮等22,23使用准弹性中子散射与计算机模拟结合的方式分析了低温下水的行为与蛋白质行为的关系:辅仁大学的Yang等24使用小角中子散射结合计算机模拟的方式系统研究了多种蛋白质在溶液中的结构等.本课题组结合无序大分子中子全散射、准弹性中子散射与计算机模拟,开发了一套解析无序大分子最可几全原子结构和动力学的方法.从Soper教授20世纪80年代开创的实验势结构精修(empirical potentialstructure refinement,EPSR)模拟无序液体全散射数据的成功经验出发16,引入GPU加速和并行计算25,从多组分小分子混合体系开始,随后在大分子的稀溶液、浓溶液、溶胀以及熔体和玻璃体系中进行了实验验证,得到了大分子的最可几全原子结构与动力学图像,分析了大分子从最可几全原子位置到整体构象的多尺度特征结构,以及从皮秒到微秒的多模式(mode)运动行为2632.该方法可以有效解释不同单体、不同原子间相互作用对多尺度空间结构和动力学的影响.本文先介绍了所用研究方法的基本原理,然后介绍了使用该方法在小分子溶液、大分子溶液、大分子熔融态和玻璃态体系中对结构与动力学行为研究得到的成果.2理论基础在中子全散射中,微分散射截面F(q),即qdd()代表单个散射体将中子散射到立体角d内的概率33,其与散射体内粒子的空间结构有如下关系34:()()c bc c b brgrqrqrrqdd()=+24()1sind(1)_incoh20+2其中q为散射矢量,、代表不同的原子类型,c代表该原子数目占体系原子总数的比值,b代表散射长度,为体系中全部原子的数密度,g(r)为径向分布函数,代表以粒子为中心时,粒子数密度对平均密度的比值沿距离r的分布情况35,计算方式为:中国科学:化学2023 年第 53 卷第 4 期679grrrrr()=1()4(2)2其中代表原子的数密度.当使用氘代技术改变体系中氢(H)与氘(D)的比例后,同一结构可得到多条不同的散射曲线,其对应的实空间结构完全相同25,即:FFwwwwHHqq().()=.(3)MNMMNN111111其中F(q)为不同氢氘比下体系的散射结果,w为粒子间相互作用的散射权重因子:()wc c b b=2(4)ijijijijH为粒子间相互作用:()Hrgrqrqrrq()=4()1sind(5)ijij0+2i与j由w矩阵与H矩阵中对应元素的列号与行号决定.根据式(1),如果各向同性的无序大分子体系中有N种不同的原子,就存在N(N+1)/2即(CN+N2)种不同的g(r),即如果想解出所有的g(r),理论上需要有N(N+1)/2个不同衬度的独立实验.在实际研究中,我们可以通过调整样品中不同组分的氘代情况,得到尽可能多的独立散射结果,同时利用计算机模拟预测补充实验数据的不足:使用同一体系不同衬度下的多个I(q)数据与对应组分的计算机模拟结果进行对比,得到体系中所有原子的最可几位置,这可以使模拟结果的散射数据验证过程更加可靠.对于小分子体系的模拟,我们主要在Soper教授的EPSR方法的基础上发展了EPSR+的程序,添加了GPU加速与并行计算来进行拟蒙特卡罗(逆蒙卡)模拟25.逆蒙卡模拟在传统蒙特卡罗(蒙卡)模拟的基础上,以散射结果为目标,通过迭代优化模拟体系中使用的势能函数,从而优化计算机模拟最终得到的结构.该模拟方法与小分子体系的中子全散射实验方法的匹配度非常好,然而在大分子体系中,由于大分子特征尺寸较大,其结构信息是通过小角散射区域内的数据体现,往往没有明显的特征峰,因此小角散射的实验数据较难通过傅里叶变换得到准确的势能函数结果,这也导致迭代过程比较不稳定.在大分子体系中使用较多的是分子动力学(mole-cular dynamics,MD)模拟,即通过设置粒子间的相互作用关系,计算其受力情况及相应的运动状态来进行体系的模拟.粒子间的相互作用参数,即力场的设置,对于模拟结果具有决定性的作用.一般来说,在分子动力学模拟中会涉及的粒子间相互作用包括成键、键角、二面角以及非键相互作用36.较常用的成键相互作用能与键角相互作用能的计算方式均为简谐式,即:EK rr=()(6)rbondeq2EK=()(7)angleeq2其中脚标eq表示平衡位置,K为作用参数,二面角相互作用的计算略有不同,为EVVV=21+cos()+21cos(2)+21+cos(3)(8)torsion123其中Vi为对应项的傅里叶系数,为二面角.粒子间的非键相互作用一般包含两项,即范德华作用力与静电力,其势能函数为:Efqq errr=+4(9)n bijijijijijijijij.2121266其中f为权重因子,q代表粒子的电荷数,e为电子电荷,r为粒子间距离,与为范德华相互作用参数.通过式(6)(9)可以计算体系中所有粒子间的相互作用,进一步可以得到粒子的受力情况:F rU rr()=d()d(10)再利用牛顿力学原理计算粒子的运动状态.粒子的初始速度可以自行设置合理数值,也可以通过模拟体系所处的温度根据泊松分布进行分配.我们可以根据已有的可靠力场进行模拟,得到体系的最可几全原子结构,然后根据式(1)进行计算,并与中子散射实验结果进行比较:如果实验结果与同一体系多种不同衬度下样品的实验结果都吻合得很好,则可以认为得到了体系的最可几全原子结构.粗粒化模拟的原则和全原子模拟基本相似,需要预先设定体系的力场参数再进行模拟,力场所包含的势能相互作用基本上与式(6)(9)相同.为了能保证粗粒化模拟的结构与中子散射结果相吻合,可以利用迭代玻尔兹曼反演方法,通过全原子模拟的结果得到粗韩泽华等:原子间相互作用对无序大分子结构和动力学影响的中子散射研究680粒化体系中的力场参数,从而保证粗粒化体系在单体尺度上的结构与全原子的结构完全相同37.该算法的基本思路为通过粗粒化结果与全原子结果间的差异迭代势能函数,最终得到更精确的势能函数:U rkTgr()=ln()(11)0targetUrU rkTgrgr()=()+ln()()(12)ii+1calculatedtarget其中U为势能函数,脚标i与i+1代表迭代中的轮次,g(r)在这里代表的是两体相互作用所涉及的分布函数,既可以是键长、键角分布,也可以是非键相互作用粒子间的相对位置分布,target为全原子模拟的结果,calcu-lated为粗粒化模拟计算得到的结果.当然,MD模拟也有自己的局限性,其中最重要的一个问题是模拟结果严重依赖于力场参数的设置,而力场参数的调试又十分复杂,甚至同一物质针对不同特性进行观察时,不同力场间的参数设置都会略有差异.当然,对于力场比较成熟的组分来说,MD模拟通常能够得到可靠的结果.3研究实例3.1小分子溶液纤维素是自然界中分布最广、产量最高的天然高分子,其绿色溶解一直是学术界和产业界关心的核心问题3840.2000年前后,武汉大学的张俐娜院士团队发现碱-尿素的水溶液在低温下可以迅速溶解纤维素4144,但该现象的机理一直未被澄清4548.我们选择与纤维素具有相似单体结构的海藻糖作为模型分子,研究了其在碱-尿素水溶液体系中的溶解机理.海藻糖结构如图1a所示,其与纤维素的单体结构类似,且是唯一一种不会被碱氧化的二糖28,49.我们配制了三类样品进行中子散射实验,分别是海藻糖的氢氧化钠-尿素水溶液、海藻糖的氢氧化钠水溶液及氢氧化钠的尿素水溶液,每一类样品中又分别配制了三种中子衬度(图1b):全氘(D2O,文中代表氢氧化钠、尿素、水全部氘代)、半氘(HDO,文中代表一半氢氧化钠、尿素、水氘代)与全氢(H2O,文中代表氢氧化钠、尿素、水全部氢化)进行实验,共得到9组实验数据,结合EPSR+模拟,把每一帧空间结构的傅里叶变换与中子散射结果进行对比.如图1b所示(其中Tr代表海藻糖,Na代表氢氧化钠,Ur代表尿素),二者吻合得很好,可以认为我们得到了体系的最可几全原子实空间结构.通过数据统计结合实空间图像,得到了在一个海藻糖分子周围溶剂分子的分布情况,如图1c所示,并对海藻糖的碱-尿素水溶液中原子间相互作用进行了分析:碱金属阳离子主要在海藻糖的氧原子周围富集,形成第一水合壳层;尿素分子较大,无法进入海藻糖分子内部,它们通过与氢氧根形成带负电的复合物(com-plexation),经由碱金属阳离子的桥接,与海藻糖间接相互作用,形成第二水合壳层.考虑到纤维素单体与海藻糖的结构类似,因此,我们推测纤维素在碱-尿素水溶液中的低温溶解机理:体积较小的碱金属阳离子通过静电相互作用进入纤维素分子内部空间位阻较高的氧原子附近区域,打破其晶体结构50,并通过与尿素-氢氧根复合物的桥接在外侧形成维持区域稳定的复合壳层结构,该结构的形成使得纤维素自身间范德华相互作用减弱,从而避免了纤维素的聚集,使其溶解.根据这个实验,我们可以预测不同类型的碱溶液对纤维素的溶解能力与碱金属阳离子的尺寸直接相关:碱金属阳离子的体积越小(LiNaNaOHKOH)30.3.2大分子溶液在大分子溶液体系中,通过联用无序大分子中子全散射和计算机模拟重建体系的最可几三维全原子结构相对简单.无序大分子中子全散射的观察尺度横跨小角中子散射和中子衍射的观察范围.由于小角中子散射的分析方法是基于粗粒化模型,而中子衍射基于全原子模型,二者恰好互补,因此,无序大分子中子全散射的方式十分适合用来建立大分子溶液体系的最可几全原子模型.同时,由于小分子的氘代试剂容易购买,通过氘代实现衬度变换实验的难度也较低.本节分别介绍在大分子稀溶液和溶胀体系中重建三维全原子结构和动力学的实验.3.2.1聚(N,N-二乙基丙烯酰胺)在乙醇水溶液中的共混不互溶机理研究共混不互溶是指大分子在两种良溶剂的混合溶剂中溶解度降低的现象,该现象明显无法使用溶度参数中国科学:化学2023 年第 53 卷第 4 期681和相似相溶的经验原则解释5155.我们以聚(N,N-二乙基丙烯酰胺)(PDEA)的乙醇水溶液为模型,研究了该现象的机理.如图2a的相图所示,当乙醇的摩尔浓度小于30%时,PDEA的乙醇水溶液是LCST体系;当乙醇的摩尔浓度大于30%时,混合溶剂是PDEA的良溶剂.因此,我们选择乙醇摩尔浓度为4.1%、8.7%、14.2%及60.5%时,于均相区、靠近相边界的20在ISIS的NIMROD谱仪上进行了中子全散射的实验56,每种浓度均配制了三种中子衬度(分别为全氘、半氘与全氢)的样品,同时使用分子动力学模拟的方式得到了体系实空间的全原子结构.模拟结果通过傅里叶变换得到的散射曲线与溶剂全氘样品的实验结果对比如图2b所示,二者完全重合,其他衬度下二者的结果也完全吻合(详见文献26的支撑材料),因此可以认为模拟结果重建了该体系在不同乙醇浓度下的最可几全原子结构.通过计算机模拟,可以得到不同乙醇浓度下PDEA的分子构象及其周围溶剂分子的分布情况(图2c),用来进行共混不互溶机理的分析.共混不互溶体系中,PDEA与两种溶剂分子间主图 1(a)海藻糖结构示意,其中不同位置的氧原子O使用标号1、2、3、4代表,C表示碳原子,M与H分别代表和碳相连与和氧相连的氢原子;(b)中子全散射实验数据(黑色虚线)、ESPR模拟结果(红线)的对比及其差异(蓝线);(c)EPSR+模拟结果中海藻糖分子周围Na离子(蓝色)与尿素上氧分子(红色)的空间分布情况28(网络版彩图)Figure 1(a)Structure of trehalose molecule.The oxygen atoms are labeled as O and tagged as 1,2,3,4 for different positions.The carbon is labeledas C.M and H are for hydrogen atoms connected to carbon and oxygen separately;(b)experimental neutron total scattering profiles(dashed blacklines),the EPSR fitted neutron structure factors(red lines)and their differences(blue lines);(c)distribution of Na+(blue)and the oxygen atoms of urea(red)around trehalose 28(color online).韩泽华等:原子间相互作用对无序大分子结构和动力学影响的中子散射研究682要的相互作用为氢键和范德华相互作用.依据图2c,我们可以分别考察体系中这两种相互作用随乙醇浓度的变化.体系中水和乙醇的逾剩焓(excess enthalpy)随着乙醇浓度的增加先减小后增大,在浓度接近30%左右时负值达到极小值,即该浓度之下水和乙醇的相互作用是随乙醇浓度的增加而增强的;随着乙醇浓度的增加,PDEA-水之间氢键减少的速率快于PDEA-乙醇间氢键的增加速率,PDEA与溶剂分子形成的氢键数量整体是减少的;PDEA与溶剂分子间范德华相互作用随着乙醇浓度的增加,有明显的先升后降趋势.结合以上现象,我们得到结论:PDEA在乙醇水溶液中的共混不互溶现象发生的必要条件是混合溶剂间的相互作用强于二者与溶质间的相互作用.图2d解释了PDEA在乙醇水溶液中共混不互溶的分子图像:乙醇浓度较低时,随着乙醇浓度的增加,PDEA与溶剂间的氢键相互作用快速变弱,而乙醇和水优先复合,浓度较低的乙醇会渗透到水四面体(tetrahedron)结构间,导致体系中的范德华相互作用没有显著增强,吉布斯自由能小于零,PDEA塌缩,发生了共混不互溶;随着乙醇浓度进一步增高,乙醇形成之字形(zigzag)结构,它们和PDEA间范德华相互作用显著变强,因此PDEA又可以溶解.3.2.2聚乙二醇溶胀体系中水的动力学研究由于氢的非弹性中子散射截面远大于氘,因此准弹性中子散射可以用来直接测量体系中氢原子的动力学过程57,即通过准弹性中子散射观察氢化高分子的氘代溶剂与氘代高分子的氢化溶剂时,得到的均是氢化部分的动力学过程.该原理可与核磁共振氢谱测量图 2(a)PDEA的乙醇水溶液体系相图,中子全散射实验条件为20下乙醇浓度分别为(1)4.1 mol%,(2)8.7 mol%,(3)14.2mol%,(4)60.5 mol%;(b)全氘溶剂样品中子全散射实验数据与模拟结果对比;(c)PDEA构象及第一水合壳层中小分子分布情况,其中红色为水分子,蓝色为乙醇分子;(d)PEDA在乙醇水溶液中共混不互溶的原子图像26(网络版彩图)Figure 2(a)Phase diagram of PDEA in water-ethanol solution.The neutron total scattering experiments were conducted at 20 C at different ethanolconcentration,e.g.,(1)4.1 mol%,(2)8.7 mol%,(3)14.2 mol%and(4)60.5 mol%;(b)neutron total scattering profiles and MD simulation results ofPDEA in fully deuterated water-ethanol solutions with different ethanol concentrations.(c)Snapshot of PDEA and its first solvation layer of water andethanol(the red dots are water molecules and the blue dots are ethanol molecules);(d)the diagram of cononsolvency of PDEA in water-ethanolsolutions 26(color online).中国科学:化学2023 年第 53 卷第 4 期683进行类比.使用准弹性中子散射结合分子动力学模拟的方式,我们研究了大分子水合壳层中结合水的动力学行为.水合壳层中水的结构与动力学行为对生物大分子正常发挥生理功能起到了重要的作用,研究显示,生物大分子从低温升至200 K时才会开始发挥生物功能,这恰好也是水合壳层发生动力学转变的温度22.有些研究认为生物大分子在200 K左右的活性变化是由水合壳层内水的运动决定的(这也是保存生物样品的冰箱温度设为200 K的原因)5861.为了分析生物大分子的转变与水合壳层的关系,首先需要研究200 K附近水合壳层的动力学行为6264.聚乙二醇(polyethylene glycol,PEG)是结构最简单的水溶性高分子,便于合成分子量和分子量分布相近的氢化、氘代聚合物,且PEG中的氧-氧间的相对位置(即氧-氧间g(r)分布与纯水相似,因此对水合壳层结构的影响最小,适合作为大分子模型来进行相关研究65.本研究中,我们使用质量浓度为85%的PEG与15%的水组成的体系作为模型.此时水含量极少,由于直接滴加水到PEG会造成水在PEG内部的不均匀分布,我们不能通过滴加水的方式制备样品).实验中,我们将完全烘干的PEG与一瓶水置于密闭容器,通过调节体系的温度来调节水的蒸气压,使水蒸气渗透进PEG,进行溶胀.这种制备方式下得到的样品仍呈固态,不具备任何流动性.首先,我们进行小角中子散射实验证明了所使用的氢化和氘代PEG具有相同的分子结构.图3是进行衬度匹配的不同浓度下PEG水溶液的小角散射结果,实验于绵阳反应堆中子源的小角谱仪“狻猊”上进行66.三种样品均使用摩尔比为1:1的h-PEG/d-PEG与0.66:0.34的D2O/H2O来进行溶剂与溶质分子间的衬度匹配(图3a).通过对图3b实验数据的分析可知单根PEG分子链的构象67.随着PEG浓度的提高,拟合结果中排斥体积参数由3/5(1 wt%)变为1/2(30 wt%,60 wt%),表明PEG在良溶剂中的排斥体积逐渐消失.小角中子散射实验的结果验证了相同分子量和分子量分布的氘代与氢化PEG具有相同的构象结构68.之后我们使用d-PEG的H2O溶胀以及h-PEG的D2O溶胀体系(水含量15%)进行了准弹性中子散射(quasi-elastic neutron scattering,QENS)实验.准弹性中子散射实验分别在美国国家标准与技术研究院中子散射中心(NCNR)的HFBS谱仪与ISIS的IRIS谱仪进行,二者的能量分辨率分别为1和17.5 eV69,70,在两台谱仪上均进行了从10 K到室温的弹性扫描以及在260、278、298、308、318 K下的准弹散射实验.如图4a,b所示,在不同的能量分辨率下,水的弹性扫描结果中转变温度并不相同:在HFBS上观察到水的转变温度为209 K,而在IRIS上观察到水的转变温度为220 K.如果结合水在升温过程中的转变是相变,则不同仪器的分辨率不应影响观测到的转变点位置;而实验中不同分辨率的仪器观测到的转变点位置不同,且能量分辨率越高的谱仪观测到的转变点温度越低,这说明在220 K附近水的转变是动力学转变:升温过程中,结合水先开始运动,结合水的运动“润滑”了PEG主链的扭转71.图 3(a)1.4 wt%的h-PEG/d-PEG水溶液随H/D比改变的小角散射强度,衬度匹配点是D2O:H2O=0.66:0.34;(b)1 wt%、30 wt%以及60 wt%的h-PEG/d-PEG水溶液在衬度匹配时的小角散射结果31(网络版彩图)Figure 3(a)SANS data of 1.4 wt%h-PEG/d-PEG in water as a function of H/D ratio.The contrast match point is D2O:H2O=0.66:0.34;(b)SANSdata of 1 wt%,30 wt%and 60 wt%h-PEG/d-PEG water solution at contrast match point 31(color online).韩泽华等:原子间相互作用对无序大分子结构和动力学影响的中子散射研究684之后,我们通过计算机模拟重建了体系的三维结构,并得到了220 K时体系的动力学过程.通过模拟结果可以计算220 K时水合壳层水与PEG间形成氢键及氢键网络的寿命,将氢键网络寿命与准弹性中子散射实验得到的特征弛豫时间进行对比,二者吻合.在此基础上,我们从过冷态结合水的全原子结构中提取了结合水的运动轨迹,分析了其动力学行为.图4c描绘了220 K时具有代表性的一个水分子的运动轨迹,轨迹为局部区域运动与跨区域间迁移的结合.因此,我们判断,PEG的结合水会围绕其主链在空间中形成许多小的“水池”(pool),结合水分子的运动包括两部分,在水池内的运动,以及脱离束缚在水池间的长距离迁移,二者结合体现了低温下结合水动力学行为的特征31,72.接着,我们与上海交通大学洪亮教授课题组32合作,总结出在200 K左右时,过冷水在水合壳层中的转变是一个普遍存在的现象,其广泛存在于多种生物大分子、有机大分子及二维材料中.3.3大分子熔融与玻璃化体系在大分子熔融和玻璃化体系中联用中子散射与计算机模拟重建最可几全原子结构的最大障碍在于计算机计算能力的限制,因此需要在全原子模拟的基础上引入粗粒化模拟及相应的粗粒原子回填方式73.本节中,我们以聚苯乙烯(polystyrene,PS)为例进行详细说明.3.3.1聚苯乙烯熔体与玻璃态的最可几全原子结构诺贝尔奖得主Anderson教授曾指出:“在固态理论中最深刻也是最有趣的未解决问题可能就是玻璃的性质与玻璃化转变,这可能成为未来研究中的下一个重大突破74.”对大分子在玻璃化形成与转变过程中的结构变化及动力学行为的研究历来是软物质领域最为重要的课题之一75.诸多相关的理论如自由体积(freevolume)76、熵(entropy)77,78、模式耦合(mode cou-pling theory,MCT)79和阻塞(jamming)理论80等充满了争议,其主要原因在于无法直接观测到高分子玻璃化的原子图像.我们使用聚合度和聚合度分布相同、氘代方式不同的PS,包括全氢化的PS-h8、氘代苯环的PS-d5以及全氘代的PS-d8配制了多种样品81,在英国ISIS的NIMROD谱仪与日本J-parc的NOVA谱仪上进行了玻璃化温度之上(453 K)降温到玻璃化温度之下(328 K)的中子全散射实验56,82,并进行了相应的分子动力学模拟(包括全原子与粗粒化),二者结合分析了PS中的特征结构及其在玻璃化转变过程中发生变化的原子尺度原因27.以PS-d8为例,其中子散射结果与计算机模拟结果的对比如图5a,b所示,二者吻合,其他组分样品的结果类似(详见文献27),因此可以认为通过计算机模拟重建了PS在不同温度下的最可几全原子结构.为了能够通过原子类型定位具体位置,我们将原子按主链碳(MCC)、主链氘(MCD)、苯环碳(PRC)和苯环氘(PRC)的方式进行分类,用以分析不同原子对间相互作用对体系最终散射结果的贡献,这有助于分析特征峰对应的原子尺度的具体结构.如图5a所示,450 K时,PS在0.6、1.4与1.8 1处有三个特征峰,对应实空间尺寸分别为10、5与3.4;且0.6 1的特征峰随温度升高而明显反常增强,1.4 1图 4PEG溶胀体系弹性扫描结果归一到10 K数据时随温度的变化趋势.(a)HFBS;(b)IRIS;(c)220 K时一个水分子的运动轨迹,总观察时长500 ns,步长20 ps,盒子单位为nm(网络版彩图)Figure 4The temperature-dependent summed elastic neutron scattering scan of PEG swelling water solution,normalized to 10 K from HFBS(a),and IRIS(b).(c)Trajectory of one water molecule at 220 K for 500 ns with the time step of 20 ps.The unit of coordinate is nm(color online).中国科学:化学2023 年第 53 卷第 4 期685的特征峰在玻璃化转变附近有较弱反常增强,1.8 1的特征峰不随温度变化.我们通过重建的PS最可几全原子结构,对三个特征峰代表的特征结构及其随温度变化的原因进行了研究与分析.根据式(1),在PS-d8中总的中子散射谱可以分解为10类(N(N+1)/2)相互作用,分别为:主链碳-主链碳、主链碳-主链氘、主链碳-苯环碳、主链碳-苯环氘、苯环碳-苯环碳、苯环碳-主链氘、苯环碳-苯环氘、主链氘-主链氘、主链氘-苯环氘以及苯环氘-苯环氘.我们具体分析了这10组相互作用随温度变化的情况,发现只有3组相互作用,即主链碳-苯环碳、苯环碳-主链氘以及苯环碳-苯环氘的强度随温度有明显变化(图5ce).在此基础上,我们分析了这3组相互作用对特征峰变化的影响以及其对应的实空间结构变化.0.6 1的特征峰主要对应PS链段间的平均距离(图6b)83.该峰强度变化受主链碳-苯环碳及苯环碳-主链氘间相互作用的影响(图5c,d),其变化反映了玻璃化过程中动力学行为的变迁,具体为苯环的绕主链运动.随着温度的降低,苯环与主链间的相互作用减弱,苯环的绕轴运动减缓,链段间平均位置靠近,因此峰的强度(负值)减小且峰位向高q移动;在发生玻璃化转变后,苯环的绕轴运动被逐步冻结,因此峰的强度在玻璃化两侧发生跃变;之后由于苯环的运动已被冻结,因此峰强与峰位随温度的变化不再明显.1.4 1的特征峰对应链内苯环与链间苯环多种相互作用(图6a,b)84,该峰强度变化受苯环碳-苯环氘间相互作用影响(图5e),反映了玻璃化过程中苯环翻转(flipping)动力学行为的变迁.随着温度的降低,苯环翻转运动的减弱导致苯环间的相对位置偏离平衡位置,反映了苯环间相互作用变弱85,因此峰的强度(负值)减小且峰位向高q移动;在发生玻璃化转变后,苯环的翻转运动同样被逐步冻结,因此峰的强度在玻璃化两侧发生跃变;之后由于苯环的运动已被冻结,峰强与峰位随温度的变化不再明显.我们有理由相信,0.6和1.4 1特征峰发生跃变的具体温度并不相同,因为它们分别对应不同尺度运动模式的冻结.1.8 1的特征峰对应的是苯环内碳和间位氢(氘)间的固有结构(图6c),该结构不随温度变化.通过中子全散射结合计算机模拟的方式,我们得到了PS的特征峰对应的实空间结构,并分析了其在玻璃化转变中发生的改变,构建了原子图像.在此基础上,我们进一步使用分子动力学模拟的结果分析了PS玻璃化过程中的动力学行为.图 5(a)PS-d8的中子全散射实验结果;(b)PS-d8的计算机模拟结果;PS-d8体系中不同原子对间相互作用对总散射曲线的贡献:(c)主链碳-苯环碳,(d)苯环碳-主链氘,(e)苯环碳-苯环氘27(网络版彩图)Figure 5(a)Neutron total scattering profiles and(b)MD simulations of PS-d8.Different pair interaction contribution to the scattering intensity ofPS-d8:(c)MCC-PRC,(d)PRC-MCD,(e)PRC-PRD 27(color online).韩泽华等:原子间相互作用对无序大分子结构和动力学影响的中子散射研究6863.3.2聚苯乙烯玻璃化过程的动力学自由体积和熵理论对玻璃化大分子的运动方式有不同的描述:自由体积理论认为体系中必须有能容纳其单体运动的空间,大分子链才能运动8688;而熵理论认为不需要有足够的自由体积,大分子们也可以协同运动(就像人们在拥挤的舞池中翩翩起舞)8991.我们的实验可以有效地分析PS在玻璃态的运动方式:计算机模拟已经得到了PS全原子结构的一帧帧图片,所有帧傅里叶变换的平均结果与中子散射的数据相符,即所有帧的图像都反映PS全原子在空间的系综平均位置,而帧与帧间的联系则能反映PS在玻璃化过程中的运动方式.如果自由体积理论正确,我们一定可以在这些图像中找到可以容纳PS单体的自由体积.然而我们无法在任一PS的实空间结构中找到哪怕一个可以容纳PS单体的“空穴”29,这说明自由体积描述的只是一个系综平均的统计参数,并非对应具体的实空间“空穴”体积.因此我们转而研究PS在玻璃化形成过程中的协同运动行为.由于PS模拟体系中原子数目庞大,全原子模拟很难提供动力学行为分析所需时间尺度(100 ns1 s)的数据,因此我们使用PS单体作为粗粒化中的基本粒子,以PS在453 K时的全原子构象为基础,通过玻尔兹曼反演方法得到粗粒化模型下此时体系的力场参数.粗粒化模拟中只能使用约化的温度、压强等参数进行研究,因此,我们设定真实温度为453 K时对应粗粒化体系的约化温度为T=1.0,在此基础上,我们进行了体系降温过程中不同温度的模拟37,并通过该降温过程研究了标度关系下玻璃化形成过程中体系的动力学行为.首先我们通过计算体系的均方位移(mean squaredisplacement,MSD)与体系的自相关函数(incoherentcorrelation function)得到体系开始发生玻璃化转变的温度是T=0.459294.我们还将其与玻璃化进一步发生时T=0.42的模拟结果进行了对比.通过计算体系的四点相关函数(four point susceptibility)4可以考察体系动力学异质性的特点,其计算方式如下29:Q a tN(,)=1e(13)iNra2=12i22()a tNQ a tQ a t(,)=(,)(,)(14)42222其中Q2为单体粒子的时间二点相关函数,N为单体粒子总数,a为预先设定的值,r为时间尺度t上粒子的位移.PS在温度T=0.45与T=0.42时四点相关函数的结果如图7a,b所示:T=0.45时,极值出现在t=0.99 s处,约为7.5;T=0.42时,极值出现在t=4.6 s处,约为50.8.4极值的位置代表体系动力学异质性最明显的时刻,此基础上,我们可以在实空间中寻找协同运动区域90,91,具有代表性的协同运动区域如图8a,b所示:在T=0.45时,协同运动区