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

硕士研究生课程《物理问题的计算机模拟方法》讲义适用专业:凝聚态物理、材料物理与化学、理论物理、光学工程学时：30—40学时参考教材：1．[德]D.W.Heermann著，秦克诚译，理论物理中的计算机模拟方法，北京大学出版社，1996。2．[荷]Frenkel&Smit著，汪文川等译，分子模拟—从算法到应用，化学工业出版社，2002。3．M.P.AllenandD.J.Tildesley,ComputerSimulationofLiquids,ClarendonPress,Oxford,1989.4.A.R.Leach,MolecularModelling:PrinciplesandApplications,AddisonWesleyLongman,England,1996.5.[德]D.罗伯著，计算材料学，化学工业出版社，2002。6.[英]B.Chopard&MichelDroz著，物理系统的元胞自动机模拟，祝玉学，赵学龙译，清华大学出版社，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．什么是计算机模拟？SimulationModelling2．为什么要进行计算机模拟？3．常用的计算机模拟方法确定性模拟方法：MD模拟MolecularDynamics随机性模拟方法：MC模拟MonteCarlo离散性模拟方法：CA模拟CellularAutomata§1.2热力学系统物理量的统计平均描述系统的坐标（自由度）：x(t)={x1(t),x2(t),…xN(t)}系统的物理量：A(x(t))1．时间平均dttAttAttt0))((10x←分子动力学（MD）模拟(1-1)2．系综平均xxxxxxdAdHfAZA)()())(()(1←蒙特卡罗(MC)模拟(1-2)))((1)(xxHfZ—分布函数（几率密度函数）(1-3)xxdHfZ))((—配分函数(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)的时间平均获得，即ttdtAttA0)]([1)(X（离散情况：kiiAktA11)(）对于平衡态：)(limtAAt实际模拟时间总是有限的，模拟时间的长短可通过判断时间的增加对平均值的影响来确定，当继续增加时间带来的平均值得变化在允许的误差范围之内时，即可认为模拟足够长了。2．计算步骤运动方程：)(,NiiiiUdtdmdtdrvvr即iNiiUdtdmFrr)(22(1-6)或mdtdii/22Fr(1-7)数值求解：用差分近似表示微分（采用不同的差分格式，可得到不同的算法）。用显示中心差分格式，将（7）式写为222)/()]()(2)([ttttttdtdiiiirrrr(1-8)由（7）和（8）式可得：mttttttiiii/)()()(2)(2Frrr(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的平均值LlNLtltltlLA121))(,),(),(A(1limrrrL的大小由继续增大L而A不变（或变化在误差范围内）来确定。§1.4蒙特卡罗（MC）方法模拟的基本思想1．基本原理以正则系综（T,V,N）为例正则分布：sEseZ1正则配分函数：)(!13NNENddddehNspr系统能量：)(2)(2NiiiNpsUmUEErpr物理量：A(rN)=A(r1,r1,…,rN)系综平均：)(1)(!1)(!1)(!1/)(/2/)(3/332NkTUNNNkTmNkTUNNNNkTENNsNNdeAQdedeAZhNddeAZhNdAhNANiiiNsrrprrprrrrpr（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）建立马尔可夫(Maρkoв)过程（链）将s个状态看作一组随机事件马尔可夫链：从状态i→j状态j(i→j)的概率pij，只与i和j有关。1jijp，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41），若kTUUije/)(4，则粒子保留在原位置，不发生ij的跃迁；若kTUUije/)(4，则发生ij的跃迁。由此进行下去，则形成一个马尔可夫链（或过程），此链的长度L(即粒子行走的步数，远大于s)，由所计算的物理量的平均值LiiLALA11lim(1-15)不再随链的加长而改变来确定。由此得到的平均值即正则平均值。一般来说，L与N,V,T有关，比如，N=32~108,L=3000~5000。归纳起来，计算系统物理量的正则系综平均值的具体步骤如下：第一步：给定系统的初始状态(粒子的初始位置)ri和每一步的改变量；第二步：选择四个随机数，其中一个代表粒子的编号i(1iN)；另外三个表示粒子空间坐标的改变x,y,z(-,=1,2,3);第三步：计算粒子i的新位置ir),,(),,(ziyixiiiiizyxzyxr第四步：计算粒子在新旧两个位置系统的能量之差),...,,...,,(),...,,...,,(2121NiNiUUUrrrrrrrr第五步：由U的大小判断粒子i是否从ri运动到ir:若0U，则riir；若0U，则再选一个随机数R(0R1),如果kTUeR/_，则riir；如果kTUeR/_，则ri不变，返回第二步。第六步：计算),...,,...,,(21NiArrrr第七步：重复上述各个步骤，直到完成L步为止，最后利用公式（15）计算A的平均值。3．粒子间相互作用势模型的选取最简单的两种模型：（1）硬球模型)()(0)(rrru(为硬球的直径)（2）L—J势126)(rBrArujiijNruU)(),...,,(21rrr§1.5元胞自动机(CA)模拟的基本思想元胞自动机：时间和空间都离散、物理参量只取有限数值集的物理系统的理想化模型cellularautomata或cellularautomaton—CA1.5.1简要的发展历程1．自繁殖系统20世纪40年代，VonNeumann,构造能解决非常复杂问题的计算机，设想模仿人脑的行为——寻求与生物过程无关的情况下自繁殖机理的逻辑抽象。根据S.Ulam的建议，VonNeumann在由元胞构成的完全离域的框架下处理这个问题，构造了一个完全离散的动力学系统——元胞自动机。第一个自复制元胞自动机——由二维方形网络组成，有数千个基本元胞构成的自繁殖结构。（1）一般机器只能构造比自己简单的客体，而采用自复制元胞机，可获得一种能产生新的、具有同样复杂性和功能的“机器”；（2）VonNeumann的元胞自动机规则具有所谓通用计算的性质，这意味着，存在一种元胞自动机的初始构型，该元胞自动机能产生任何计算机算法的解。通用计算的性质指：用元胞自动机演化规则能够模拟任何计算机流程（逻辑选择器开关）。2．生命游戏机1970年，数学家JohnConway生命游戏机的概念，寻找能导致复杂行为的简单规则。设想一个类似于棋盘的二维方形网格，每个元胞可能的状态是活（状态1）或死（状态0），其更新规则是：有三个活元胞包围的一个死元胞恢复为活元胞；由两个以下或三个以上活元胞包围的活元胞因孤立或拥挤而死亡。结果表明，生命游戏机有出乎意料的丰富行为，从原“汤”中显示出来的复杂结构，演变发展成为某些特殊的技艺，例如，可能形成所谓的滑翔机——紧邻元胞的特殊排列，这些元胞具有沿直线弹道穿越空间运动的特性。生命游戏机也是具有计算通用性的元胞自动机。3．模拟物理系统（1）20世纪70年代，Hardy,Pomeau和Pazzis建立了所谓的HPP格子气体模型，用以在质量和动量守恒的情况下在方形网格上模拟粒子的碰撞行为。（2）1986年，Frisch,Hasslacher和Pomeau提出了著名的FHP模型，这是在六边形网格上模拟二维流体动力学的第一个严格模型——全离散计算机模型替代风洞试验。HPP和FHP通常称之为格子气自动机（LGA—LatticeGasAutomata）（3）Ising自旋动力学模型，20世纪80年代末，Vichniac提出Q2R规则。（4）格子Boltzmann方法与多粒子模型格子Boltzmann方法或模型（LBM