中国矿业大学分子模拟课件第十一章

分子模拟牛继南[e-mail]2011.3第十一章分子动力学计算实例11.1水分子系统的计算水分子是最常见的溶剂分子，许多系统中都含有水分子，因此在分子动力学的发展过程中，水分子系统的计算一直是一个重点。原因有两点：一，有充分的实验资料可以检验计算的重要性，二，水分子本身具有多原子分子的运动模式，分子间又有特殊的氢键存在，若计算成功则可方法推广到其它复杂系统。（1）简单的水分子力场常见的水分子简单力场由以下几种：中心力力场centralforce，CF可转移的分子间势能transferableintermolecularpotential，TIP简单点电荷力场simplepointcharge，SPC中心力力场最早又Stillinger提出，其形式为：力场中水分子本身为刚性分子。•Stillinger和Rahman用中心力力场，对216个水分子系统（密度1g/cm3）在29.5°C条件下进行分子动力学计算，得到了H-H间的径向分布函数，结果发现，除了1.5Å处的峰外，2.25处也出现了明显的峰，这和实际的XRD结果十分吻合。可转移的分子间势能TIP力场是由Jorgensen利用蒙地卡罗方法和精确的量子力学计算结果比较所推导出来的。除了利用了各种热力学性质，最重要的是比较了水分子二聚体的结构。此力场将水分子视为刚体。力场形式为：a,b为不同水分子上的原子，O电荷-0.8，H电荷0.4，AO2=0.58kcal/mol，AH2=0，CO2=0.58kcal/mol，CH2=0.•此力场由于形式简单，并且不考虑分子中各种高频率的振动，因此通常与其它力场合并应用到各种体系中。•为了提高力场的分布函数性能，Jorgensen将原来力场（TIPS）中水分子的3个作用点改为4个（TIPS2），力场形式和原来形式相同。其参数为:O电荷0，H电荷0.535，质量中心的电荷-1.07，AO2=0,00695kcal/mol，AH2=AM2=0，CO2=600kcal/mol，CH2=CM2=0.此力场得到的线性二聚体结构能为-6.20kcal/mol,接近于实验值-5.44±0.7kcal/mol随后，Jorgensen等陆续又发展出了三作用点的TIP3P和四作用点的TIP4P力场。简单点电荷力场SPC力场是由Berendsen等所提出的，此力场也将水分子视为刚性分子。力场形式为：力场中仅仅O-O间有VDW作用，而H和其它分子间只有库伦作用。SPC力场可以获得合理的RDF分布和热力学性质，但扩散系数有较大误差（200%-300%）。•各个力场的参数：•其中应用最广泛的是TIP3P力场。一些著名软件如CHARMM和AMBER等都采用此力场。（2）abinition水分子力场Nieser、Corongiu和Clementi等在1990年利用精确的量子从头算方法，发展出了基于液体二聚体和三聚体的力场，简称NCC力场。此力场和TIP4P力场一样，H带负电荷，HH连线中心带正电。各原子定义如图：•力场形式为:利用NCC力场可精确计算水分子图案体的各种性质。•随后Corongiu改进了NCC力场，加入了分子内部振动项，形式为：式子中可以计算振动光谱（3）超额自由能的计算超额自由能（excessfreeenergy）的定义为真实溶液的自由能和理想的无分子间作用的溶液的自由能差。超额自由能可以由下式计算:U为系统势能，λ为耦合参数，耦合参数的意义为：其为1时对应液体分子受全力场的作用，为0时对应理想溶液的状态。•为了计算ΔA，将此积分分为Nsim个相等的区间（10-20个），对应于λ从0到1的变化，分别求出每一个λ的积分项求得ΔA。•式子中Δλ=1/Nsim.•以SPC力场为例，Quintana将其化为：•其对108个水分子，密度为1g/cm3，分别在298k和373k进行动力学计算。先执行λ=1的运算，然后逐步降低λ至0，降低的幅度为0.1，系统达到平衡后，再计算2ps求出的值。•结果和实验结果十分类似。11.2正戊烷构象的研究•聚合物分子的许多构象主要由围绕单键的二面角旋转所形成，研究构象时可以在0-360°范围内旋转每个二面角。但如果分子含n个二面角，每次转动Δφ，则共需要综合考察个构象。以含10个二面角、每次转动60°为例，共需考察610个构象，设每个构象需30s的研究时间，则总共需要60年，这样的研究是不切实际的。•实际当中通常旋转较大的角来减少考察次数。如正戊烷中有两个C-C-C-C二面角，每次旋转120°，则共有9个构象需要考察。•定义构型二面角呈180°为反式（t），呈60°为正邻位交叉式（g+），呈-60°为负邻位交叉式（g-）。从60°为起点进行旋转。•表中的构象有几个是等效的，如2和3,4和5,8和9，而2和4,3和5，6和7是对映构型，因此实际上这9个构型仅有4中不同的情形：•对于这四种构型，可先按照能量最小化的方法找出附近的最优能量结构。•Mencareli用MM+力场，SD+对角块NR方法求得的最优结构如下：•Mencareli进而利用这些结构作为初始构型，进行MD计算，他将分子动力学过程中的每一步都作为分子的一种构象，按照玻尔兹曼分布定理，平衡体系中任一构象出现的概率为：•式中为吉布斯自由能，Hi约等于构型能。•可以看出当温度升高到600K时，各种构象出现的概率与实验值非常接近，表明该温度的分子动力学计算可以成功探测正戊烷的所有构象异构体，在低温条件下（300K）虽然趋势一致，但分布不同，表示分子动力计算的时间不够长，未达到研究所有构象的要求。11.3链状分子的计算•烷类分子是最最常见的链状分子，液态正十六烷（n-hexadecane）是柴油的主要成分，在其中添加α-烷基萘可以提高柴油的流动性或降低其固化点。Ying等利用MS对正十六烷和两种添加剂的混合物进行研究，他们发现加入α-甲基萘后，扩散系数提高了3倍，与实验符合的很好。11.4链状双性体的分子动力计算链状双形体分子为一根或数根碳氢链脸上头端的极性官能团所形成的分子。极性头端具有强烈的亲水性，而碳氢烷链具有疏水性，因此这种分子被称为双性。•双性分子最重要的特性是可以形成延伸的层状结构，在水/空气界面形成单层的langmuir膜，在固体表面形成Langmuir-Blogdett膜（LB膜）。LB膜可以为单层或多层结构。•实验室中可以控制膜层的厚度和双性分子的排列，作为半导体材料的绝缘体、过滤组件、抗反射镀膜等。•双性分子在生物系统中叶有重要作用，如细胞膜即为双性的双层脂质分子组成，在水溶液中可形成微胞。•链状双性分子的排列通常由XRD或NMR实验测量。NMR测量的碳链秩序性可由秩序参数表示：式子中，θi为第i根碳链轴与参考轴的夹角，参考轴的方向为所有碳链的平均方向，通常设为z轴，δij为Kroneckerdelta参数，若i=j，其为1，若不等，则为0.•Essex利用联合原子模型对双性分子进行了动力学计算。由于联合原子模型对力场进行简化，因此对秩序参数的计算方法做一定的修正。他对分子轴向量的定义如下：z向量为通过Cn+1和Cn-1的向量y向量为垂直z轴且过Cn、Cn+1和Cn-1的向量x向量为垂直于y轴和z轴的向量则分子动力计算的秩序参数和NMR测得的持续参数的关系为：以这种方法，Essex获得了和实验一致的结果。•Kim等对LB薄膜进行了分子动力学计算，系统为石墨表面的LB膜。因为没有可用的力场，他根据实验数据设计出特定的势能函数，石墨表面的作用势与原子α离表面的距离zα相关，其形式为：•式中ρ为表面的数值密度，为VDW势能参数，双性分子头端与其镜像间的库伦作用为：•Kim等以全原子模型对64根链在石墨表面的LB进行了计算，计算的结果显示，双性分子约呈20°的倾斜。•他们通过计算还发现了这种倾斜是因为头端的面积增加所引起的。11.5聚合物的特殊运动•安全玻璃中含有大量的苯环，苯环可以绕着单键转动而不改变整体结构，所以一般认为安全玻璃受到撞击后系统所接受的能量被用作苯环的旋转。由于运动消耗了接受的能量，因此不会导致玻璃的破碎。•Schaefer等由固态NMR测量结果认为苯环的主要运动为180度的旋转，称为π翻转。但是对于苯环旋转消耗多少能量，不同的文献有不同的推测，一般认为翻转能量在3至12kcal/mol之间。•陈正隆首先提出用分子动力学方法研究安全玻璃中的苯环运动。左图为一个苯环扭转角的运动轨迹，此苯环的扭转角在13-28ps的时间范围内改变了180°，正是了π运动的发生。计算还显示在300k时仅有两个苯环发生翻转，但在410k时共有9个苯环发生π翻转。这表明温度越高，外界提供的能量越多，越容易发生π翻转。•为了定义苯环的翻转运动，对翻转坐标进行定义：•坐标原点为发生翻转的苯环的中心，z轴为沿着苯环键的方向，x轴位于L1-L2-C4平面上。•按照此坐标系，苯环绕z轴的角动量为：•p为沿角标方向的动量。•定义顺时针方向绕z轴的旋转为正值。•苯环旋转的动能为:•Izz为z方向的转动惯量，即：Tasi模拟的一个苯环翻转时的动能和角动量变化有正有负，不是向同一方向转动相反反向翻转的能量障碍•下表为不同温度下苯环翻转的时间和对应的翻转能量障碍。苯环的翻转能量障碍范围在2.8-12.8kcal/mol之间，与实验值符合的很好。