┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ 装 ┊ ┊ ┊ ┊ ┊ 订 ┊ ┊ ┊ ┊ ┊ 线 ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊
二甲醚在橡胶中扩散特性的分子动力学模拟
?1???????mkvk.vk?，用来保持总动能守恒，为了把温度调到期v.F上式中?G???jj????j??k望的温度下，需要把初始的温度设定为期望值。 (3)Nose-Hoover方法
1984年Nose提出了“Nose动力学”，这种方法最有力的方法是扩展体系Hamilton 量的做法，不但能够用于恒温体系，而且可以启发延用到恒压体系或恒温一恒压体系。 Nose动力学的直接结果是实现了正则分布。之后，Hoover对这个方法进行了修正，在运动方程中引入了一个摩擦项和热储项来延展系统的哈密顿量，摩擦力和每个粒子速度和摩擦参数?（热浴参数）的乘积成比例，这个摩擦参数是一个动态量。Hoover的运动方程能够用“实时间”的等间距抽样得到正则分布，这对非平衡的模拟显得更加重要。
Hoover方程中，粒子的运动方程被替代为:
d2riFidr??i (2-20) 2?dtmidt其中，热浴参数?的运动方程是：
d?1??T?T0? (2-21) dtQT为系统的当前瞬时温度，T0为参考温度。藕合的力度由常数Q(储热的“质量参数”)和参考温度共同决定。
用质量参数描述耦合强度有些困难，为了保持藕合强度，可以使Q和T0相 联系起来:
2?TT0 Q? (2-22) 24??T是系统和热储之间的动能振荡周期，它独立于系统的大小和参考温度。
（4）Berendsen弱耦合方法
设体系放置在一个温度为T0的大热浴内，体系温度为T(t)。体系和热浴之间的导热过程在唯象上服从Fourier定律:
dT?t?1??T0?T?t?? (2-23) d?t??比例系数是1/?,?的量纲是时间，叫做弛豫时间。导热弛豫时间越小表明导热速度越快。把Tt??t?2?在t-?t2处Taylor展开，再根据Fourier定律就可以得到:
???tT0T?t??t2??1???1????2 (2-24)
T?t??t2??T?t??t2???
第 18 页 共44页
┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ 装 ┊ ┊ ┊ ┊ ┊ 订 ┊ ┊ ┊ ┊ ┊ 线 ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊
化简后可得到
二甲醚在橡胶中扩散特性的分子动力学模拟
???tT0?T?T??t ??1?? (2-25) ?1??1??0??T???T?t??t2???或者
??T?t??t2?vi?t??t2? (2-26) ?T?t??t2?vi?t??t2??是标度因子，上式在推导过程中认为温度是体系粒子动能之和。即：
N31 NKBT?t???mivi?t??vi?t?
2i?12 (2-27)
当体系粒子速度的绝对值vi??t2偏小时，则相对应的温度Tt??t?2?低于T，
0从公式(2-25)可以得到?>1，把?>1代到公式(2-26)得到的vi??t于把Tt??t2就会变大，等同
?2?提高，反之亦然。所以通过这种方法可以调节温度。
Berendsen弱耦合算法把一个系统弛豫到目标温度的效果很好。但是当系统达到平衡后，弱耦合算法不能探索到一个正确的正则系综。 2.3.9能量最小化
为了使模拟计算结果更准确需要对分子结构进行优化。当一个无规则晶胞生成时，分子有可能不是等价地分布在晶胞中的，这样就会造成真空区。为了矫正这个，也要对晶胞优化。对分子的结构和晶胞进行优化的过程实质上就是对两者的能量最小化的过程。对于能量的优化有很多种方法，如最陡下降法、共轭梯度法、牛顿一拉斐孙方法等。
最陡下降法是最简单的一种结构优化算法。在最低能量搜索的过程中，最陡下降法对能量函数进行微分，计算梯度，每次沿能量下降最多的方向前进。当搜索位置离能量极小点比较远时，用这种方法能够迅速向极小点靠近。当接近极小点时，会产生振荡，收敛速度变慢。
共轭梯度法在最低能量搜索过程中和最陡下降法一样对能量函数进行微分并计算梯度，但是在选择搜索方向时，不但考虑当前的梯度，而且还考虑原来的搜索方向，经过综合决定下一步的搜索方向。共轭梯度法的收敛速度快，但是容易陷入能最局部极小点，而且构象能远离极小点时共轭梯度法稳定性较差。
牛顿一拉斐孙方法利用一阶微分确定搜索方向，同时用二阶微分确定梯度在什么地方改变方向。这种方法计算精度高，可以迅速收敛，但是计算量非常大。
最陡下降法和共轭梯度法可以用于优化较大的系统，但是计算精度较差。而牛顿
第 19 页 共44页
┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ 装 ┊ ┊ ┊ ┊ ┊ 订 ┊ ┊ ┊ ┊ ┊ 线 ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊
二甲醚在橡胶中扩散特性的分子动力学模拟
一拉斐孙方法不适合用于处理较大的系统。本研究使用的能量最小化方法Materials Studio中的Smart Minimizer方法，它集成了最陡下降法、共轭梯度法和牛顿法。在进行能量优化时，如果系统初始构型离能量极小点较远时，首先采用最陡下降法进行优化，使体系迅速向极小点靠近。当能量迭代达到1000kca1/mol时，采用共轭梯度法。最后当能量迭代小于l0kcal/mol时采用牛顿一拉斐孙方法，可以提高计算的精度。Smart Minimizer方法结合以上三种方法的优势，能够大大加快构型优化的速度和精度。
2.4扩散系数的计算
分子扩散系数是表征物质分子扩散能力的物理量。根据斐克定律，组分A在组分B中的分子扩散系数，等于该物质在单位时间内，单位浓度梯度下、经单位面积沿扩散方向传递的物质量。
分子动力学模拟(MD)近年来成为模拟传递过程的一种重要手段。分子动力学计算扩散系数有两种方法:一种是通过计算体系的各种速度相关函数并对其积分，即Green-Kubo法;一种是通过计算体系平均平方位移并对其微分，即Einsitein法。 （1）Green-Kubo法
从微观角度，自扩散系数可以通过粒子的速度相关函数获得
1???? D????vi?0????vi?t??dt （2-27）
30?i??i?如果扩散物质的浓度比较低，粒子之间的相互作用力为短程力，方程（2-27） 可以改写为方程（2-28）
??1 D???vi?0???vi?t??dt （2-28）
30公式（2-27），（2-28）即为平衡系统时间相关函数求得的扩散系数Green-Kubo式，其中vi?t?表示原子i在t时刻的速度。 （2）Einstein法
动力学将所模拟分子进行标记，计算标记粒子的浓度分布，假定在时间t=o, 标记粒子被集中在坐标原点。为了计算浓度分布的时间演绎，须联立Fick定律 和物质守恒方程得到:
?2c?r,t??D?2c?t??0 (2-29) 2?t在边界条件c?r,0????r?下，(2-29)得到
第 20 页 共44页
┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ 装 ┊ ┊ ┊ ┊ ┊ 订 ┊ ┊ ┊ ┊ ┊ 线 ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊
二甲醚在橡胶中扩散特性的分子动力学模拟
r2?t???c?r,t?r2dr (2-30)
在式(2-30)中，令
?c?r,t?r2dr?1 (2-31)
将式 (2-31)乘以r2并在全空间上积分，得到方程：
(2-32)左边等于
?2rc?r,t?dr?D?r2?2c?r,t?dr (2-32) ??t?r2?t??t。右边进行分布积分，得到(2-33)，它将扩散系数D与
浓度分布联系起来，D是宏观传递系数，r2?t?是微观解释。
?r2?t??t?D?r2?2c?r,t?dr?2dD (2-33)
对于每个粒子i，测量在时间t内经过一定的距离，得到均方距离对时间的函数 关系如式 (2-34)
?r?t?21n2???ri?t? (2-34) Nt?112ri?t??ri?0? (2-35) 6t16NtN基于Einstein方程得到无限稀释扩散系数为: D??limt??自扩散系数为：
Dself?limt???r?t??r?0?iii?12 (2-36)
式中极限的部分近似用均方位移对时间的变化率来代替，就是说用均方位移随时间变化的曲线斜率a来代替。因此上式简化为
D?a （2-37） 6理论上，通过Green-Kubo法和Einstein方法得到的扩散系数是相等的，但对于
不同的体系两种方法的使用范围有所不同。本文使用这两种算法分别计算二甲醚在橡胶内的的扩散系数。
第 21 页 共44页