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

硕士研究生课程 硕士研究生课程 《物理问的计算机模拟方法》讲义 适用专业: 凝聚态物理、材料物理与化学、理论物理、光学工程 学 时：30—40学时 参考教材： 1．[德]D.W.Heermann 著，秦克诚译，理论物理中的计算机模拟方法，北京大学出版社，1996。 2．[荷] Frenkel & Smit 著，汪文川 等译，分子模拟—从算法到应用，化学工业出版社，2002。 3．M.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 序言 1.2 热力学系统物理量的统计平均 1.3 分子动力学方法模拟的基本思想 1.4 蒙特卡罗方法模拟的基本思想 1.5 元胞自动机模拟的基本思想 1.5.1 简要的发展历程 1.5.2 简单元胞自动机：奇偶规则 1.5.3 元胞自动机的一般定义 第2章​ 确定性模拟方法—分子动力学方法（MD） 2.1 分子动力学方法 2.2 微正则系综分子动力学方法 2.3 正则系综分子动力学方法 2.4 等温等压系综分子动力学方法 第3章​ 随机性模拟方法—蒙特卡罗方法（MC） 3.1 预备知识 3.2 布朗动力学（BD） 3.3 蒙特卡罗方法 3.4 微正则系综蒙特卡罗方法 3.5 正则系综蒙特卡罗方法 3.6 等温等压系综蒙特卡罗方法 3.7 巨正则系综蒙特卡罗方法 第4章​ 离散性模拟方法—原胞自动机（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．时间平均 ← 分子动力学（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 , 2Δt, 3Δt, ……, LΔt (t0 = 0) 计算物理量 A(r1(lΔt), r2(lΔt), ……, rN(lΔt)) (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个相等的原胞（s>>1） 注意：由于积分中不含动量，所以我们只需要在位置空间积分，而不需要在相空间中积分。 当系统的代表点落入第i个原胞时，则认为系统处在状态i，因此，s为系统可能的微观状态数目。 于是，积分（10）和（11）可近似表述为 (1-12) (1-13) （2）​ 建立马尔可夫(Maρkoв)过程（链） 将s个状态看作一组随机事件 马尔可夫链：从状态i→j状态j(i→j)的概率pij ，只与i和j 有关。 ，i=1, 2, …, s 若i 经历n步到达j ，其概率表示为 ，存在极限概率 ( j=1, 2, …, s) uj为系统处在状态j的概率。 于是，沿无限长的马尔可夫链，物理量A的平均值可写为 (1-14) 选取 ，则（14）式为A的正则系综平均值。 （3）抽样方法 采用怎样的抽样方法所构成的马尔可夫链能得到上述平均值？ 粒子位置坐标： 粒子编号： =1, 2, …N 坐标的三个方向： = 1, 2, 3 系统状态：i = 1, 2, …, s 给定粒子位置坐标的变化量 (小于系统体积的限度) 给定系统的初态i，随机选定4个随机数，其中三个 （ =1， 2，3），且 –1 1，一个表示粒子编号 =1, 2, …N，由此随机确定粒子位置的变化： （ 确保 ） 若 ，则 运动到新的位置，即系统由状态i过渡到状态j； 若 ，则再选一个随机数 4（0 4 1）， 若 ，则粒子保留在原位置，不发生ij的跃迁； 若 ，则发生i j的跃迁。 由此进行下去，则形成一个马尔可夫链（或过程），此链的长度L(即粒子行走的步数，远大于s)，由所计算的物理量的平均值 (1-15) 不再随链的加长而改变来确定。由此得到的平均值即正则平均值。一般来说，L与N, V, T 有关，比如，N=32~108, L=3000~5000。 归纳起来，计算系统物理量的正则系综平均值的具体步骤如下： 第一步：给定系统的初始状态(粒子的初始位置)ri和每一步的改变量； 第二步：选择四个随机数，其中一个代表粒子的编号i ( 1 i N )；另外三个表示粒子空间坐标的改变x, y, z ( - , =1, 2, 3); 第三步：计算粒子i的新位置 第四步：计算粒子在新旧两个位置系统的能量之差 第五步：由 的大小判断粒子i是否从ri运动到 : 若 ，则ri ； 若 ，则再选一个随机数R(0 < R < 1), 如果 ，则ri ； 如果 ，则ri 不变，返回第二步。 第六步：计算 第七步：重复上述各个步骤，直到完成L步为止，最后利用公式（15）计算A的平均值。 3．​ 粒子间相互作用势模型的选取 最简单的两种模型： （1）硬球模型 ( 为硬球的直径) （2）L—J 势 § 1.5 元胞自动机(CA)模拟的基本思想 元胞自动机：时间和空间都离散、物理参量只取有限数值集的物理系统的理想化模型 cellular automata 或 cellular automaton — CA 1.5.1 简要的发展历程 1．自繁殖系统 20世纪40年代，Von Neumann, 构造能解决非常复杂问题的计算机，设想模仿人脑的行为——寻求与生物过程无关的情况下自繁殖机理的逻辑抽象。 根据S. Ulam的建议，Von Neumann 在由元胞构成的完全离域的框架下处理这个问题，构造了一个完全离散的动力学系统——元胞自动机。 第一个自复制元胞自动机——由二维方形网络组成，有数千个基本元胞构成的自繁殖结构。 （1）一般机器只能构造比自己简单的客体，而采用自复制元胞机，可获得一种能产生新的、具有同样复杂性和功能的“机器”； （2）Von Neumann 的元胞自动机规则具有所谓通用计算的性质，这意味着，存在一种元胞自动机的初始构型，该元胞自动机能产生任何计算机算法的解。通用计算的性质指：用元胞自动机演化规则能够模拟任何计算机流程（逻辑选择器开关）。 2．生命游戏机 1970年，数学家John Conway 生命游戏机的概念，寻找能导致复杂行为的简单规则。 设想一个类似于棋盘的二维方形网格，每个元胞可能的状态是活（状态1）或死（状态0），其更新规则是：有三个活元胞包围的一个死元胞恢复为活元胞；由两个以下或三个以上活元胞包围的活元胞因孤立或拥挤而死亡。 结果表明，生命游戏机有出乎意料的丰富行为，从原“汤”中显示出来的复杂结构，演变发展成为某些特殊的技艺，例如，可能形成所谓的滑翔机——紧邻元胞的特殊排列，这些元胞具有沿直线弹道穿越空间运动的特性。 生命游戏机也是具有计算通用性的元胞自动机。 3．模拟物理系统 （1）20世纪70年代，Hardy, Pomeau 和 Pazzis 建立了所谓的HPP格子气体模型，用以在质量和动量守恒的情况下在方形网格上模拟粒子的碰撞行为。 （2）1986年，Frisch, Hasslacher 和 Pomeau 提出了著名的FHP模型，这是在六边形网格上模拟二维流体动力学的第一个严格模型——全离散计算机模型替代风洞试验。 HPP和 FHP 通常称之为格子气自动机（LGA—Lattice Gas Automata） （3）Ising自旋动力学模型，20世纪80年代末，Vichniac提出Q2R规则。 （4）格子Boltzmann 方法与多粒子模型 格子Boltzmann 方法或模型（LBM）：网格上定义的物理模型，在这个网格上与每个格位相关联的变量是平均粒子数，或具有一定速度粒子出现的概率。该模型可以用均化或因子分解方法由元胞自动机动力学推导出来，或者自定义，而与特定的实现无关。 格子Boltzmann 模型保持了元胞自动机方法的微观水平解释，但忽略了多体的相关函数，但这种方法已经成为目前模拟物理系统中最有前途的方法之一。 在严格的元胞自动机方法与较灵活的格子Boltzmann 模型之间，有一种目前处于发展中的模型——多粒子模型。这种模型保留量化状态的概念，但接受无限数值集，因此既保证了数值稳定性（与LBM相反），又考虑了多体相关性。在模拟物理系统时，大量的可能状态提供更多的灵活性，并产生小的统计噪声。但多粒子动力学更难设计，且在数值计算上比格子Boltzmann 方法更慢。 1.5.2 简单元胞自动机：奇偶规则 简单元胞自动机的演化规则（20世纪70年代，Edward Fredkin 提出，定义在二维方形网上）： 网格的每一个格位是一个元胞，以其位置r =(i, j) 来标记，其中i和j 为行和列的标号。函数 描述每个元胞在时间t的状态，其值为0或1。从时间 t = 0 的初始条件及网格上给定的构形值 开始，t = 1 时状态按下列步骤求得： （1）对于每个格位r，都计算出其位于东、南、西、北4个最近邻格位 的 值之和。应使系统在i和 j 两个方向循环（如同在环面上），从而确定出所有格位的 计算值。 （2）如果这个和值为偶数，则新状态 为0（白色），否则为1（黑色）。 重复上述步骤，得到t = 2, 3, 4, …的状态。 这个元胞自动机的奇偶规则可表示为： 式中，符号代表“异或”逻辑运算，也即模2和：11=00=0，10=01=1。 反复迭代这个规则时，可得到非常精美的几何图形。 1.5.3 元胞自动机的一般定义 定义： （1）​ 规整的元胞网格覆盖d 维空间的一部分； （2）​ 归属于网格的每个格位r 的一组布尔变量(r, t) = {1(r, t), 2(r, t), …, m(r, t)}给出每个元胞在时间t = 0, 1, 2, …的局部状态； （3）​ 演化规则 R = { R1, R2, …, Rm}按下列方式指定状态(r, t)的时间演化过程： j(r, t+1) = Rj[(r, t), (r + 1, t), (r + 2, t),…, (r + q, t)] 式中r + q 指定从属于元胞r的给定邻居。 邻居： 二维元胞自动机的两种邻居：（1）Von Neumann邻居，有一个中心元胞（要演化的元胞）和四个位于其近邻东西南北方位的元胞组成；（2）Moore邻居，除了前面涉及的最近邻元胞外，还包括次近邻的4个元胞，共九个元胞。 还有一种有用的邻居称为Margolus邻居，将空间划分成22元胞的邻接单元块，这个规则对位于单元块内的位置即左上、右上和左下、右下很敏感。 三种邻居如下图所示。 边界条件： （1）​ 周期性边界条件； （2）​ 固定边界条件； （3）​ 绝热边界条件； （4）​ 映射边界条件。 备注： 确定型元胞自动机：演化规则确定，给定的初始状态将始终演化出同样的式样。 概率型自动元胞机：演化规则包含一定的随机性，给定的初始状态可能演化出不同的式样。 第二章 确定性模拟方法—分子动力学方法（MD） § 2.1 分子动力学方法 用分子动力学方法模拟气体、液体或固体系统 1．分子动力学的描述形式： 哈密顿描述 拉格朗日描述 牛顿运动方程描述 2．划分MD元胞 选取适当大小的 V=L3，太大耗费计算时间，太小不能准确反映系统的性质。 目的：维持系统恒定的密度。 对平衡态系统，液体和气体原胞的形状无关紧要，但对晶体则不然。 3．​ 周期性边界条件 目的：减少表面效应 4．​ 相互作用势、最小影像约定和切断距离 n = (n1, n2, n3 ) 元胞内的粒子间 元胞内的粒子与影像 的相互作用 粒子的相互作用 最小影像约定(minimum image convention): 最小影像 元胞中的一个粒子只与元胞中的另外N-1个粒子中的每一个或其最近邻影像发生相互作用。 切断距离（cutoff distance）：rc L/2 （切断半径） 上述约定相当于用rc 来切断位势，即rij > rc 的相互作用忽略不计，代价是忽略了背景。 5．​ 积分格式 于是，运动方程可写为 (2-1) （1）递推公式 在第一章中求解运动方程时，我们是直接求解的关于粒子位置坐标的二阶微分方程，得到的递推公式需要知道最初的两点位置才能启动计算，但在实际计算中，我们常常是给出最初的位置和速度，于是，我们可通过选取一定的差分格式，有运动方程(16)得到关于位置和速度的递推公式。 或 已知 递推公式的具体形式取决于差分各式的选取。 （2）时间步长 选择时间步长的原则：在保证计算精度的前提下，尽量节省计算时间。 实例：Ar原子系统，L—J势 时间步长取为 (=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 （分子动力学方法的自然选择） 1．Verlet 算法 运动方程的显示中心差分格式 由可得到如下递推公式 （2-2） 此即Verlet 算法（A2）, 其特点是： （1）​ 要启动此算法，必须已知两点 ，所以又称为二步法； （2）​  和 不能同步算出，第n+1步算出的第n步的速度。 2．Verlet 算法的速度形式 由此可得如下递推公式： （2-3） 此即Verlet 算法的速度形式（A3）, 其特点是： （1）​ 启动此算法，必须已知两点 ； （2） 和 同步计算。（ 的计算需要知道 ） 此算法比A2算法更稳定。 3．趋恒阶段的能量调整 由于系统初始状态的能量离我们要模拟系统的能量（或温度）有一定差异，这就需要在系统趋于平衡的阶段进行调整。最常用的方法之一就是进行速度标度： 标度因子： （2-4） Tref 为设定的系统的温度值。 能量设定值： 能量参考值： 要注意的问题： （1）​ 一般不需要每一步都进行标度，可每隔若干步（比如50步）调整一次； （2）​ 当能量达到所设定的值后，停止标度，在以下的模拟时间内能量保持不变。 4．​ 实例 氩原子系统， N=256 L—J势: 作用力： 取时间单位： 取长度单位： = 0.3405nm 取质量单位：m=6.6338210-26 kg(氩原子质量) 使用上述单位，可得到约化单位下的作用力表达式： 标度因子： 约化温度： 取时间步长h = 210-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．速度标度方法 算法A4: (1)​ 初始位置 和初始速度 ； (2)​ 计算 ； (3)​ 计算 ； (4)​ 对速度进行标度： ，返回（2）。 说明： （1）​ 在前面的微正则系综的模拟中也对速度进行了标度，但其目的是使总能量达到所需要得值，而且不必每步进行速度标度，可隔若干步标度一次，一旦总能量达到所设定的值，即可停止标度；而在正则系综的模拟中，每步都必须进行标度，以始终保持温度恒定。 （2）​ 通过一个广义位势引进能量涨落，连同细致耦合的一种特殊选择，会导致速度标度机制。（见书中的证明，公式3.59有误） （3）​ 由于控制温度的反馈环内的时间延迟，温度将在一定程度上涨落。 2．随机方法（Anderson 热浴） 体系与一强加了温度的热浴相耦合，与热浴的耦合由偶尔作用于随机选择的粒子上的脉冲力表示。 在开始模拟前首先要选择系统与热浴的耦合强度，这一耦合强度由随机碰撞的频率 决定。 模拟步骤： （1）规定初始位置 和初始速度 ； （2）计算 ； （3）粒子i在时间间隔h内发生碰撞的几率为h，产生随机数，如果< h, 则粒子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） 定义偏速度： 递推公式可写为： 为了计算第n +1步的压强，需要知道第n +1的速度，为了解决这个问题，利用偏速度近似有： （没有严格的证明） 压强的递推公式为： （2-12） 计算时，用偏速度 代替 。于是，可把上述算法表述为如下形式： 算法A5（NPH系综） （1）​ 规定初始位置和初始速度： ； （2）​ 规定一个同所要求的密度不矛盾的 ； （3）​ 规定 （比如， ）； （4）​ 由递推公式计算 ； （5）​ 计算偏速度 ； （6）​ 计算 ； （7）​ 利用偏速度计算 ； （8）​ 计算 ； （9）​ 计算 。 将上述算法同速度标度结合起来，就可得到等温等压系综算法，即 算法A6（NPT系综） （1）-（9）和算法A5相同； （10）​ 速度标度 。