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

硕士研究生课程

《物理问题的计算机模拟方法》讲义

适用专业: 凝聚态物理、材料物理与化学、理论物理、光学工程

学时：30—40学时

参考教材：

1．[德] 著，秦克诚译，理论物理中的计算机模拟方法，北京大学出版社，1996。

2．[荷] Frenkel & Smit 著，汪文川等译，分子模拟—从算法到应用，化学工业出版社，2002。

3． and Computer Simulation of Liquids, Clarendon Press, Oxford, 1989.

4. Molecular Modelling: Principles and Applications, Addison

Wesley Longman, England, 1996.

5. [德] D.罗伯著，计算材料学，化学工业出版社，2002。

6. [英] B. Chopard & Michel Droz 著，物理系统的元胞自动机模拟，

祝玉学，赵学龙译，清华大学出版社，2003。

目录

第一章计算机模拟方法概论

序言

热力学系统物理量的统计平均

分子动力学方法模拟的基本思想

蒙特卡罗方法模拟的基本思想

元胞自动机模拟的基本思想

1.5.1 简要的发展历程

1.5.2 简单元胞自动机：奇偶规则

1.5.3 元胞自动机的一般定义

第二章确定性模拟方法—分子动力学方法（MD）

分子动力学方法

微正则系综分子动力学方法

正则系综分子动力学方法

等温等压系综分子动力学方法

第三章随机性模拟方法—蒙特卡罗方法（MC）预备知识

布朗动力学（BD）

蒙特卡罗方法

微正则系综蒙特卡罗方法

正则系综蒙特卡罗方法

等温等压系综蒙特卡罗方法

巨正则系综蒙特卡罗方法

第四章离散性模拟方法—原胞自动机（CA）引言

元胞自动机模拟

*元胞自动机模拟的应用

第一章 计算机模拟方法概论

§ 序言

1．什么是计算机模拟

Simulation Modelling 2．为什么要进行计算机模拟 3．常用的计算机模拟方法

确定性模拟方法：MD 模拟 Molecular Dynamics 随机性模拟方法：MC 模拟 Monte Carlo 离散性模拟方法：CA 模拟 Cellular Automata

§ 热力学系统物理量的统计平均

描述系统的坐标（自由度）：x(t)={x 1(t),x 2(t),…x N (t)} 系统的物理量：A (x (t)) 1．时间平均

dt t A t t A t

t t ⎰

-=

))((1

x ← 分子动力学（MD ）模拟 (1-1)

2．系综平均

⎰⎰

Ω

Ω

=>=

x x x

x x d A d H f A Z

A )()( ))(()(1ρ ← 蒙特卡罗(MC)模拟 (1-2)

))((1

)(x x H f Z

=ρ— 分布函数（几率密度函数） (1-3)

⎰Ω=x x d H f Z ))(( — 配分函数

(1-4)

Ω—相空间

H (x )—系统的哈密顿函数 对于处于平衡态的系统，可以证明：

∞>=

对于实际的有限时间内的平均，则有

A A >≈<

实际模拟的系统大小也是有限的：有限的粒子数N 或有限的系统限度L 对统计平均结果有影响。

§ 分子动力学(MD)方法模拟的基本思想

1．基本原理

系统：N 个粒子，体积V ，粒子质量为m

描述一个粒子运动状态的自由度：（r i , p i ） （p i =m v i ） 相空间：6N 维，相空间中的一点的坐标 X N =[r N , (m v N )] r N =(r 1, r 2, …, r N ), v N =(v 1, v 2, …, v N )

粒子间的相互作用势：U(r N

)=U(r 1, r 2, …, r N )=∑

j

i ij u )(r

决定系统相轨迹X N （t ）的运动方程：

⎪⎪

⎪⎪⎩

⎪⎪⎪⎪⎨

⎧==-∇== ( )0(N ..., 2, ,1)

( , 0）加上边界条件（周期性初始条件））

（N N N i i i i X X i U dt d m dt

d r v

v r (1-5)

物理量A 的宏观值，由A(X N ) 的时间平均获得，即

⎰''=t

t d t A t t A 0

)]([1

)(X （离散情况：∑==k i i A k t A 11)(）

对于平衡态：)(lim

t A A t ∞

→= 实际模拟时间总是有限的，模拟时间的长短可通过判断时间的增加对平均值的影响来确定，当继续增加时间带来的平均值得变化在允许的

误差范围之内时，即可认为模拟足够长了。

2．计算步骤

运动方程：

)( , N i i i i U dt

d m dt d r v

v r -∇== 即

i N i i

U dt

d m F r r =-∇=)(22 (1-6)

或

m dt d i i

/2

2F r = (1-7)

数值求解：用差分近似表示微分 （采用不同的差分格式，可得到不同

的算法）。

用显示中心差分格式，将（7）式写为

222)/()]()(2)([t t t t t t dt

d i i i i

∆∆-+-∆+≈r r r r (1-8)

由（7）和（8）式可得：

m t t t t t t i i i i /)()()(2)(2F r r r ∆+∆--=∆+ (1-9) 第一步：由（9）式计算第i 个粒子在t +Δt 时刻的位置坐标。 要启动计算，我们必须要知道最初两点r i (0)和r i (Δt )