外力作用下高分子链在表面吸附性质的蒙特卡罗模拟方法与流程

本发明属于高分子的理论计算与模拟领域，具体涉及外力作用下高分子链在表面吸附性质的蒙特卡罗模拟方法。

背景技术

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

伴随高分子物理实验研究的发展，高分子的理论计算与模拟也发展成为了高分子学科的一种重要研究手段。粗粒化分子动力学可以模拟研究高分子链的构象性质和生物大分子的能量及结构。外力对高分子链构象影响也可以采用模拟方法研究。

技术实现要素：

本发明的目的是为了克服现有技术存在的缺点和不足，而提供一种外力作用下高分子链在表面吸附性质的蒙特卡罗模拟方法。

本发明所采取的技术方案如下：外力作用下高分子链在表面吸附性质的蒙特卡罗模拟方法，包括以下步骤：

步骤一、构建高分子链模型，采用三维的简立方格模型，采用rosenbluth-rosenbluth链生长方式随机地生成一条链长为n的高分子链，其中链长为n的高分子链由n个单体组成，在高分子链的生成过程中，每个单体有26个可选择的矢量方向，高分子链中相邻的单体通过键长可涨落的键相连，其键长取值为1、高分子链的两端均为自由，在z＝0和z＝d处分别放置一个均质且不可穿透的表面，高分子链被放置在两个平行表面之间，其中表面间距d>n^v，其中v为flory指数0.588，在z＝0处的表面对高分子链单体具有吸引作用，而在z＝d处的表面对高分子链不存在吸引作用，高分子链在x和y方向将满足周期性边界条件，在两个表面间施加一个平行于x轴正方向的均匀外力场，在整个模拟过程中高分子链上的每个单体一直受到该恒定的外力场驱动作用，当构象发生变化时，力的方向始终保持不变，单体与单体之间的相互关系如下：(a)所有单体均满足自回避的条件，即2个单体不能占据同一个格点；(b)键与键之间不允许交叉；(c)非键相邻单体之间仅考虑体积排斥作用；

