分子动力学入门

模拟的数学方法
取 N 102~ 106 ， 前比值为0.2~0.01。取前值，模拟 粗糙；取后值，模拟计算量太大 处理方法：使用周期性边界条件 –周期性边界条件 a




两个不同粒子在x或y方向上的最大分离距离为a/2
a ? 2
a/2
a
–最小像约定：两粒子分离距离>最大分离距离， 相互作用力可以忽略，而加入其中像粒子之一相 互作用力来考虑
arccos(2(random 0.5))
sign (random 0.5) arccos [2(random 0.5)]
0 180 ; 180 180
分量 ：
z
vx v sin cos


y
vy v sin sin
vz v cos
动力学模拟：预言反应过程、验证理论猜想、理解实验观测现象。
/nano/IWGN.Research.Directions/
Top 500 Supercomputers in the world
Four orders of magnitude in 15 years A “small” PC cluster today
• 假设
–分子为球，惰性，分子间的作用只取决于分子 间的距离
• 分子动力学 (Molecular Dynamics，MD)
MD的应用
• 领域：物理、化学、生物、材料等
• MD方法能实时将分子的动态行为显示到计算机屏 幕上, 便于直观了解体系在一定条件下的演变过程
• MD含温度与时间, 因此还可得到如材料的玻璃化转 变温度、热容、晶体结晶过程、输送过程、膨胀过 程、动态弛豫(relax)以及体系在外场作用下的变化过 程等
分子间势能及相互作用
• N个粒子系统的总势能
V V (r12 ) V (r13 ) V (r1N ) V (r23 ) V (r2 N ) V (rN -1,N )
i j 1
V (r )
ij
N
V
V
V d
r
V (r ) 0
12 6
分子间势能及相互作用
• 一些气体的参数
Neon (nm) /kB(K) 0.275 36 Argon 0.3405 119.8 Krypon 0.360 171 Xenon Nitrogen 0.410 221 0.370 95
kB=1.38x1023(J/K): Boltzmann常数
2 xn vn t 1 a t 2 n
模拟的数学方法
• Leap-frog算法：
vn 1/2 vn 1/2 tan xn 1 xn tvn 1/2 vn vn 1/2
t 2Βιβλιοθήκη anx的截断误差为 O(t 4 ) ，v的截断误差 O(t 2 )
模拟的数学方法
V(r), F(r)
V(r) F(r)
排斥力
吸引力
能量尺度； 长度尺度
为方便，时常归一化：
1 12 1 6 V (r ) 4 r r
r
V (r )

记 V / V ； r / r
1 1 4 r / r /
多体系统—分子动力学
Multi-body System Molecular Dynamics
热力学统计物理数值实验
Science Research
Theory
Experiment
Computation
计算机模拟已经与理论与实验并列，成为三种基本的科学研究手段之一
Scientific Computations
F 0.02, Twall 1.4
CASE(1) - Couette Flow
z
time
CASE(1) - Couette Flow
z
CASE(2) - Contact Angle Simulation
Mass: m[1]=1, m[2]=8, m[3]=0.8 L=25.05, W=6.56 H=10.29 T=1.2
x
模拟
• 微观量
–温度 根据统计热力学，平衡态下经典系统的能量中的 每一个二次项具有平均值kBT/2,即
d 1 Nk BT mi vi2 2 2 i
空间维数 粒子个数 ：取时间平均
注意：上式在系统质心速度为0时适用
模拟
问题：如何给定系统的初始条件，得到所需要的平 衡态温度Teq？ 解决方法之一：速度标定法 任给初始条件，模拟到平衡，得到系统平衡态温 度T。一般TTeq。令
f Teq / T
用速度
f
1/ d
vi vi
再模拟直到平衡，若所得温度仍不等于Teq，再进 行上述过程
给定初始条件：xi，vi
计算到平衡态
Y
计算结束
|T - Teq|≤ ε
N
f = Teq/ T vi = vi f 1/d
其它方法：Gaussian热浴法（约束温度调节方法） 其基本原理在运动方程中加入“摩擦力”项， 并将其与粒子速度联系起来。
计算量随体系大小急剧增长
引 言
• 物质基本构成—分子、原子
–在分子、原子这个微观水平上来考察物质：多 体世界 –查清楚微观世界，宏观就清楚了
• 从微观考虑问题的现实可行性 • 从微观考虑问题的必要性
–物性的观测性参数：热传导、温度、压力、粘 性、... …
• 微观处理的前提
–已知微观粒子间的相互作用
模拟的数学方法
考虑的粒子总数不变
• 初始条件
条 件 一 ： 规 则 给 法 条 件 二 ： 随 机 给 法
–随机初始条件给法之一 要求 | v | Vmax 大小：
v Vmax (2 random1)
模拟的数学方法
模拟的数学方法
random：随机数产生函数，产生(0，1)之间的随 机数。 方向(按球坐标给法)：
• 边界条件
–模拟能力限制，不能模拟大量分子，只能模拟有 限空间中的有限个分子：有限空间边界 –固体（刚性）边界条件
• 不仅仅有分子间的相互作 用，还引入了壁面的作用 • 分子量大时，壁面作用可 忽略不计
总分子数 N a3
和壁面作用分子数 壁面积
和壁面作用分子数 壁面积 6a 2 1 1 3 3 总分子数 体积 a a N
1 1 0.2 1 1 1
6 ij 12 ij VLJ (r ) ij ij r r
ij 1
0.2 1
i, j=1,2,3, 1 – red fluid, 2 -- wall, 3 – green fluid
1 T dNkB
m v
i
2 i i
模拟
–热容 定义热容
E：系统总能
E Cv T V 计算系统在温度T和T+T时的总能ET、ET +T，
ET T ET Cv T
模拟
–压强
• 对壁面的压强 刚 性 壁
t t
时刻，速度为 vi
t 时刻，速度为 vi
fi mai mvi
平衡态时，系统温度不变，因此dEk/dt=0
v a
i
i i
0
fv m v
i i
i i 2 i
宏观性质的统计
• 系统的势能
• 系统的动能
Ep
1i j N

