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

1、硕士研究生课程物理问题的计算机模拟方法讲义适用专业: 凝聚态物理、材料物理与化学、理论物理、光学工程学 时：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:

2、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 元胞自动机模拟的基本思想 简要的发展历程 简单元胞自动机：奇偶规则 元胞自动机的一般定义第二章 确定

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为什么

4、要进行计算机模拟？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)系统的哈密顿函数对于处于平衡态的系统，可以证明： 对于实际的有限时间内的平均，

5、则有 实际模拟的系统大小也是有限的：有限的粒子数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) 的时间平均获得，即 （离散情况：）对于平衡态： 实际模拟时间

6、总是有限的，模拟时间的长短可通过判断时间的增加对平均值的影响来确定，当继续增加时间带来的平均值得变化在允许的误差范围之内时，即可认为模拟足够长了。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), , r

7、N(lt) (l = 1, 2, , N)第三步：计算物理量A的平均值 L的大小由继续增大L 而<A>不变（或变化在误差范围内）来确定。§ 1.4 蒙特卡罗（MC）方法模拟的基本思想1 基本原理以正则系综（T, V, N）为例正则分布：正则配分函数：系统能量：物理量：A(rN) = A(r1, r1, , rN )系综平均： （1-10） （位形积分） （1-11）用MC方法计算上述多维积分。2 计算步骤（1） 划分原胞 N个粒子3N个空间自由度，3N维空间划分成s个相等的原胞（s>>1）注意：由于积分中不含动量，所以我们只需要在位置空间积分，而不需要在相空间

8、中积分。当系统的代表点落入第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）抽样方法 采用怎样的抽样方法所构成的马

9、尔可夫链能得到上述平均值？ 粒子位置坐标：粒子编号：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），若，则粒子保留在原位置，不发生i®j的跃迁；若，则发生i &#

10、174; j的跃迁。由此进行下去，则形成一个马尔可夫链（或过程），此链的长度L(即粒子行走的步数，远大于s)，由所计算的物理量的平均值 (1-15)不再随链的加长而改变来确定。由此得到的平均值即正则平均值。一般来说，L与N, V, T 有关，比如，N=32108, L=30005000。归纳起来，计算系统物理量的正则系综平均值的具体步骤如下：第一步：给定系统的初始状态(粒子的初始位置)ri和每一步的改变量d；第二步：选择四个随机数，其中一个代表粒子的编号i ( 1£ i £ N )；另外三个表示粒子空间坐标的改变dx, dy, dz ( -d £ da £

11、;d , a=1, 2, 3);第三步：计算粒子i的新位置 第四步：计算粒子在新旧两个位置系统的能量之差 第五步：由的大小判断粒子i是否从ri运动到: 若，则ri®； 若，则再选一个随机数R(0 < R < 1), 如果 ，则ri®； 如果 ，则ri 不变，返回第二步。第六步：计算第七步：重复上述各个步骤，直到完成L步为止，最后利用公式（15）计算A的平均值。3 粒子间相互作用势模型的选取最简单的两种模型： （1）硬球模型 (s 为硬球的直径) （2）LJ 势§ 1.5 元胞自动机(CA)模拟的基本思想元胞自动机：时间和空间都离散、物理参量只取有限数值

12、集的物理系统的理想化模型cellular automata 或 cellular automaton CA 简要的发展历程1自繁殖系统 20世纪40年代，Von Neumann, 构造能解决非常复杂问题的计算机，设想模仿人脑的行为寻求与生物过程无关的情况下自繁殖机理的逻辑抽象。 根据S. Ulam的建议，Von Neumann 在由元胞构成的完全离域的框架下处理这个问题，构造了一个完全离散的动力学系统元胞自动机。第一个自复制元胞自动机由二维方形网络组成，有数千个基本元胞构成的自繁殖结构。（1）一般机器只能构造比自己简单的客体，而采用自复制元胞机，可获得一种能产生新的、具有同样复杂性和功能的“机

13、器”；（2）Von Neumann 的元胞自动机规则具有所谓通用计算的性质，这意味着，存在一种元胞自动机的初始构型，该元胞自动机能产生任何计算机算法的解。通用计算的性质指：用元胞自动机演化规则能够模拟任何计算机流程（逻辑选择器开关）。2生命游戏机1970年，数学家John Conway 生命游戏机的概念，寻找能导致复杂行为的简单规则。设想一个类似于棋盘的二维方形网格，每个元胞可能的状态是活（状态1）或死（状态0），其更新规则是：有三个活元胞包围的一个死元胞恢复为活元胞；由两个以下或三个以上活元胞包围的活元胞因孤立或拥挤而死亡。结果表明，生命游戏机有出乎意料的丰富行为，从原“汤”中显示出来的复杂

