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

硕士研究生课程物理问题的计算机模拟方法讲义适用专业 :凝聚态物理、材料物理与化学、理论物理、光学工程学时：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 Modelling 2为什么要进行计算机模拟？3常用的计算机模拟方法确定性模拟方法： MD 模拟Molecular Dynamics 随机性模拟方法： MC 模拟Monte Carlo 离散性模拟方法： CA 模拟Cellular Automata 1.2 热力学系统物理量的统计平均描述系统的坐标（自由度） ：x(t)=x1(t),x2(t),xN(t) 系统的物理量： A(x(t) 1时间平均dttAttAttt0)(10 x 分子动力学（ MD）模拟(1-1) 2系综平均xxxxxxdAdHfAZA)()()()(1 蒙特卡罗 ( MC) 模拟(1-2) )(1)(xxHfZ 分布函数（几率密度函数）(1-3) xx dHfZ)( 配分函数(1-4) 相空间H(x)系统的哈密顿函数对于处于平衡态的系统，可以证明：AA对于实际的有限时间内的平均，则有AA实际模拟的系统大小也是有限的：有限的粒子数 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)=Njiiju)(r决定系统相轨迹 XN（t）的运动方程：()0(N.,2, 1)(,0）加上边界条件（周期性初始条件）（NNNiiiiXXiUdtdmdtdrvvr(1-5) 物理量 A 的宏观值，由 A(XN) 的时间平均获得，即tt dtAttA0)(1)(X（离散情况：kiiAktA11)(）对于平衡态：)(limtAAt实际模拟时间总是有限的， 模拟时间的长短可通过判断时间的增加对平均值的影响来确定，当继续增加时间带来的平均值得变化在允许的误差范围之内时，即可认为模拟足够长了。2计算步骤运动方程：)(,NiiiiUdtdmdtdrvvr即iNiiUdtdmFrr)(22(1-6) 或mdtdii/22Fr(1-7) 数值求解：用差分近似表示微分（采用不同的差分格式，可得到不同的算法） 。用显示中心差分格式，将（7）式写为222)/()()(2)(ttttttdtdiiiirrrr(1-8) 由（7）和（ 8）式可得：mttttttiiii/)()()(2)(2Frrr(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 的平均值LlNLtltltlLA121)(,),(),(A(1limrrrL 的大小由继续增大L 而不变（或变化在误差范围内）来确定。 1.4 蒙特卡罗（ MC ）方法模拟的基本思想1基本原理以正则系综 （T, V , N）为例正则分布：sEseZ1正则配分函数：)(!13NNENddddehNspr系统能量：)(2)(2NiiiNpsUmUEErpr物理量： A(rN) = A(r1, r1, , rN ) 系综平均：)(1)(!1)(!1)(!1/ )(/2/)(3/332NkTUNNNkTmNkTUNNNNkTENNsNNdeAQdedeAZhNddeAZhNdAhNANiiiNsrrprrprrrrpr（1-10）NkTUNdeQNrr/)(（位形积分）（1-11）用 MC 方法计算上述多维积分。2计算步骤（1） 划分原胞N 个粒子 3N 个空间自由度， 3N 维空间划分成 s个相等的原胞（ s1）注意：由于积分中不含动量，所以我们只需要在位置空间积分，而不需要在相空间中积分。当系统的代表点落入第i 个原胞时，则认为系统处在状态i，因此，s为系统可能的微观状态数目。于是，积分（ 10）和（ 11）可近似表述为ikTUiNieAQA/1(1-12) ikTUNieQ/(1-13) （2） 建立马尔可夫 (Mako)过程（链）将 s个状态看作一组随机事件马尔可夫链： 从状态 ij 状态 j (ij) 的概率 pij，只与i和j有关。1jijp，i=1, 2, , s若i经历 n 步到达j，其概率表示为)(nijp，存在极限概率jnijnup)(lim( j=1, 2, , s) jjjuu1, 0uj为系统处在状态 j 的概率。于是，沿无限长的马尔可夫链，物理量A 的平均值可写为)(ijijiiiipAuAA (1-14) 选取kTUNiieQu/1，则（ 14）式为 A 的正则系综平均值。（3）抽样方法采用怎样的抽样方法所构成的马尔可夫链能得到上述平均值？粒子位置坐标：)(ix粒子编号： =1, 2, N 坐标的三个方向： = 1, 2, 3 系统状态： i = 1, 2, , s给定粒子位置坐标的变化量 (小于系统体积的限度 ) 给定系统的初态 i，随机选定 4 个随机数，其中三个（ =1， 2，3） ，且 1 1，一个表示粒子编号 =1, 2, N，由此随机确定粒子位置的变化：)()()(ijixxx（ 确保)()(jixx）若ijUU，则运动到新的位置，即系统由状态i 过渡到状态 j；若ijUU，则再选一个随机数4（0 41） ，若kTUUije/)(4，则粒子保留在原位置，不发生ij 的跃迁；若kTUUije/ )(4，则发生 i j 的跃迁。由此进行下去，则形成一个马尔可夫链（或过程），此链的长度L(即粒子行走的步数，远大于 s)，由所计算的物理量的平均值LiiLALA11lim(1-15) 不再随链的加长而改变来确定。由此得到的平均值即正则平均值。一般来说，L 与 N, V, T 有关，比如， N=32108, L=30005000。归纳起来，计算系统物理量的正则系综平均值的具体步骤如下：第一步：给定系统的初始状态(粒子的初始位置 )ri和每一步的改变量；第二步：选择四个随机数， 其中一个代表粒子的编号i ( 1 i N )；另外三个表示粒子空间坐标的改变x, y, z ( - , =1, 2, 3); 第三步：计算粒子i 的新位置ir),(),(ziyixiiiiizyxzyxr第四步：计算粒子在新旧两个位置系统的能量之差),.,.,(),.,.,(2121NiNiUUUrrrrrrrr第五步：由U 的大小判断粒子 i 是否从 ri运动到ir: 若0U，则 riir；若0U，则再选一个随机数R(0 R rc的相互作用忽略不计， 代价是忽略了背景。5积分格式jiijijjiiNiirruU)()()()()(FrF于是，运动方程可写为jiijiiirdtdmdtd)()(,Fvvr(2-1) （1）递推公式在第一章中求解运动方程时，我们是直接求解的关于粒子位置坐标的二阶微分方程，得到的递推公式需要知道最初的两点位置才能启动计算，但在实际计算中，我们常常是给出最初的位置和速度，于是，我们可通过选取一定的差分格式，有运动方程(16)得到关于位置和速度的递推公式。)1()()1()(ninininivvrr或)()()()(htvtvhtrtr已知00)0()0(vvrr递推公式的具体形式取决于差分各式的选取。（2）时间步长选择时间步长的原则：在保证计算精度的前提下，尽量节省计算时间。实例： Ar 原子系统， LJ势时间步长取为sh14210)(10约化时间步长(=100ps=10fs) （3）约束条件保证能量、动量或角动量守恒。（4）减小计算误差的技巧数值计算不可避免有误差，与误差有关的因素主要有：差分格式时间步长切断半径最小影像约定等等6. 计算热力学量dttVttAttAttNt)(),(, )()(lim00vr（1） 若系统的 NVE 恒定，则NVEAA（微正则系综平均 microcannonical ensemble average ）（2） 若系统的 NVT 恒定，则NVTAA（正则系综平均 cannonical ensemble average ）在实际计算中， 往往还需要要求系统的总动量P=0 或恒定，因整个系统不受外界的作用。（3） 温度与能量均分定理)3(21ckNNkTEN总粒子数， Nc约束条件个数 , 一般情况下， N NcP=0 Nc = 3 （4） 位势切断误差与修正g(r) 对关联函数（ pair correlation function）VN 粒子数密度rr dg)(当原点位置上一个粒子时，在r 附近 dr 内发现有一个粒子的几率。对气体或液体，)()(rr（各向同性）平均每个粒子的能量cccrrrVVdrrrgruNUdrrrgrudrrrgrudrrrgruddrrrgrudrgruNU02202022)()(2)()(2)()(2)()(2)()(21)()(21r切断误差修正为：crcdrrrgruU2)()(26计算机模拟的组织（1） 初始化给定粒子的初始位置、初始速度等；（2） 趋衡从初始转台开始，按运动方程要求的规律，从非平衡态趋向平衡态；（3） 投产计算物理量的统计平均值。 2.2 微正则系综分子动力学方法微正则系综： NVE 恒定，且 P = 0 （分子动力学方法的自然选择）1Verlet 算法运动方程的显示中心差分格式mttttttiiii/)()()(2)(2Frrrttttttiii2/)()()(rrv由可得到如下递推公式2/ )(/211211hrrvmhFrrrninininininini（2-2）此即 Verlet 算法（A2）, 其特点是：（1） 要启动此算法，必须已知两点10,iirr，所以又称为二步法；（2）nir和niv不能同步算出，第n+1 步算出的第 n 步的速度。2Verlet 算法的速度形式2222)(21)()()()(21)()()(tmttttdttdtdttdtttiiiiiiiFvrrrrr2)()()()()()()()()()(tmtttttttmttttttmttttFFvvFvvFvv（向后差分）（向前差分）由此可得如下递推公式：2/ )(2/1121mFFhvvmFhhvrrnininininininini（2-3）此即 Verlet 算法的速度形式 （A3）, 其特点是：（1） 启动此算法，必须已知两点00,iivr；（2）1nir和1niv同步计算。（1niF的计算需要知道1nir）此算法比 A2 算法更稳定。3趋恒阶段的能量调整由于系统初始状态的能量离我们要模拟系统的能量（或温度）有一定差异，这就需要在系统趋于平衡的阶段进行调整。最常用的方法之一就是进行速度标度：11ninivv标度因子：2/12/) 33(irefmvkTN（2-4）Tref为设定的系统的温度值。能量设定值：refrefkTNE)33(21能量参考值：221imvE2EEref2)(21irefvmE要注意的问题：（1） 一般不需要每一步都进行标度，可每隔若干步（比如50 步）调整一次；（2） 当能量达到所设定的值后，停止标度，在以下的模拟时间内能量保持不变。4实例氩原子系统，N=256 LJ势: 6124)(ijijijrrru作用力：814221)(48)()(ijijjiiijijxrrxxxrurF取时间单位：sm122/1210348取长度单位：= 0.3405nm 取质量单位： m=6.63382 10-26 kg(氩原子质量 ) 8142814221482148)(ijijjiijijjiijxrrxxmmrrxxrF使用上述单位，可得到约化单位下的作用力表达式：8*14*11)()(ijijjiijxrrxxrF标度因子：2/12*2/12*2/12*222/1216/)1(16/)1(48/)13(/)33(iirefirefirefvTNvkTNvmmkTNmvkTN约化温度：kTTkTrefref/*Kk8.119取时间步长 h = 2 10-14s = 20fs, 相当于约化时间0.064 约化温度： T* = 2.53 303K, T* = 0.722 86.5K 约化数密度：* = N/L3, N = 256 * = 0.636 L = 7.83 , * = 0.83134 L = 6.75 当 N = 64 时， * = 0.83134 L = 4.25, 所以，取 rc = 2.5 是错的，因为要求rc L/2 = 2.13。 （习题 3.5）计算源程序见 p129, 程序 PL1 2.3 正则系综分子动力学方法正则系综： NVT 恒定，且 P = 0 如何保持温度 T 恒定？1速度标度方法11ninivv2/12/)33(irefmvkTN算法 A4: (1) 规定初始位置0ir和初始速度0iv；(2) 计算2/ )(2/1121mFFhvvmFhhvrrnininininininini；(3) 计算2/12/)33(irefmvkTN；(4) 对速度进行标度：11ninivv，返回（ 2） 。说明：（1） 在前面的微正则系综的模拟中也对速度进行了标度，但其目的是使总能量达到所需要得值，而且不必每步进行速度标度，可隔若干步标度一次，一旦总能量达到所设定的值，即可停止标度；而在正则系综的模拟中，每步都必须进行标度，以始终保持温度恒定。（2） 通过一个广义位势引进能量涨落，连同细致耦合的一种特殊选择，会导致速度标度机制。（见书中的证明，公式3.59 有误）（3） 由于控制温度的反馈环内的时间延迟，温度将在一定程度上涨落。2随机方法（ Anderson 热浴）体系与一强加了温度的热浴相耦合，与热浴的耦合由偶尔作用于随机选择的粒子上的脉冲力表示。在开始模拟前首先要选择系统与热浴的耦合强度，这一耦合强度由随机碰撞的频率决定。模拟步骤：（1）规定初始位置0ir和初始速度0iv；（2）计算2/)(2/1121mFFhvvmFhhvrrnininininininini；（3）粒子 i 在时间间隔 h内发生碰撞的几率为h ，产生随机数，如果 h , 则粒子 i 经历一次随即碰撞，否则，返回（2） ；（4）粒子 i 从温度为 T 的麦克斯韦玻尔兹曼分布中获得一个新的速度，返回（ 2） 。 2.4 等温等压系综分子动力学方法1等压等焓系综N P H 恒定， P = 0 (P:压强， P：总动量 ) VPEHE(PE为外压强，平衡时， PE =P) VPEHE若系统与外界绝热，则VPEE，即0H（等焓过程）要保持压强不变，体积则必须变化，为此，将体积 V 作为一个新的动力学变量，构造一个新的拉氏函数：VPVMUmUTUTVVEiiVV2221)(21),(rrrrL引入标度坐标：LrVriii/3/1或iiLriiiiiiLLLmpvr)(（0L，体积变化缓慢）iiiiimLrmrL2，VMV（共轭动量）VPVMLUmLVVEii22221)(21),(L系统的哈密顿量为：VPVMLUmLVVHEii22221)(21),(（2-5）运动方程为：iiHmLdtd)(2（对粒子）（2-6）VHVMdtd)(（对活塞）（2-7）将（2-5）代入（ 2-6）得iiiiiiiLFrrUUmLLmL22LLLmFiii2VVdtdVdVdVdtdVdtdLL3/23/13/131VVVVLL31313/13/2所以，VVLmFdtdiii3222（2-8）将（2-5）代入（ 2-7）得EEPPPVUVM所以，MPPdtVdE22（2-9）由维里定理?jiijijiijiijijFrmLKPV312132313222FrjiijijiiFrVmLP31312(2-10) 归纳起来：jiijijiiEiiiFrVmLPMPPdtVdVVLmFdtd31313222222这些方程给出一个恒定的平均压强。下面推出递推公式，利用泰勒展开保留到二阶项得：222122212121dtVdhdtdVhVVdtdhdtdhnnnnnnnn将（2-8）和（ 2-9）式代入上式得：)(21312121221EnnnnnnnnnnnnPPhMVhVVVVhmLFhh（2-11）定义偏速度：)(21312111EnnnnnnnnnnPPhMVVVVhmLFh递推公式可写为：1111nnnnnnVhVVh为了计算第 n +1 步的压强，需要知道第n +1 的速度，为了解决这个问题，利用偏速度近似有：)(2131211211111211211EnnnnnnnnnnPPhMVhVhVVhmLFhhh（没有严格的证明）压强的递推公式为：jinijnijinninnFrVmLP111211131)(31（2-12）计算时，用偏速度1n代替1n。于是，可把上述算法表述为如下形式：算法 A5（NPH 系综）（1） 规定初始位置和初始速度：0000,iiiirr；（2） 规定一个同所要求的密度不矛盾的0V；（3） 规定0V（比如，00V） ；（4） 由递推公式计算11,nnV；（5） 计算偏速度11,nnV；（6） 计算111,nnijnFrF；（7） 利用偏速度计算1nP；（8） 计算1nV；（9） 计算1n。将上述算法同速度标度结合起来，就可得到等温等压系综算法，即算法 A6（NPT 系综）（1）- （9）和算法 A5 相同；（10）速度标度11nn。