V (rij )
• 系统的总能 E = Ep+Ek • 系统的温度
p2 Ek i 2mi
• 空间尺度：电子 结构 • 时间尺度：动力 学
methods
Time
properties
systems
/nano/IWGN.Research.Directions/
预言材料性质、验证理论猜想、理解实验观测现象。
/nano/IWGN.Research.Directions/
• 相互作用
–标量形式：
F V (r ) m
f V (r ) m r
–直角坐标：
fy fx V (r ) x V (r ) y , m r r m r r 至此，各粒子间相互作用已知，可进行模拟了
分子间势能及相互作用
模拟的数学方法
• Euler法和Euler-Cromer方法？
1 p dA i vi vi mi t
t时间里作用在单位面积壁上的压力
–粒子速度分布
N (v )
v
v
v v v 2v
选速度间隔v，模拟nt个时间步，记录在每个速 度间隔中的粒子数，最后归一化。
模拟
模拟
• 气、液状态方程
维里定理(Virial Theorem)
MD的基本原理
• 用牛顿经典力学计算许多分子在相空间中 的轨迹
–求解系统中的分子或原子间作用势能和系统外 加约束共同作用的分子或原子的牛顿方程。 –模拟系统随时间推进的微观过程。 –通过统计方法得到系统的平衡参数或输运性质 –计算程序较为复杂，占用较多内存
MD的主要步骤
• 选取要研究的系统及其边界，选取系统内 粒子间的作用势能模型 • 设定系统中粒子的初始位置和初始动量 • 建立模拟算法，计算粒子间作用力及各粒 子的速度和位置 • 当体系达到平衡后，依据相关的统计公式， 获得各宏观参数和输运性质
–不能用：不能保持总能量守恒 –Verlet算法：速度形式
xn 1 vn 1 vn 1 a an t 2 n 1 vn 1 vn 1 2 an t 2 1 v v 1 n 1 2 an t n 2
p k BT

1 b
a 2
气体密度
确定系数a和b
CASE(1) - Couette Flow
• Size of domain is:12.51x7.22x16.71
E 4 ij ( ij
6 ij 12 ij
r
6

r
6
12
)
i, j represent different species
0.5 μm
Fig. 2. The effect of converging geometry obtained by MD simulation of one million particles in the microscale. Dzwinel, W., Alda, W., Pogoda, M., and Yuen, D.A., 2000, Turbulent mixing in the microscale: a 2D molecular dynamics simulation, Physica D, Vol. 137, pp. 157-171.
1 pV Nk BT d
压强
ri Fi
粒子i所受到 的其它粒子 的合相互作 粒子i的位置矢量 用力 在温度的模拟 基础上再模拟 此项
i
体积
温度的模 拟 可得此项
模拟
例：用此可确定高密度气体和液体状态方程(van der Waals方程)中的系数 理想气体状态方程在高密度情况下不可用
ij ij dE F 4 ij (6 ij 7 12 12 ) dr r r
12
If δ is less than 0, then the two species are immiscible.
ll ww wl 1, ll 1, ww 1.5, wl 0.6,
V
CASE(3)- Rayleigh-Taylor instability generation(重力场)
Fig. 1. The snapshots of MD simulations of the R–T instability for two particle systems: (A) closed; (B) open (the heavy fluid is coloured in light grey while the lighter one is dark grey, the black part of the figure is empty, i.e., it does not contain particles).
d
rd rd
r
r
d V (r ) ( ) r
rd V (r ) d ( ) rd r
-
刚球模型
斥力力心点模型
Southerland模型
分子间势能及相互作用
• Lennard-Jones势能
12 6 V (r ) 4 r r