14、结构，演变发展成为某些特殊的技艺，例如，可能形成所谓的滑翔机紧邻元胞的特殊排列，这些元胞具有沿直线弹道穿越空间运动的特性。生命游戏机也是具有计算通用性的元胞自动机。3模拟物理系统（1）20世纪70年代，Hardy, Pomeau 和 Pazzis 建立了所谓的HPP格子气体模型，用以在质量和动量守恒的情况下在方形网格上模拟粒子的碰撞行为。（2）1986年，Frisch, Hasslacher 和 Pomeau 提出了著名的FHP模型，这是在六边形网格上模拟二维流体动力学的第一个严格模型全离散计算机模型替代风洞试验。HPP和 FHP 通常称之为格子气自动机（LGALattice Gas Auto

15、mata）（3）Ising自旋动力学模型，20世纪80年代末，Vichniac提出Q2R规则。（4）格子Boltzmann 方法与多粒子模型格子Boltzmann 方法或模型（LBM）：网格上定义的物理模型，在这个网格上与每个格位相关联的变量是平均粒子数，或具有一定速度粒子出现的概率。该模型可以用均化或因子分解方法由元胞自动机动力学推导出来，或者自定义，而与特定的实现无关。格子Boltzmann 模型保持了元胞自动机方法的微观水平解释，但忽略了多体的相关函数，但这种方法已经成为目前模拟物理系统中最有前途的方法之一。在严格的元胞自动机方法与较灵活的格子Boltzmann 模型之间，有一种目前处于

16、发展中的模型多粒子模型。这种模型保留量化状态的概念，但接受无限数值集，因此既保证了数值稳定性（与LBM相反），又考虑了多体相关性。在模拟物理系统时，大量的可能状态提供更多的灵活性，并产生小的统计噪声。但多粒子动力学更难设计，且在数值计算上比格子Boltzmann 方法更慢。 简单元胞自动机：奇偶规则简单元胞自动机的演化规则（20世纪70年代，Edward Fredkin 提出，定义在二维方形网上）：网格的每一个格位是一个元胞，以其位置r =(i, j) 来标记，其中i和j 为行和列的标号。函数描述每个元胞在时间t的状态，其值为0或1。从时间 t = 0 的初始条件及网格上给定的构形值开始，t

17、= 1 时状态按下列步骤求得：（1）对于每个格位r，都计算出其位于东、南、西、北4个最近邻格位的值之和。应使系统在i和 j 两个方向循环（如同在环面上），从而确定出所有格位的计算值。（2）如果这个和值为偶数，则新状态为0（白色），否则为1（黑色）。重复上述步骤，得到t = 2, 3, 4, 的状态。这个元胞自动机的奇偶规则可表示为：式中，符号Å代表“异或”逻辑运算，也即模2和：1Å1=0Å0=0，1Å0=0Å1=1。 反复迭代这个规则时，可得到非常精美的几何图形。 元胞自动机的一般定义定义：（1） 规整的元胞网格覆盖d 维空间的一部分；（2）

18、归属于网格的每个格位r 的一组布尔变量F(r, t) = F1(r, t), F2(r, t), , Fm(r, t)给出每个元胞在时间t = 0, 1, 2, 的局部状态；（3） 演化规则 R = R1, R2, , Rm按下列方式指定状态F(r, t)的时间演化过程：Fj(r, t+1) = RjF(r, t), F(r + d1, t), F(r + d2, t), F(r + dq, t)式中r + dq 指定从属于元胞r的给定邻居。邻居：二维元胞自动机的两种邻居：（1）Von Neumann邻居，有一个中心元胞（要演化的元胞）和四个位于其近邻东西南北方位的元胞组成；（2）Moore邻

19、居，除了前面涉及的最近邻元胞外，还包括次近邻的4个元胞，共九个元胞。还有一种有用的邻居称为Margolus邻居，将空间划分成2´2元胞的邻接单元块，这个规则对位于单元块内的位置即左上、右上和左下、右下很敏感。三种邻居如下图所示。边界条件：（1） 周期性边界条件；（2） 固定边界条件；（3） 绝热边界条件；（4） 映射边界条件。备注： 确定型元胞自动机：演化规则确定，给定的初始状态将始终演化出同样的式样。 概率型自动元胞机：演化规则包含一定的随机性，给定的初始状态可能演化出不同的式样。第二章 确定性模拟方法分子动力学方法（MD）§ 2.1 分子动力学方法 用分子动力学方法模拟

