物理问题的计算机模拟方法(1)—分子动力学.doc

硕士研究生课程物理问题的计算机模拟方法讲义适用专业: 凝聚态物理、材料物理与化学、理论物理、光学工程学 时：3040学时参考教材： 1德D.W.Heermann 著，秦克诚译，理论物理中的计算机模拟方法，北京大学出版社，1996。 2荷 Frenkel & Smit 著，汪文川 等译，分子模拟从算法到应用，化学工业出版社，2002。 3M.P.Allen and D.J.Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1989. 4. A.R.Leach, Molecular Modelling: Principles and Applications, Addison Wesley Longman, England, 1996. 5. 德 D.罗伯 著，计算材料学，化学工业出版社，2002。 6. 英 B. Chopard & Michel Droz 著，物理系统的元胞自动机模拟，祝玉学，赵学龙 译，清华大学出版社，2003。目 录第一章 计算机模拟方法概论1.1 序言1.2 热力学系统物理量的统计平均1.3 分子动力学方法模拟的基本思想1.4 蒙特卡罗方法模拟的基本思想1.5 元胞自动机模拟的基本思想1.5.1 简要的发展历程1.5.2 简单元胞自动机：奇偶规则1.5.3 元胞自动机的一般定义第二章 确定性模拟方法分子动力学方法（MD）2.1 分子动力学方法2.2 微正则系综分子动力学方法2.3 正则系综分子动力学方法2.4 等温等压系综分子动力学方法第三章 随机性模拟方法蒙特卡罗方法（MC）3.1 预备知识3.2 布朗动力学（BD）3.3 蒙特卡罗方法3.4 微正则系综蒙特卡罗方法3.5 正则系综蒙特卡罗方法3.6 等温等压系综蒙特卡罗方法3.7 巨正则系综蒙特卡罗方法第四章 离散性模拟方法原胞自动机（CA）4.1 引言4.2 元胞自动机模拟*4.3元胞自动机模拟的应用第一章 计算机模拟方法概论 1.1 序言1 什么是计算机模拟？ Simulation Modelling2为什么要进行计算机模拟？3常用的计算机模拟方法确定性模拟方法：MD模拟 Molecular Dynamics随机性模拟方法：MC模拟 Monte Carlo离散性模拟方法：CA模拟 Cellular Automata 1.2 热力学系统物理量的统计平均描述系统的坐标（自由度）：x(t)=x1(t),x2(t),xN(t)系统的物理量：A(x(t)1时间平均 分子动力学（MD）模拟 (1-1)2系综平均 蒙特卡罗(MC)模拟 (1-2) 分布函数（几率密度函数） (1-3) 配分函数 (1-4) 相空间 H(x)系统的哈密顿函数对于处于平衡态的系统，可以证明： 对于实际的有限时间内的平均，则有 实际模拟的系统大小也是有限的：有限的粒子数N或有限的系统限度L对统计平均结果有影响。 1.3 分子动力学(MD)方法模拟的基本思想1 基本原理系统：N个粒子，体积V，粒子质量为m描述一个粒子运动状态的自由度：（ri, pi） （pi=mvi）相空间：6N维，相空间中的一点的坐标 XN=rN, (mvN) rN=(r1, r2, , rN), vN=(v1, v2, , vN)粒子间的相互作用势：U(rN)=U(r1, r2, , rN)=决定系统相轨迹XN（t）的运动方程： (1-5)物理量A的宏观值，由A(XN) 的时间平均获得，即 （离散情况：）对于平衡态： 实际模拟时间总是有限的，模拟时间的长短可通过判断时间的增加对平均值的影响来确定，当继续增加时间带来的平均值得变化在允许的误差范围之内时，即可认为模拟足够长了。2 计算步骤运动方程： 即 (1-6) 或 (1-7)数值求解：用差分近似表示微分 （采用不同的差分格式，可得到不同的算法）。用显示中心差分格式，将（7）式写为 (1-8)由（7）和（8）式可得： (1-9)第一步：由（9）式计算第i个粒子在t+t时刻的位置坐标。 要启动计算，我们必须要知道最初两点ri(0)和ri(t)第二步：对不同时刻 t = t , 2t, 3t, , Lt (t0 = 0) 计算物理量A(r1(lt), r2(lt), , rN(lt) (l = 1, 2, , N)第三步：计算物理量A的平均值 L的大小由继续增大L 而不变（或变化在误差范围内）来确定。 1.4 蒙特卡罗（MC）方法模拟的基本思想1 基本原理以正则系综（T, V, N）为例正则分布：正则配分函数：系统能量：物理量：A(rN) = A(r1, r1, , rN )系综平均： （1-10） （位形积分） （1-11）用MC方法计算上述多维积分。2 计算步骤（1） 划分原胞 N个粒子3N个空间自由度，3N维空间划分成s个相等的原胞（s1）注意：由于积分中不含动量，所以我们只需要在位置空间积分，而不需要在相空间中积分。当系统的代表点落入第i个原胞时，则认为系统处在状态i，因此，s为系统可能的微观状态数目。于是，积分（10）和（11）可近似表述为 (1-12) (1-13)（2） 建立马尔可夫(Mako)过程（链）将s个状态看作一组随机事件马尔可夫链：从状态ij状态j(eiej)的概率pij ，只与ei和ej 有关。 ，i=1, 2, , s若ei 经历n步到达ej ，其概率表示为，存在极限概率 ( j=1, 2, , s)uj为系统处在状态j的概率。于是，沿无限长的马尔可夫链，物理量A的平均值可写为 (1-14) 选取 ，则（14）式为A的正则系综平均值。（3）抽样方法 采用怎样的抽样方法所构成的马尔可夫链能得到上述平均值？ 粒子位置坐标：粒子编号：g =1, 2, N坐标的三个方向：a = 1, 2, 3 系统状态：i = 1, 2, , s给定粒子位置坐标的变化量d (小于系统体积的限度)给定系统的初态i，随机选定4个随机数，其中三个xa （a =1， 2，3），且 1 xa 1，一个表示粒子编号 g =1, 2, N，由此随机确定粒子位置的变化： （ 确保 ）若，则 g 运动到新的位置，即系统由状态i过渡到状态j；若，则再选一个随机数x 4（0 x 4 1），若，则粒子保留在原位置，不发生ij的跃迁；若，则发生i j的跃迁。由此进行下去，则形成一个马尔可夫链（或过程），此链的长度L(即粒子行走的步数，远大于s)，由所计算的物理量的平均值 (1-15)不再随链的加长而改变来确定。由此得到的平均值即正则平均值。一般来说，L与N, V, T 有关，比如，N=32108, L=30005000。归纳起来，计算系统物理量的正则系综平均值的具体步骤如下：第一步：给定系统的初始状态(粒子的初始位置)ri和每一步的改变量d；第二步：选择四个随机数，其中一个代表粒子的编号i ( 1 i N )；另外三个表示粒子空间坐标的改变dx, dy, dz ( -d da d , a=1, 2, 3);第三步：计算粒子i的新位置 第四步：计算粒子在新旧两个位置系统的能量之差 第五步：由的大小判断粒子i是否从ri运动到: 若，则ri； 若，则再选一个随机数R(0 R rc 的相互作用忽略不计，代价是忽略了背景。5 积分格式 于是，运动方程可写为 (2-1) （1）递推公式在第一章中求解运动方程时，我们是直接求解的关于粒子位置坐标的二阶微分方程，得到的递推公式需要知道最初的两点位置才能启动计算，但在实际计算中，我们常常是给出最初的位置和速度，于是，我们可通过选取一定的差分格式，有运动方程(16)得到关于位置和速度的递推公式。 或 已知递推公式的具体形式取决于差分各式的选取。（2）时间步长选择时间步长的原则：在保证计算精度的前提下，尽量节省计算时间。实例：Ar原子系统，LJ势 时间步长取为 (=100ps=10fs)（3）约束条件保证能量、动量或角动量守恒。（4）减小计算误差的技巧数值计算不可避免有误差，与误差有关的因素主要有： 差分格式 时间步长 切断半径 最小影像约定 等等6. 计算热力学量 （1） 若系统的NVE恒定，则 （微正则系综平均microcannonical ensemble average）（2） 若系统的NVT恒定，则 （正则系综平均cannonical ensemble average）在实际计算中，往往还需要要求系统的总动量P=0或恒定，因整个系统不受外界的作用。（3） 温度与能量均分定理 N总粒子数，Nc约束条件个数, 一般情况下，N Nc P=0 Nc = 3（4） 位势切断误差与修正g(r) 对关联函数（pair correlation function） 粒子数密度当原点位置上一个粒子时，在r 附近dr内发现有一个粒子的几率。对气体或液体， （各向同性）平均每个粒子的能量切断误差修正为：6 计算机模拟的组织（1） 初始化给定粒子的初始位置、初始速度等；（2） 趋衡从初始转台开始，按运动方程要求的规律，从非平衡态趋向平衡态；（3） 投产计算物理量的统计平均值。 2.2 微正则系综分子动力学方法微正则系综：NVE恒定，且P = 0 （分子动力学方法的自然选择）1Verlet 算法运动方程的显示中心差分格式由可得到如下递推公式 （2-2）此即Verlet 算法（A2）, 其特点是：（1） 要启动此算法，必须已知两点，所以又称为二步法；（2） 和不能同步算出，第n+1步算出的第n步的速度。 2Verlet 算法的速度形式 由此可得如下递推公式： （2-3）此即Verlet 算法的速度形式（A3）, 其特点是：（1） 启动此算法，必须已知两点；（2）和同步计算。（的计算需要知道）此算法比A2算法更稳定。3趋恒阶段的能量调整由于系统初始状态的能量离我们要模拟系统的能量（或温度）有一定差异，这就需要在系统趋于平衡的阶段进行调整。最常用的方法之一就是进行速度标度： 标度因子： （2-4）Tref 为设定的系统的温度值。 能量设定值： 能量参考值： 要注意的问题：（1） 一般不需要每一步都进行标度，可每隔若干步（比如50步）调整一次；（2） 当能量达到所设定的值后，停止标度，在以下的模拟时间内能量保持不变。4 实例氩原子系统， N=256LJ势:作用力： 取时间单位：取长度单位：s = 0.3405nm取质量单位：m=6.6338210-26 kg(氩原子质量)使用上述单位，可得到约化单位下的作用力表达式：标度因子： 约化温度： 取时间步长h = 210-14s = 20fs, 相当于约化时间0.064约化温度：T* = 2.53 303K, T* = 0.722 86.5K约化数密度：r* = N/L3, N = 256r* = 0.636 L = 7.83 , r* = 0.83134 L = 6.75当 N = 64 时，r* = 0.83134 L = 4.25, 所以，取rc = 2.5 是错的，因为要求 rc L/2 = 2.13。（习题3.5）计算源程序见p129, 程序PL1 2.3 正则系综分子动力学方法正则系综：NVT恒定，且P = 0 如何保持温度T恒定？1速度标度方法 算法A4:(1) 规定初始位置和初始速度；(2) 计算；(3) 计算；(4) 对速度进行标度：，返回（2）。 说明：（1） 在前面的微正则系综的模拟中也对速度进行了标度，但其目的是使总能量达到所需要得值，而且不必每步进行速度标度，可隔若干步标度一次，一旦总能量达到所设定的值，即可停止标度；而在正则系综的模拟中，每步都必须进行标度，以始终保持温度恒定。（2） 通过一个广义位势引进能量涨落，连同细致耦合的一种特殊选择，会导致速度标度机制。（见书中的证明，公式3.59有误）（3） 由于控制温度的反馈环内的时间延迟，温度将在一定程度上涨落。2随机方法（Anderson 热浴）体系与一强加了温度的热浴相耦合，与热浴的耦合由偶尔作用于随机选择的粒子上的脉冲力表示。在开始模拟前首先要选择系统与热浴的耦合强度，这一耦合强度由随机碰撞的频率n 决定。模拟步骤：（1）规定初始位置和初始速度；（2）计算；（3）粒子i在时间间隔h内发生碰撞的几率为hn，产生随机数x，如果x hn, 则粒子i 经历一次随即碰撞，否则，返回（2）；（4）粒子i 从温度为T的麦克斯韦玻尔兹曼分布中获得一个新的速度，返回（2）。 2.4 等温等压系综分子动力学方法1等压等焓系综 N P H恒定，P = 0 (P:压强，P：总动量) (PE为外压强，平衡时，PE =P) 若系统与外界绝热，则 ，即（等焓过程）要保持压强不变，体积则必须变化，为此，将体积V作为一个新的动力学变量，构造一个新的拉氏函数： 引入标度坐标： 或 （，体积变化缓慢） ， （共轭动量）系统的哈密顿量为： （2-5）运动方程为： （对粒子） （2-6） （对活塞） （2-7）将（2-5）代入（2-6）得所以， （2-8）将（2-5）代入（2-7）得所以， （2-9）由维里定理 (2-10)归纳起来：这些方程给出一个恒定的平均压强。下面推出递推公式，利用泰勒展开保留到二阶项得： 将（2-8）和（2-9）式代入上式得： （2-11）定义偏速