外力驱动作用下高分子链在表面吸附性质的计算机模拟∗

李洪 艾倩雯 汪鹏君† 高和蓓 崔毅 罗孟波

1)(温州大学数理与电子信息工程学院,温州 325035)

2)(温州职业技术学院信息系,温州 325035)

3)(浙江大学物理系,杭州 310027)

1 引 言

近年来,随着高分子在化学工业和生物领域的应用[1,2]日益增加,研究控制高分子性质的物理机制对高分子的应用具有重要的意义[3].从自然界中的天然橡胶、蛋白质到工业应用的胶体、色谱层析法等都与高分子的吸附现象相关,高分子链在吸引表面上吸附特性的研究受到人们的广泛关−注[4−12].以往对高分子链在具有吸引作用界面上的研究,主要集中在高分子链在表面上的临界吸附温度和构象在吸附过程中的变化[9−13].高分子链与表面之间相互作用强度是影响高分子吸附的主要因素,如果表面存在强相互吸引作用,那么高分子链在表面的吸附量增加,高分子链趋于在表面附近形成较薄的吸附层;而当高分子链与表面相互吸引作用较弱时,表面上的高分子链伸展,在溶液中形成分子刷[14].在高分子溶液中增加外力驱动作用将影响高分子的吸附及吸附态的构象[15].高分子链的尾端或中间某一个单体受到外力作用,其构象将受外力大小和方向的影响[16−18].因此,高分子链在具有吸引作用的表面附近受外力驱动作用下的吸附性质研究具有重要的意义和应用前景.

伴随高分子物理实验研究的发展,高分子的理论计算与模拟也发展成为高分子学科的一种重要研究手段[19−21].粗粒化分子动力学可以模拟研究高分子链的构象性质[22]和生物大分子的能量及结构[23].外力对高分子链构象的影响也可以采用模拟方法研究[16−18].高分子链由于分子热运动而具有各种不断变化的构象,并且在某个时间的构象是完全随机的,构象数量非常大,然而蒙特卡罗方法可以直接模拟高分子链的随机性问题.单个高分子链的构象统计也是一个复杂的计算问题,Metropolis的重要性抽样方法可以有效地模拟计算高分子链的性质[24−27].本文采用键长涨落模型和协同运动算法模拟在吸引表面附近的高分子链受到恒定的外力场驱动作用,研究了外力驱动作用对高分子链吸附性质的影响以及通过构象形变来交叉校验临界吸附点.

2 模拟模型

高分子链模型采用三维的简立方格模型,其中链长为N的高分子链由N个单体组成,高分子链中相邻的单体通过键长可涨落的键相连,其键长取值为高分子链的两端均为自由即非接枝的高分子链,高分子链被放置在两个平行表面之间,表面间距D>Nv,其中v=0.588为Flory指数.因此,在Z=0和Z=D处分别放置一个均质且不可穿透的表面,在Z=0处的表面对高分子链单体具有吸引作用,而在Z=D处的表面对高分子链不存在吸引作用,目的是使高分子链不会远离具有吸引作用的表面,在模拟中取D=100.高分子链在X和Y方向将满足周期性边界条件.在两个表面间施加一个平行于X轴正方向的均匀外力场,在整个模拟过程中高分子链上的每个单体一直受到该恒定的外力场驱动作用,当构象发生变化时,力的方向始终保持不变.单体与单体之间的相互关系如下:1)所有单体均满足自回避的条件,即两个单体不能占据同一个格点;2)键与键之间不允许交叉;3)非键相邻单体之间仅考虑体积排斥作用.

模拟过程如下: 首先采用Rosenbluth-Rosenbluth链生长方式[28]随机地生成一条链长为N的高分子链,在高分子链的生成过程中,每个单体有26个可选择的矢量方向;然后,让高分子链做随机的布朗运动,在整条高分子链中随机选择一个单体进行运动,单体在运动时有6个矢量方向{(1,0,0),(−1,0,0),(0,1,0),(0,−1,0),(0,0,1),(0,0,−1)},并通过键长涨落算法[26,27]和协同运动算法[29−31]产生新的构象.当选定的单体在尝试运动后有三种可能的情况,如图1所示:1)单体k尝试运动后,前后两端的键都超出键长范围,如图1(a)和图1(b)所示,尝试运动失败,保留原始位置;2)单体i尝试运动后,两端的键仍然满足键长范围,如图1(a)所示的单体i按箭头方向运动后如图1(c)所示,这种情况下采用键长涨落算法产生新的构象;3)单体i尝试运动后,一端的键仍然在允许的键长范围内,但另一端的键超出了键长范围,如图1(d)所示,单体i按箭头方向运动,此时按协同运动算法,超出键长范围这端的近邻单体依次向前一单体位置运动,直到遇到可以满足键长允许条件的单体j为止,从而产生新的构象.运用Metropolis重要性抽样方法来确定新构象的接受概率P,即P=min{1,exp(−∆E/KBT)}(KB为玻尔兹曼常数,∆E为新旧构象产生的能量差).其能量变化从吸附能和外力驱动能两个方面来考虑:1)单体与下表面之间存在相互吸引作用,其作用强度为ε=−1,高分子链每次运动前后接触能变化记为∆Es=ε∆M,其中∆M为运动前后高分子链在表面接触数的变化量;2)单体均受到沿x轴正方向的外力F驱动作用,高分子链每次运动受外力驱动产生的能量变化记为其中∆xk为协同运动团簇中的每个单体k在X方向上发生的位移;3)每次运动前后产生的能量变化为∆E=∆Es+∆EF.