20、气体、液体或固体系统1分子动力学的描述形式： 哈密顿描述 拉格朗日描述 牛顿运动方程描述2划分MD元胞选取适当大小的 V=L3，太大耗费计算时间，太小不能准确反映系统的性质。目的：维持系统恒定的密度。对平衡态系统，液体和气体原胞的形状无关紧要，但对晶体则不然。3 周期性边界条件目的：减少表面效应4 相互作用势、最小影像约定和切断距离 n = (n1, n2, n3 ) ß ß 元胞内的粒子间 元胞内的粒子与影像的相互作用 粒子的相互作用最小影像约定(minimum image convention): 最小影像 元胞中的一个粒子只与元胞中的另外N-1个粒子中的每一个或其最近

21、邻影像发生相互作用。切断距离（cutoff distance）：rc £ L/2 （切断半径）上述约定相当于用rc 来切断位势，即rij > rc 的相互作用忽略不计，代价是忽略了背景。5 积分格式 于是，运动方程可写为 (2-1) （1）递推公式在第一章中求解运动方程时，我们是直接求解的关于粒子位置坐标的二阶微分方程，得到的递推公式需要知道最初的两点位置才能启动计算，但在实际计算中，我们常常是给出最初的位置和速度，于是，我们可通过选取一定的差分格式，有运动方程(16)得到关于位置和速度的递推公式。 或 已知递推公式的具体形式取决于差分各式的选取。（2）时间步长选择时间步长的原

22、则：在保证计算精度的前提下，尽量节省计算时间。实例：Ar原子系统，LJ势 时间步长取为 (=100ps=10fs)（3）约束条件保证能量、动量或角动量守恒。（4）减小计算误差的技巧数值计算不可避免有误差，与误差有关的因素主要有： 差分格式 时间步长 切断半径 最小影像约定 等等6. 计算热力学量 （1） 若系统的NVE恒定，则 （微正则系综平均microcannonical ensemble average）（2） 若系统的NVT恒定，则 （正则系综平均cannonical ensemble average）在实际计算中，往往还需要要求系统的总动量P=0或恒定，因整个系统不受外界的作用。（3）

23、 温度与能量均分定理 N总粒子数，Nc约束条件个数, 一般情况下，N>> Nc P=0 Þ Nc = 3（4） 位势切断误差与修正g(r) 对关联函数（pair correlation function） 粒子数密度当原点位置上一个粒子时，在r 附近dr内发现有一个粒子的几率。对气体或液体， （各向同性）平均每个粒子的能量切断误差修正为：6 计算机模拟的组织（1） 初始化给定粒子的初始位置、初始速度等；（2） 趋衡从初始转台开始，按运动方程要求的规律，从非平衡态趋向平衡态；（3） 投产计算物理量的统计平均值。§ 2.2 微正则系综分子动力学方法微正则系综：NVE

24、恒定，且P = 0 （分子动力学方法的自然选择）1Verlet 算法运动方程的显示中心差分格式由可得到如下递推公式 （2-2）此即Verlet 算法（A2）, 其特点是：（1） 要启动此算法，必须已知两点，所以又称为二步法；（2） 和不能同步算出，第n+1步算出的第n步的速度。 2Verlet 算法的速度形式 由此可得如下递推公式： （2-3）此即Verlet 算法的速度形式（A3）, 其特点是：（1） 启动此算法，必须已知两点；（2）和同步计算。（的计算需要知道）此算法比A2算法更稳定。3趋恒阶段的能量调整由于系统初始状态的能量离我们要模拟系统的能量（或温度）有一定差异，这就需要在系统趋于平

25、衡的阶段进行调整。最常用的方法之一就是进行速度标度： 标度因子： （2-4）Tref 为设定的系统的温度值。 能量设定值： 能量参考值： 要注意的问题：（1） 一般不需要每一步都进行标度，可每隔若干步（比如50步）调整一次；（2） 当能量达到所设定的值后，停止标度，在以下的模拟时间内能量保持不变。4 实例氩原子系统， N=256LJ势:作用力： 取时间单位：取长度单位：s = 0.3405nm取质量单位：m=6.63382´10-26 kg(氩原子质量)使用上述单位，可得到约化单位下的作用力表达式：标度因子： 约化温度： 取时间步长h = 2´10-14s = 20fs,

26、相当于约化时间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） 由于控制温度的