物理问题的计算机模拟方法(1)—分子动力学
- 格式：doc
- 大小：401.50 KB
- 文档页数：22
硕士研究生课程
《物理问题的计算机模拟方法》讲义
适用专业: 凝聚态物理、材料物理与化学、理论物理、光学工程
学 时：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.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))((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)([ttttttdtdiiiirrrr (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(1limrrr
L的大小由继续增大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个相等的原胞（s>>1）
注意：由于积分中不含动量，所以我们只需要在位置空间积分，而不需要在相空间中积分。
当系统的代表点落入第i个原胞时，则认为系统处在状态i，因此，s为系统可能的微观状态数目。
于是，积分（10）和（11）可近似表述为
ikTUiNieAQA/1 (1-12)
ikTUNieQ/ (1-13)
（2） 建立马尔可夫(Maρkoв)过程（链）
将s个状态看作一组随机事件
马尔可夫链：从状态i→j状态j(i→j)的概率pij ，只与i和j 有关。
1jijp，i=1, 2, …, s
若i 经历n步到达j ，其概率表示为)(nijp，存在极限概率
jnijnup)(lim ( j=1, 2, …, s)
jjjuu1 ,0
uj为系统处在状态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，则粒子保留在原位置，不发生ij的跃迁；
若kTUUije/)(4，则发生i j的跃迁。
由此进行下去，则形成一个马尔可夫链（或过程），此链的长度L(即粒子行走的步数，远大于s)，由所计算的物理量的平均值
LiiLALA11lim (1-15)
不再随链的加长而改变来确定。由此得到的平均值即正则平均值。一般来说，L与N, V,
T 有关，比如，N=32~108, L=3000~5000。
归纳起来，计算系统物理量的正则系综平均值的具体步骤如下：