图1 高分子链中单体尝试运动的示意图 (a)尝试运动前的状态;(b)单体k尝试往右运动,相邻两键均断开,尝试运动失败;(c)单体i尝试向下运动,满足键长涨落条件;(d)单体i尝试向左运动,然而单体i到j之间的所有单体进行协同运动Fig.1.Schematic diagram of the monomer in polymer attempt to move:(a)The state of polymer before movement;(b)the monomer k tries to move right;then two bonds between monomer k and its adjacent monomers are disconnected so that the trying movement is failure;(c)the monomer i tries to move downward,two bonds both meet the bond fluctuation conditions;(d)the monomer i tries to move left,however,all the monomers from i to j try to move cooperatively.

在模拟过程中,高分子链中的单体不断地进行布朗运动,我们把蒙特卡罗步(MCS)作为一个时间计量单位,在每一个蒙特卡罗步中高分子链的所有单体试图平均运动一次.在每个温度下,高分子链都将经历时间τ=2.5×N2.13MCS[9]来达到一个平衡状态,在后续的100τ MCS记录高分子构象样本.对于链长为N的高分子链,均产生1000个独立构象用来对结构求平均,以确保研究结果的准确性.

3 结果与讨论

3.1 无外力驱动作用时高分子链的吸附特性

高分子链随着温度降低会产生从脱附状态到吸附状态的转变,该相变点称为临界吸附点,在此温度称为临界吸附温度Tc.高分子链的吸附特性可以通过不同温度下平均表面接触数M来表示.图2描述了高分子链在不同温度的平均表面接触数M与链长N的关系.我们可以通过分析在不同温度下高分子链的M来估计无限长高分子链的临界吸附点,即采用有限尺寸标度方法来计算高分子链的临界吸附温度Tc,其标度关系可以表示为[10,32,33]

(1)式中t=(T−Tc)/Tc标度温度,ϕ为交叉指数和1/δ为另一个临界指数.(1)式清晰地表达了M的值在不同温度下有着不同的变化行为.当温度T从T>Tc到T

在标度理论中,临界吸附温度Tc和交叉指数ϕ可以描述高分子链的吸附特性[9].它们可以从平均表面接触数M与链长N的指数关系计算得到.为更精确地计算临界吸附温度Tc,在Tc附近其他温度的平均表面接触数M可以通过模拟数据的二次插值计算来获取.通过以上方法,我们计算得到了在外力F=0时高分子链的临界吸附温度Tc=1.95和指数ϕ=1.其临界吸附温度大于未采用协同运动算法所得到的Tc=1.65[9],两者差异的原因是高分子模型中采用不同的运动方式.

图2 在临界吸附温度附近M-N的双对数关系其中M为平均表面接触数;N为高分子链的链长,N=40—400;外力F=0;临界吸附温度Tc为1.95;交叉指数ϕ为1Fig.2.The double logarithmic plot of M-N near the critical adsorption temperature,whereMis the average number of surface contacts and N is the chain length.The polymer chain is changed from N=40 to 400.The critical adsorption temperature Tcis 1.95 and the cross-index ϕ is 1.

3.2 外力驱动作用下高分子链的吸附特性

图3 在不同温度下高分子链的平均表面接触数M与链长N的关系,其中外力F=0.3,链长N=40—400Fig.3.The relationship between the average number of surface contactsMand the chain length N at different temperatures T,where the external force is F=0.3,and the chain length is from N=40 to 400.

图4 平均表面接触数M与温度T的关系Fig.4.The relationship between the average number of surface contactsMand the temperature T.

图5 高分子链的温度T与外力F的伪相图 其中链长N=200,DS为脱附态和AS为吸附态;插图的三种构象分别为(a)F=0.1,T=1.0,(b)F=10,T=1.0和(c)F=10,T=0.2Fig.5.The pseudo–phase diagram of the polymer chain between the desorbed state(DS)and the adsorbed state(AS)for the temperature T and the external force F,in which the chain length is N=200.Three conformations of the insets are(a)F=0.1,T=1.0,(b)F=10,T=1.0 and(c)F=10,T=0.2.

