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

分 子 模 拟 牛继南 njn0516@cumt.edu.cn 2011.3 3.1 周期性边界条件 执行分子模拟时,通常选取一定数目N的 分子,将其放在一个立方状的盒子里. 设盒子 边长为L,则其体积为V=L3，若分子质量为 m, 那么盒子的密度为: 第三章 周期性边界条件和静电作用求和 3 L Nm =ρ 体系的密度要等于实际实验中所测得的 密度. 例如, 现在要进行1000个水分子的分子 模拟, 已知水的密度为1g/cm3，则可以求出 需要构建的盒子大小: Α=×=⇒ ×× = − 04.311004.31 )1002.6/18(1000 /1 8 3 23 3 cmL L g cmg • 以二维简单系统为例,如图所示为盒中粒子的排列 和移动方向. 其中位于中间的是要模拟的盒子,周 围的盒子与它具有相同的排列和运动, 被称为周期 性镜像(periodic mirror image)系统. 当计算系统中一 个粒子运动到外部,则必有另一个粒子从相对的方 向进入, 如粒子2. 这样的限制条件就使得系统中的 粒子数目保持恒定,密度不变. 实际计算当中,为了使系统的 密度保持恒定,通常采用在周 围增加同样的模拟单元构成 周期性体系的 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 来处理: • 在计算的盒子周围增加同样的模拟单元 （镜像）构成周期性体系来维持计算的进 行，使得最内部的盒子更贴近于实际的内 部液体或固体结构，此即为周期性边界条 件(periodic boundary condition). 周期性边界条件不仅仅是为了保证粒子 数平衡,它的引入也使得这个模型更贴近实 际. 模拟的中心盒子实际上是对现实物质世 界的一个简化模型, 从宏观结构中选取的一 个典型单元. • 分子模拟的核心是计 算分子或原子上受到 的力, 那么处于周期性 边界条件中的这些粒 子该如何计算??? 以粒子1为例: 先考虑最原始的盒子,必须计算它和2,3,4,5等 粒子的作用. 此外,我们还需要计算它和镜像A,B,C,D,E,F… 中所有粒子的作用. 整个系统是周期性的,因而这些作用项是无限 的, 所以直接计算是不可能的! 根据周期性边界条件,可以有下面两种方法 计算系统中分子间作用力,: 1，最小镜像变换 2，直接镜像变换 最小镜像变换(Minimum image convention) • 给定一个和原始盒子 尺寸一样的范围，让 每个粒子只和这范围 内的粒子反应一次. • 从图上看原始盒子只 和周围最近的这些镜 像盒子反应,因此又称 为最近邻镜像变换. • 为了提高计算速度,进 一步提出了,球形截距, 减少了反应的粒子数. • 在最小镜像变换中,势能的截距不能超过盒 子尺寸的一半,否则势能就会发生重叠. • 最小镜像变换的核心思想是每个盒子都是 等同的,一个盒子受到力的作用相当于每个 盒子都受到作用. • 为了保证计算的正确,盒子的尺寸要大于所 用力场截距的二倍,所以这种方法一般用来 计算较大的体系.有些要求不高的,对库伦能 也进行截断处理. • 缺点:不利于库伦力的精确计算. 直接镜像变换(Explicit image convention) • 不再限制截断半径的 范围,原始盒子中的每 个粒子都能和截距内 的其它所有粒子作用. 截距可以小于盒子长 度,也可大于盒子长度. • 相当于直接和周围截 断半径内的镜像盒子 作用,因此称为直接镜 像变换. 直接镜像变换的截断半径可以取的很大, 因此可以具有较高的计算精度.但随着向外 的扩展包含的盒子也变多,计算的粒子数目 也增多,计算代价增高. 最小镜像变换和直接镜像变换的异同 • 共同点: 都是仅考虑截距内的其它粒子对原始盒子 内粒子的作用, 然后将原始盒子的粒子状况镜像到 其它盒子. • 不同点: 和周围镜像盒子的关系不同. 最小镜像变换中,原始盒子和周围盒子都是"自己", 给一个盒子中粒子一个力,其它盒子上这个粒子的 镜像也会有同样的力, 因为都是一个"人". 直接镜像变换中,原始盒子和周围盒子"谁是谁", 原 始盒子粒子可以接受周围盒子中这个粒子的镜像 给予的力,不过周围盒子老跟着原始盒子"学". 3.2 截断半径 • 为了减少计算量, 在精度允许的情况下通常 会对势能进行截断, 如前面的两种镜像变换 处理中都用到了截断半径.通常的截断方法 有以下几种: 1. 直接截断 2. 利用开关函数 3. 利用常量移动 直接截断 • 即直接选取径向截断半径.通常选取rc时,有 两点原则(以下面这张vdw势能的典型曲线为 例 ): 1.势能接近于0; 2.为了提高计算速度, 截断半径在合理情况下 尽量缩短. • 但是势能在截断半径 处并不真正的为0，在 一些较为精确的计算 当中, 如果直接取为0 会造成能量不连续,会 形成假的低势能点. 为 了解决这种情况, 通常 引入开关函数 (switching function). 利用开关函数 r < Rs r > Rc sw函数 利用常量移动 • 除了利用开关函数,还可以用常量移动的方 法来减小势能面的断层. 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 P o t e n t i a l E n e r g y ( K c a l / m o l ) r (A) U实际=U原来-U(rc) rc 缺点:导数不连续 3.3 静电力的计算 在一个空间中, 原子间作用力随距离下降的 速度不大于r--d(指这个空间的维数)时, 这个 作用力被称为长程力.(例如三维空间d=3) 按照这个分类, 库伦力是典型的长程力. 计算机模拟时, 这些库伦力将是一个很大的 问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 , 因为它们的作用范围将远远大于盒子 的宽度的一半(以典型的含500个分子的体系 而言). ij q r qq U 0 21 4πε = • 例如: 在空间中两个带+0.5e和-0.5e的两个电 荷, 它们的之间的库伦能为: )Kcal/mol( 5.05.0 332 332 )J/mol( 4 21 0 21 ij ij A ij q r r zz N r qq U ×− = = = πε 0 200 400 600 800 1000 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 C o u l o m b E n e r g y ( K c a l / m o l ) r ij (angstrom) rij=100Å Uq=0.83Kcal/mol rij=1000Å Uq=0.083Kcal/mol • 通常的要求是: 库伦能要减少到小数点后第 二位才算是可以忽略, 但这时的盒子的长度 要宽于1000×2=2000Å，是极不现实的! • 原因: 通常的盒子长度在几十Å左右, 2000Å 的盒子将包括极大的原子数目, 已经远远超 出了大多数现代计算机的范围, 即便是超级 计算机!!! • time≈N2 or L6 解决方法 1. Ewald求和 2. 作用场(Reaction field) 3. 溶剂介电模型 直接的势能截断 × Ewald求和 Ewald求和法为Ewald在1921年 所提出,最初是用来计算离子 晶体的能量.这种方法的是先 选定一个原始盒子, 盒子中的 粒子除了与盒子中的其它粒 子作用外, 还和所有镜像中的 粒子作用. 按照周围盒子的顺 序, 在镜像盒子的数目达到极 限时, 整个区域如同一球体, 称为Ewald球体. Ewald求和就是建立在Ewald球体的基础之上(收敛要求), 对周围的介电常数也有要求. • 设原始盒子(中心盒子)为立方体, 边长为L, 含有N个粒子, 则镜像系统盒子的位置可以 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示为(±iL,±jL,±kL),其中(i,j,k=0,1,2,3…). • 这时原始盒子中所有粒子间的库伦作用为: • 与原始盒子最近的盒子共有6个, 其位置 分别为(L,0,0), (-L,0,0), (0,L,0), (0,-L,0), (0,0,L), (0,0,-L), 原始盒子与这6个盒子中粒 子的库伦作用为: • 周围镜像盒子的位置向量可以用 来表示, 其中nx,ny,nz均为整数, 则以此类推,可 以计算原始盒子和周围所有盒子中粒子的 库伦作用: • 如果加上原始盒子中所有粒子自身的相互 作用, 则上式可以变为: 其中, 第一个加和号右上角的一撇表示在 时, 也就是在原始盒子 中, 排除了粒子自身与自身的反应. 但是如果按照这个式子,对库伦能进行加和,则收敛的速度非常非常慢, how do it? • 对库仑项进行拉普拉斯变换，这时库仑项 分成两项，一项可以在实空间内快速收敛 ，其表达式为: real • erfc()为互补误差函数, 其表达式为: 误差函数的收敛速度非常快, 以至于库伦求和只需计算几个|n|的值就可以收敛. 求和的 速度取决于式子中的α, 其值越大, 收敛速度越快. 在实际计算中, 实际上只需选择较大 的α,值, 可使|n|=0, 即只需计算原始晶胞内的电荷作用. • 第二项可以在倒易空间进行收敛, 其表达式 为: 式中 为原始盒子内粒子的倒易矢量, 表达式为 这个式子的收敛速度也远快于原始的库伦加和式, 但其中的 正比于α, 因此需要和 实空间的计算速度进行平衡.通常选取α=5/L, 这时约有100~200个 的加成项. rep • 此外, 计算倒易空间加和的时间并没有扣除 电荷自身的作用, 因此需要额外扣除: self 此项不依赖于原子的位置, 因此给定原子后, 是一个常数. • 另外一项与环绕Ewald球的介质性质有关,若 介质外为导体,则不要此项, 如果是真空,则 需要加入此项: • 因此, 利用Ewald法计算长程库伦势能的表达 式为: Ewald加和的物理意义 • Ewald方法是各种计算静电库伦作用方法中 最为可行的方法, 通常用于高带电性的体系. 含Ewald加和的计算要比一般的分子模拟计 算要费时, 因为要计算很多加成项, 其速度 约为无静电作用时的十分之一. • 此外, 人们为了提高Ewald的计算速度, 还由 此改进出了PME(particle-mesh ewald), PPPM(particle-particle particle-mash Ewald)等 速度较高的静电加和方法. 作用场方法 • 作用场(reaction field)是一种简便的计算库伦能的 方法, 此方法计算截断半径内所有电荷间的库伦静 电作用, 而将截断半径外视为介电常数为εs的均匀 介质, 介质将会在截断半径的内部形成电场, 通过 计算出电场对内部分子偶极矩的作用, 就可以得到 分子的势能. 作用场方法的缺点是如果截断半径内的分子数目发生改 变,则计算的能量会出现不连续的情况.要克服这种情况, 经常用开关函数来控制截断半径附近情况, 避免能量出 现这种情况. 按照作用场方法, 均匀电介质所形成的电势为: 其中, 为作用半径内围着的i分子的偶极. 图中 黑色粗箭头为这些分子电偶极的矢量和, 因此i分 子所受作用场的作用势则可以表示为: i 名词解释: 偶极矩 • 分子是由带正电荷的原子核和带负电荷的电子所 组成, 正负电荷的总值相等, 所以整个分子呈电中 性, 但电荷中心可以重合也可以不重合. 前者称为 非极性分子, 如H2，CH4，C6H6等, 后者称为极性 分子,如HCl, H2O等. • 分子的极性通常用偶极矩来表示, 偶极矩的概念由 德拜在1912年提出: 两个带电荷分别为+q和-q的质 点, 相距d, 则之间的偶极矩为:μ=qd • 偶极矩是一个矢量, 方向规定从正到负. 分子的偶 极矩为所有键的偶极矩的矢量和. 溶剂介电模型 • 为了简化溶液系统中带电质点间相互的库 伦作用, Smith和Pettitt提出了所谓的溶剂介 电模型. 在这个模型中, 介电常数ε是一个有 效介电常数, 其大小和质点间的距离有关. 距离越近, 溶剂分子的效应越低, 介电常数 值与真空介电常数值越接近, 距离越远, 溶 剂分子的效应越显著, 介电常数接近与溶液 体系的介电常数. • 其中, εr是溶液体系的介电常数, εeff为有效介电常 数, 最小值为1，S为参数, 随溶液而不同, 通常的取 值范围在0.15~0.3Å之间. 右图为水溶液的有效介电常数的变化, 当距离小于2Å时, 有效介电常数开始向1靠近, 当距离大于30Å时, 有效介 电常数趋近于一般的介电常数值80. • 这种与距离相关的介电常数模型适用于像水溶液 等介电常数很大的系统, 因为介电常数越大, 带电 质点间的库伦作用也越小, 远程力也越不重要. 例 如分别带±0.5e电量的两个质点, 当在水溶液中相 距30Å时, 之间的库伦作用势为: 因此利用这个方法可以大量简化系统中远程作用力的计算, 如从上图可以看到, 当质点间距离超过15Å以后, 有效介电常数已经超过70，它们之间的库伦作用 已经显得非常微弱. 但需要注意的是: 这种方法缺乏相应的理论依据, 除非系统过于庞大, 否则 一般不要轻易使用. 3.4 作用列表的生成 • 我们构建的初始模型当中, 不同的原子之间 有着不同的作用, 那么这些作用计算机是如 何识别的呢? 成键作用 非键作用 成键作用 • 原子间如果存在共价键作用, 则这些原子间 的作用靠如下两种方式来进行识别: 1.通过事先给定的连接来生成列表(这个列表 在计算过程中是固定的) 2.在计算开始的时候搜索一定范围来自动生 成列表(这个列表在生成可能会发生改变) 非键作用 • 主要是通过直接搜索截距内原子生成列表. 两点需要注意: 1.应排除已经成键的原子(或者说成1-2,1-3甚 至1-4反应的原子) 2.计算开始生成的列表会在过程中不断变化( 通常是在每隔几步的间隔就升级一次列表) 上次课补充 • 描述反应的方式有以下几个: 1.给出反应列表 2.原子序数 3.元素名称 4.原子名字 5.原子电荷 主要取决于所用的软件的功能, 这在一定程度上甚至也决定了软件的易用性