步骤二、让高分子链做随机的布朗运动，在整条高分子链中随机选择一个单体进行运动，单体在运动时有6个矢量方向{(1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(0,0,1),(0,0,-1)}，当选定的单体在尝试运动后有三种可能的情况，(1)单体尝试运动后，前后两端的键都超出键长范围，尝试运动失败，保留原始位置；(2)单体尝试运动前后两端的键仍然满足键长范围，这种情况下采用键长涨落算法产生新的构象；(3)单体尝试运动后一端的键仍然在允许键长范围，但另一端的键超出键长范围，按协同运动算法，超出键长范围这端的近邻单体依次向前一单体位置运动，直到遇到可以满足键长允许条件的单体为止，从而产生新的构象；

步骤三、计算新构象的接受概率p，运用metropolis重要性抽样方法来确定，即p＝min{1,exp(-△e/kbt)},kb为玻尔兹曼常数，△e为新旧构象产生的能量差，△e从吸附能和外力驱动能两个方面来考虑：(1)单体与下表面之间存在相互吸引作用，其作用强度为ε＝-1，高分子链每次运动前后接触能变化记为△es＝ε△m，其中△m为运动前后高分子链在表面接触数的变化量；(2)单体均受到沿x轴正方向的外力f驱动作用，高分子链每次运动受外力驱动产生的能量变化记为其中协同运动团簇中的每个单体k在x方向上发生的位移△xk；(3)每次运动前后产生的能量变化为△e＝△es+△ef；

步骤四、以蒙特卡罗步mcs作为一个时间计量单位，在每一个蒙特卡罗步中高分子链的所有单体试图平均运动一次，在每个温度下，高分子链都将经历时间τ＝2.5×n^2.13mcs来达到一个平衡状态，在后续的100τmcs记录高分子构象样本。

对于链长为n的高分子链，随机产生1000个独立的初始构象，然后随机布朗运动，达到平衡态后对高分子形态和性质求统计平均。

本发明的有益效果如下：本发明为高分子链在具有吸引作用的表面附近受外力驱动作用下的吸附性质研究提供一种外力作用下高分子链在表面吸附性质的蒙特卡罗模拟方法，可以用来研究了外力驱动作用对高分子链吸附性质的影响以及通过构象形变来交叉校验临界吸附点,对高分子的应用研究具有重要意义。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，根据这些附图获得其他的附图仍属于本发明的范畴。

图1为高分子链中单体尝试运动的示意图，(a)尝试运动前的状态，(b)单体k尝试往右运动，相邻两键均断开，尝试运动失败，(c)单体i尝试向下运动，满足键长涨落条件和(d)单体i尝试向左运动，然而单体i到j之间的所有单体进行协同运动；

图2为在临界吸附温度附近<m>-n的双对数关系，其中<m>为平均表面接触数；n为高分子链的链长，从n＝40到n＝400；外力f＝0。临界吸附温度tc为1.95，交叉指数为1；

图3为高分子链的平均表面接触数<m>与链长n在不同温度t的关系，其中外力f＝0.3，链长从n＝40到n＝400；

图4为平均表面接触数<m>与温度t的关系，其中链长n＝200，外力f＝0,0.1,0.3,0.5,1.0和10；

图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.2；

图6为高分子链均方回转半径<rg²>以及x，y，z方向分量与外力f的关系，其中链长n＝200，温度t＝1；、

图7为均方回转半径在y、z方向上的分量<rg²>y和〈rg²>z与温度t的关系。其中链长n＝200；外力f＝0，0.5，1和10；

图8为(a)高分子链平均表面接触数<m>和(b)均方回转半径在z方向上的分量与外力f的关系，其中链长为n＝200，温度为t＝1.2和0.2；

图9为本发明的流程示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作进一步地详细描述。

如图9所示，一种外力作用下高分子链在表面吸附性质的蒙特卡罗模拟方法，包括以下步骤：

步骤一、构建高分子链模型，采用三维的简立方格模型，采用rosenbluth-rosenbluth链生长方式随机地生成一条链长为n的高分子链，其中链长为n的高分子链由n个单体组成，在高分子链的生成过程中，每个单体有26个可选择的矢量方向，高分子链中相邻的单体通过键长可涨落的键相连，其键长取值为1、高分子链的两端均为自由，在z＝0和z＝d处分别放置一个均质且不可穿透的表面，高分子链被放置在两个平行表面之间，其中表面间距d>n^v，其中v为flory指数0.588，在z＝0处的表面对高分子链单体具有吸引作用，而在z＝d处的表面对高分子链不存在吸引作用，高分子链在x和y方向将满足周期性边界条件，在两个表面间施加一个平行于x轴正方向的均匀外力场，在整个模拟过程中高分子链上的每个单体一直受到该恒定的外力场驱动作用，当构象发生变化时，力的方向始终保持不变，单体与单体之间的相互关系如下：(a)所有单体均满足自回避的条件，即2个单体不能占据同一个格点；(b)键与键之间不允许交叉；(c)非键相邻单体之间仅考虑体积排斥作用；

步骤二、让高分子链做随机的布朗运动，在整条高分子链中随机选择一个单体进行运动，单体在运动时有6个矢量方向{(1,0,0),(-1,0,0),(0,1,0),(0,-1,0),(0,0,1),(0,0,-1)}，当选定的单体在尝试运动后有三种可能的情况，(1)单体尝试运动后，前后两端的键都超出键长范围，尝试运动失败，保留原始位置；(2)单体尝试运动前后两端的键仍然满足键长范围，这种情况下采用键长涨落算法产生新的构象；(3)单体尝试运动后一端的键仍然在允许键长范围，但另一端的键超出键长范围，按协同运动算法，超出键长范围这端的近邻单体依次向前一单体位置运动，直到遇到可以满足键长允许条件的单体为止，从而产生新的构象；

步骤三、计算新构象的接受概率p，运用metropolis重要性抽样方法来确定，即p＝min{1,exp(-△e/kbt)},kb为玻尔兹曼常数，△e为新旧构象产生的能量差，△e从吸附能和外力驱动能两个方面来考虑：(1)单体与下表面之间存在相互吸引作用，其作用强度为ε＝-1，高分子链每次运动前后接触能变化记为△es＝ε△m，其中△m为运动前后高分子链在表面接触数的变化量；(2)单体均受到沿x轴正方向的外力f驱动作用，高分子链每次运动受外力驱动产生的能量变化记为其中协同运动团簇中的每个单体k在x方向上发生的位移△xk；(3)每次运动前后产生的能量变化为△e＝△es+△ef；

步骤四、以蒙特卡罗步mcs作为一个时间计量单位，在每一个蒙特卡罗步中高分子链的所有单体试图平均运动一次，在每个温度下，高分子链都将经历时间τ＝2.5×n^2.13mcs来达到一个平衡状态，在后续的100τmcs记录高分子构象样本。

对于链长为n的高分子链，随机产生1000个独立的初始构象，然后随机布朗运动，达到平衡态后对高分子形态和性质求统计平均，以确保研究结果的准确性。

下面通过上述方法研究高分子链在不同大小的外力作用下的吸附特性。

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

高分子链随着温度降低会产生从脱附状态到吸附状态的转变，该相变点称为临界吸附点，在此温度称为临界吸附温度tc。高分子链的吸附特性可以通过不同温度下平均表面接触数<m>来表示。图2描述了高分子链在不同温度的平均表面接触数<m>与链长n的关系。可以通过分析在不同温度下高分子链的<m>来估计无限长高分子链的临界吸附点，即采用有限尺寸标度方法来计算高分子链的临界吸附温度tc，其标度关系可以表示为：

在公式(1)中t＝(t-tc)/tc标度温度，为交叉指数和1/δ为另一个临界指数。在公式(1)中清晰地表达了<m>的值在不同温度下有着不同的变化行为。当温度t从t>tc到t<tc变化时，变量tn^1/δ的符号会发生改变。当t＝tc时，平均表面接触数<m>与链长n之间存在一个幂律关系：

在标度理论中临界吸附温度tc和交叉指数可以描述高分子链的吸附特性。它们可以从平均表面接触数<m>与链长n的指数关系计算得到。为更精确计算临界吸附温度tc，在tc附近其它温度的平均表面接触数<m>可以通过模拟数据的二次插值计算来获取。通过以上方法，计算得到了在外力f＝0时高分子链的临界吸附温度tc＝1.95和指数其临界吸附温度大于未采用协同运动算法所得到的tc＝1.65，两者差异的原因是高分子模型中采用不同的运动方式。

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

图3中展示了高分子链的平均表面接触数<m>与链长n的双对数关系，其中外力f＝0.3，温度为t＝3，2，1.6，1.2，1和0.5。从图3中，观察到受外力驱动作用下高分子链的平均表面接触数<m>与链长n之间都呈线性关系，且不同温度下的曲线几乎平行。由此可知，采用有限尺寸标度方法不能计算高分子链在外力驱动作用下的临界吸附温度tc，即不同温度下高分子链在<m>与n的关系中具有相同的标度指数。因此，下面将主要以较长的链长n＝200的高分子链为研究对象。

链长为n＝200的高分子链平均表面接触数<m>与温度t的关系，如图4所示，其中外力f＝0,0.1,0.3,0.5,1和10。高温时，高分子链平均表面接触数<m>→0。随温度t降低<m>发生伪相变，且其伪相变温度随着外力f增大而减小。当温度t→0时，高分子链几乎全部被吸附到了表面上，即<m>→n。

在<m>-t的关系中，对高分子链的表面接触数的中间数据进行线性拟合，其延长线与横轴的交叉点来粗略估计在外力f驱动作用下的临界吸附温度tc。得到tc与f的关系，如图5所示。tc随着f的增大而减小，但是当f>f^*时临界吸附温度tc不再变化，即tc^*＝0.41，其中f^*≈2(目测估计值)。由外力f和温度t构造了高分子链从脱附状态(desorbedstate，标记为ds)到吸附状态(adsorbedstate，标记为as)的伪相图，如图5所示，其插图的三种构象分别为(a)吸附态f＝0.1,t＝1.0，(b)脱附态f＝10,t＝1.0和(c)吸附态f＝10,t＝0.2。

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

高分子链在表面的吸附会伴随着其构象的变化，均方回转半径<rg²>可以描述高分子链的构象，其计算公式为：

其中，n为高分子链的链长，为高分子链中第i个单体的位置矢量，为高分子的质心位置矢量及其计算公式为：

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

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

高分子链在表面的吸附性质与温度相关，根据图5的伪相图中选择两个温度t＝1.2(>tc^*)和0.2(<tc^*)说明高分子链吸附和构象受外力f驱动的不同作用。高分子链的平均表面接触数<m>和均方回转半径在z方向上的分量<rg²>z与外力f的关系，如图8所示。当t＝1.2时，<m>随外力f增大而减小，而且其<m>值逐渐趋于0，即高分子链发生了由吸附状态向脱附状态的转变；当t＝0.2时，<m>随外力f增大而减小，但<m>仍然接近n，即高分子链中虽存在部分单体脱附但整条链仍处于吸附状态。当f＝0时，当t＝1.2和0.2时<rg²>z→0，高分子链处于接近表面的吸附状态。而随着f逐渐增大的过程中，温度t＝1.2时<rg²>z单调递增直到f>0.4后保持不变，此时<rg²>z明显大于0，其中fc＝0.4与图5中的相应值一致；而温度t＝0.2时即使外力f>2后<rg²>z仍然接近于0，其值明显小于温度t＝1.2时的<rg²>z，说明高分子链仍然在吸引表面附近。因此，在温度t＝1.2时，当外力f足够大使高分子链发生从吸附状态到脱附状态的相变，而在温度t＝0.2时随外力增加没有发生脱附相变。在两个不同温度t>tc^*和t<tc^*，高分子链吸附性质和构象性质受外力驱动作用的不同影响，其性质与图4的伪相图相一致。

通过本发明的方法，发现高分子链的临界吸附温度tc随外力f增大而减小，据此构建了链长n＝200的高分子链的脱附状态和吸附状态相对于温度t和外力f的伪相图。我们从高分子链的构象角度交叉校验了高分子链的临界吸附点，发现从高分子链的<rg²>的极小值与从其y和z分量的变化估计的临界吸附点相一致。最后讨论了温度t>tc^*和t<tc^*两种情况的高分子链吸附性质和构象性质受外力驱动作用的不同现象，当tc^*<t<tc时处于吸附态的高分子链随着外力f驱动作用的增加发生脱附，然而当t<tc^*时高分子链不会受外力f驱动作用发生脱附现象，这与高分子链吸附和脱附的伪相图性质相一致。

以上所揭露的仅为本发明较佳实施例而已，当然不能以此来限定本发明之权利范围，因此依本发明权利要求所作的等同变化，仍属本发明所涵盖的范围。