3.3 外力驱动对高分子链构象的影响

其中,N为高分子链的链长,ri为高分子链中第i个单体的位置矢量,rcm为高分子的质心位置矢量,其计算公式为

为了研究吸引表面附近的高分子链构象受到外力驱动的影响,我们模拟了高分子链的均方回转半径及X,Y,Z 方向的分量与外力F 的关系,如图6所示,其中温度T=1.当外力F较小时,高分子链的及其分量随F的增大而保持不变,即此时的外力F不足以使高分子链构象发生变化.当外力F进一步增大时的X和Y分量出现了分叉,随着外力F的增加逐渐增大,垂直外力的Y方向变小,而且回转半径随外力F的增加而减小,即外力驱动作用使高分子链构象发生形变.直到F=Fc时,达到极小值,此时减小到极小值而增大到极大值,而且它们的值相等后随外力F的增加都几乎不再变化,此时在垂直外力方向的Y和Z分量相等,说明没有受到表面的限制即高分子链处于脱附状态.当外力增大到一定值后,高分子链的及其三个分量随F的增大而几乎不再变化.因此,高分子链的临界吸附点也可以从高分子链的极小值,或Y和Z分量的变化(即高分子链构象形变)来粗略估计.

为了从高分子链构象角度进一步说明高分子链在表面吸附受外力的影响,我们分别计算了均方回转半径在Y,Z方向上的分量和与温度T的关系,如图7所示.从图7可以观察到,高温时高分子链这两个分量的值相等,然后在T=Tc出现分叉,随温度T的减小而增大,随温度T的减小而减小.其原因是在临界吸附温度以下高分子链受到表面吸引作用而靠近表面,又因为体积排斥作用使高分子链沿表面发生伸展.分叉点的温度值随着外力F 的增加而减小,且与其临界吸附温度Tc相一致.

图6 高分子链均方回转半径以及X,Y,Z方向分量与外力F的关系Fig.6.The mean square radius of gyration and its components at the different external force F.

图7 均方回转半径在Y,Z方向上的分量和与温度T的关系,其中链长N=200,外力F=0,0.5,1.0和10Fig.7.The components of mean square radius of gyration in the Y,Z direction and at different temperature T,where the chain length is N=200 and the external forces are F=0,0.5,1.0 and 10.

高分子链在表面的吸附性质与温度相关,我们根据图5的伪相图,选择两个温度和来说明高分子链吸附和构象受外力F驱动的不同作用.高分子链的平均表面接触数M和均方回转半径在Z方向上的分与外力F的关系如图8所示.当T=1.2时,M随外力F的增大而减小,而且其M值逐渐趋于0,即高分子链发生了由吸附状态向脱附状态的转变;当T=0.2时,M随外力F的增大而减小,但M仍然接近N,即高分子链中虽存在部分单体脱附但整条链仍处于吸附状态.当F=0时,T=1.2和0.2时→ 0,高分子链处于接近表面的吸附状态.而随着F的逐渐增大,温度T=1.2时单调递增直到F>0.4后保持不变,此时明显大于0,其中Fc=0.4与图5中的相应值一致;而温度T=0.2时即使外力F>2后仍然接近于0,其值明显小于温度T=1.2时的,说明高分子链仍然在吸引表面附近.因此,在温度T=1.2时,当外力F足够大时,使高分子链发生从吸附状态到脱附状态的相变,而在温度T=0.2时随外力增加没有发生脱附相变.在两个不同温度T>和T<,高分子链吸附性质和构象性质受外力驱动作用的不同影响,其性质与图4的伪相图相一致.

图8 (a)高分子链平均表面接触数M和(b)均方回转半径在Z方向上的分量与外力F的关系Fig.8.(a)The average number of surface contactMand(b)the component of the mean square radius of gyration in the Z direction at different external force F.

4 结 论

本文采用键长涨落模型和协同运动算法模拟在吸引表面附近并受到平行表面的外力驱动作用的高分子链的热力学性质,研究了高分子链在不同大小的外力作用下的吸附特性.模拟发现高分子链的临界吸附温度Tc随外力F的增大而减小,据此构建了链长N=200的高分子链的脱附状态和吸附状态相对于温度T和外力F的伪相图.我们从高分子链的构象角度交叉校验了高分子链的临界吸附点,发现从高分子链的的极小值与从其Y和Z分量的变化估计的临界吸附点相一致.最后讨论了温度T>和T<两种情况的高分子链吸附性质和构象性质受外力驱动作用的不同现象,当