文章目录
1.分子模拟技术
分子模拟的分类
- Monte Carlo 模拟:1949年,利用随机行走产生微观态。通过系综平均,得到体系的平衡性质。
- 分子动力学模拟:1957年,利用运动轨迹产生微观态。通过时间平均,得到体系的动力学性质。
- 混合方法模拟:QM/MM方法,QM+MD方法,MC+MD方法等。
2.分子动力学的基本概念
分子动力学
- 初始化条件:势函数,相互作用,系综,外界条件,初始构型
- 动力学方程:Newton力学方程
- 一次信息:原子的坐标,三维结构,原子的速度和加速度,原子运动
- 二次信息:热力学性质,动力学性质,光学性质
初始化条件——模拟盒子
- 边界条件:
在小体系中,边界效应一般不能忽略。
如果研究的体系受环境影响较小,可采用自然边界条件;如果研究的体系会受到环境的影响,则最好采用周期边界条件。 - 周期边界条件:
中心盒子在x,y,z方向无限扩展,以消除人为边界效应。
当某个粒子运动出模拟盒子的某一边界时,一个映像粒子从另一对立边界以相同的状态进入到盒子中。
μ = μ − ⌊ μ L μ ⌋ ( μ = x , y , z ) ⌊ ⌋ 为下取整 \mu = \mu -\left \lfloor \frac{\mu }{L_{\mu } } \right \rfloor\quad\left ( \mu =x,y,z \right ) \quad\left \lfloor \right \rfloor 为下取整 μ=μ−⌊Lμμ⌋(μ=x,y,z)⌊⌋为下取整
利用周期边界条件可以基于少量粒子模拟的同时考虑周围环境分子多研究体系的影响。 - 势函数(力场):
描述分子内相互作用和分子间相互作用的数学模型。
分子内相互作用:
(1)键相互作用,例如bond_style命令
U b o n d ( r ) = k r ( r − r 0 ) 2 U_{bond}\left ( r \right )=k_{r} \left ( r-r_{0} \right )^{2} Ubond(r)=kr(r−r0)2
(2)角相互作用,例如angle_style命令
U a n g l e ( θ ) = k θ ( θ − θ 0 ) 2 U_{angle}\left ( \theta \right )=k_{\theta } \left ( \theta -\theta _{0} \right )^{2} Uangle(θ)=kθ(θ−θ0)2
(3)二面角相互作用,例如dihedral_style命令
U d i h e d r a l ( φ ) = k φ [ 1 ± c o s ( n φ ) ] U_{dihedral}\left ( \varphi \right )=k_{\varphi}\left [ 1\pm cos\left ( n\varphi \right ) \right ] Udihedral(φ)=kφ[1±cos(nφ)]
(4)离面角相互作用,例如improper_style命令
U i m p r o p e r ( χ ) = k χ ( χ − χ 0 ) 2 U_{improper}\left (\chi \right )=k_{\chi }\left ( \chi -\chi _{0} \right )^{2} Uimproper(χ)=kχ(χ−χ0)2
分子间相互作用
(1)van der Waals作用,例如pair_style lj命令
U ( r i j ) = 4 ε [ ( σ r i j ) 12 − ( σ r i j ) 6 ] U\left (r_{ij}\right)=4\varepsilon \left [ \left ( \frac{\sigma }{r_{ij} } \right ) ^{12} -\left ( \frac{\sigma }{r_{ij} } \right )^{6} \right ] U(rij)=4ε[(rijσ)12−(rijσ)6]
(2)静电相互作用,例如pair_style coul命令
U ( r i j ) = k q i q j ϵ r i j U\left (r_{ij}\right)=k\frac{q_{i}q_{j} }{\epsilon r_{ij} } U(rij)=kϵrijqiqj
(3)三体相互作用,例如pair_style tersoff命令
U i = 1 2 ∑ j ≠ i f C ( r i j ) [ A e − λ 1 r i j − B ( 1 + β n ζ i j n ) − 1 2 n e − λ 2 r i j ] U_{i}=\frac{1}{2} \sum_{j\neq i}^{}f_{C} \left ( r_{ij} \right ) \left [ Ae^{-\lambda _{1}r_{ij} }-B\left ( 1+ \beta ^{n}\zeta _{ij}^{n} \right ) ^{-\frac{1}{2n} } e^{-\lambda _{2}r_{ij} } \right ] Ui=21j=i∑fC(rij)[Ae−λ1rij−B(1+βnζijn)−2n1e−λ2rij]
其中
ζ i j = ∑ k ≠ i , j f C ( r i k ) g ( θ i j k ) e λ 3 m ( r i j − r i k ) m \zeta _{ij}=\sum_{k\neq i,j}^{} f_{C} \left ( r_{ik} \right ) g\left ( \theta _{ijk} \right ) e^{\lambda _{3}^{m}\left ( r_{ij}- r_{ik} \right)^{m} } ζij=k=i,j∑fC(rik)g(θijk)eλ3m(rij−rik)m
(4)多体相互作用,例如pair_style eam命令
U i = F α ( ∑ j ≠ i ρ β ( r i j ) ) + 1 2 ∑ j ≠ i ϕ α β ( r i j ) U_{i}=F_{\alpha } \left ( \sum_{j\neq i}^{} \rho _{\beta }\left ( r_{ij} \right ) \right )+\frac{1}{2}\sum_{j\neq i}^{} \phi _{\alpha \beta }\left ( r_{ij} \right ) Ui=Fα⎝ ⎛j=i∑ρβ(rij)⎠ ⎞+21j=i∑ϕαβ(rij)
初始化条件——短程势函数的处理
当分子间相互作用使积分
∫
r
c
u
t
∞
U
(
r
)
4
π
r
2
d
r
\int_{r_{cut}}^{\infty} U\left ( r \right ) 4\pi r^{2}dr
∫rcut∞U(r)4πr2dr收敛,则为短程相互作用,否则为长程相互作用。
-
最小映像原理:
中心盒子中的一个粒子只与此盒子中的其他N-1个粒子,或它们的最近映像发生相互作用。
μ i j = μ i j − [ μ i j L μ ] ( μ = x , y , z ) [ ] 为就近取整 \mu_{ij} = \mu_{ij} -\left [ \frac{\mu_{ij} }{L_{\mu } } \right ]\quad\left ( \mu =x,y,z \right ) \quad\left [\quad \right ] 为就近取整 μij=μij−[Lμμij](μ=x,y,z)[]为就近取整 -
势能截断:
不可能将粒子与它们的全部映像间的相互作用都进行计算,势函数可以在一定距离r_{cut}处进行截断,一般势函数的截断半径不大于模拟盒子的最短边长的一半。
(1)简单截断
势能在截断处不连续(分子跨过截断半径时总能量不守恒),力在截断处为无穷大(MD过程不稳定)
U ~ ( r ) = { U ( r ) ( r ≤ r c u t ) 0 ( r > r c u t ) \tilde{U} \left( r\right) =\left\{\begin{matrix}U\left ( r \right ) \quad\left ( r\le r_{cut} \right ) \\0 \quad \quad \quad \left ( r> r_{cut} \right ) \end{matrix}\right. U~(r)={U(r)(r≤rcut)0(r>rcut)
(2)位移截断
势能在截断处连续,力在截断处为有限值但不连续。
U ~ ( r ) = { U ( r ) − U ( r c u t ) ( r ≤ r c u t ) 0 ( r > r c u t ) \tilde{U} \left( r\right) =\left\{\begin{matrix}U\left ( r \right )-U\left ( r_{cut} \right ) \quad\left ( r\le r_{cut} \right ) \\0 \quad \quad \quad\quad\quad\quad\quad \left ( r> r_{cut} \right ) \end{matrix}\right. U~(r)={U(r)−U(rcut)(r≤rcut)0(r>rcut)
(3)位移—力截断
势能在截断处连续,力在截断处连续
U ~ ( r ) = { U ( r ) − U ( r c u t ) − d U ( r ) d r ∣ r c u t ( r − r c u t ) ( r ≤ r c u t ) 0 ( r > r c u t ) \tilde{U} \left( r\right) =\left\{\begin{matrix}U\left ( r \right )-U\left ( r_{cut} \right ) -\frac{dU\left ( r \right ) }{dr}\mid _{r_{cut} } \left ( r-r_{cut} \right ) \quad\left ( r\le r_{cut} \right ) \\0 \quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \left ( r> r_{cut} \right ) \end{matrix}\right. U~(r)={U(r)−U(rcut)−drdU(r)∣rcut(r−rcut)(r≤rcut)0(r>rcut)
初始化条件——长程势函数的处理
- Ewald加和:
长程相互作用(如静电相互作用)需要将粒子与它们的全部映像间的相互作用都进行计算。
U ( r ) = 1 2 ∑ n = − ∞ + ∞ ∑ i , j q i q j 4 π ϵ 0 ( r i j → + n L → ) 其中盒子的边长为 L , n = 0 为中心盒子,否则为映像盒子 U\left ( r \right ) =\frac{1}{2} \sum_{n=-\infty }^{+\infty }\sum_{i,j}^{} \frac{q_{i}q_{j}}{4\pi \epsilon _{0}\left ( \overrightarrow{r_{ij}}+n\overrightarrow{L} \right ) } 其中盒子的边长为L,n=0为中心盒子,否则为映像盒子 U(r)=21n=−∞∑+∞i,j∑4πϵ0(rij+nL)qiqj其中盒子的边长为L,n=0为中心盒子,否则为映像盒子
直接计算静电相互作用收敛慢,故采用Ewald加和算法:在每个点电荷周围加一个各向同性、数值相同、符号相反的Gauss型电荷(屏蔽电荷) ρ i ( r ) = q i α 3 π 3 2 e − α 2 r 2 \rho _{i}\left ( r \right )=\frac{q_{i}\alpha ^{3}}{\pi ^{\frac{3}{2} }} e^{-\alpha ^{2}r^{2}} ρi(r)=π23qiα3e−α2r2
将静电相互作用分为一个短程力和一个长程力。短程力在实空间做截断,长程力在倒易空间做截断。- 实空间部分
U 1 ( r ) = 1 2 ∑ n = − ∞ + ∞ ∑ i , j q i q j 4 π ϵ 0 e r f c ( α ∣ r i j → + n L → ∣ ) ∣ r i j → + n L → ∣ α 越大收敛越快 U_{1}\left ( r \right ) =\frac{1}{2} \sum_{n=-\infty }^{+\infty}\sum_{i,j}\frac{q_{i}q_{j}}{4\pi \epsilon _{0}}\frac{erfc\left ( \alpha\left | \overrightarrow{r_{ij}}+n\overrightarrow{L} \right | \right ) }{\left | \overrightarrow{r_{ij}}+n\overrightarrow{L} \right | }\quad\quad\alpha 越大收敛越快 U1(r)=21n=−∞∑+∞i,j∑4πϵ0qiqj∣ ∣rij+nL∣ ∣erfc(α∣ ∣rij+nL∣ ∣)α越大收敛越快 - 倒易空间部分
U 2 ( r ) = 1 2 ∑ k → ≠ 0 1 ϵ 0 L 3 1 k 2 e − k 2 4 α 2 ∣ ∑ j q j e i k → r j → ∣ 2 α 越小收敛越快, k = 2 π n L 2 U_{2}\left(r\right)=\frac{1}{2}\sum_{\overrightarrow{k}\neq0}\frac{1}{\epsilon _{0}L^{3}}\frac{1}{k^{2}}e^{-\frac{k^{2}}{4\alpha^{2}}}\left|\sum_{j}q_{j}e^{i\overrightarrow{k}\overrightarrow{r_{j}}}\right|^{2}\quad\quad\alpha 越小收敛越快,k=\frac{2\pi n}{L^{2}} U2(r)=21k=0∑ϵ0L31k21e−4α2k2∣ ∣j∑qjeikrj∣ ∣2α越小收敛越快,k=L22πn - 自能项
U 3 ( r ) = − α π ∑ i q i 2 4 π ϵ 0 U_{3}\left(r\right)=-\frac{\alpha}{\sqrt{\pi}}\sum_{i}\frac{q_{i}^{2}}{4\pi \epsilon _{0}} U3(r)=−παi∑4πϵ0qi2 - 校正项
U 4 ( r ) = 2 π 3 L 3 ∣ ∑ i q i r i → 4 π ϵ 0 ∣ 2 U_{4}\left(r\right)=\frac{2\pi}{3L^{3}}\left|\sum_{i}\frac{q_{i}\overrightarrow{r_{i}}}{4\pi \epsilon _{0}}\right|^{2} U4(r)=3L32π∣ ∣i∑4πϵ0qiri∣ ∣2
- 实空间部分
3.初等分子动力学(NVE)
运动方程的有限差分算法
1. Verlet算法:
将粒子的位置展开为Taylor级数
r
i
→
(
t
+
△
t
)
=
r
i
→
(
t
)
+
v
i
→
(
t
)
△
t
+
1
2
a
i
→
(
t
)
△
t
2
+
⋯
\overrightarrow{r_{i}}\left(t+\bigtriangleup t\right)=\overrightarrow{r_{i}}\left(t\right)+\overrightarrow{v_{i}}\left(t\right)\bigtriangleup t+\frac{1}{2}\overrightarrow{a_{i}}\left(t\right) \triangle t^{2}+\cdots
ri(t+△t)=ri(t)+vi(t)△t+21ai(t)△t2+⋯
r
i
→
(
t
−
△
t
)
=
r
i
→
(
t
)
−
v
i
→
(
t
)
△
t
+
1
2
a
i
→
(
t
)
△
t
2
−
⋯
\overrightarrow{r_{i}}\left(t-\bigtriangleup t\right)=\overrightarrow{r_{i}}\left(t\right)-\overrightarrow{v_{i}}\left(t\right)\bigtriangleup t+\frac{1}{2}\overrightarrow{a_{i}}\left(t\right) \triangle t^{2}-\cdots
ri(t−△t)=ri(t)−vi(t)△t+21ai(t)△t2−⋯
上述两式相加,得
r
i
→
(
t
+
△
t
)
=
−
r
i
→
(
t
−
△
t
)
+
2
r
i
→
(
t
)
+
a
i
→
(
t
)
△
t
2
\overrightarrow{r_{i}}\left(t+\triangle t\right)=-\overrightarrow{r_{i}}\left(t-\triangle t\right)+2\overrightarrow{r_{i}}\left(t\right)+\overrightarrow{a_{i}}\left(t\right)\triangle t^{2}
ri(t+△t)=−ri(t−△t)+2ri(t)+ai(t)△t2
即
r
i
→
n
+
1
=
−
r
i
→
n
+
1
+
2
r
i
→
n
+
h
2
a
i
→
n
\overrightarrow{r_{i}}^{n+1}=-\overrightarrow{r_{i}}^{n+1}+2\overrightarrow{r_{i}}^{n}+h^{2}\overrightarrow{a_{i}}^{n}
rin+1=−rin+1+2rin+h2ain
上述两式相减,得
v
i
→
(
t
)
=
r
i
→
(
t
+
△
t
)
−
r
i
→
(
t
−
△
t
)
2
△
t
\overrightarrow{v_{i}}\left(t\right)=\frac{\overrightarrow{r_{i}}\left(t+\triangle t\right)-\overrightarrow{r_{i}}\left(t-\triangle t\right)}{2\triangle t}
vi(t)=2△tri(t+△t)−ri(t−△t)
即
v
i
→
n
=
1
2
h
(
r
i
→
n
+
1
−
r
i
→
n
−
1
)
\overrightarrow{v_{i}}^{n}=\frac{1}{2h}\left(\overrightarrow{r_{i}}^{n+1}-\overrightarrow{r_{i}}^{n-1}\right)
vin=2h1(rin+1−rin−1)
(初始化):
r
i
→
1
=
−
r
i
→
0
+
h
v
i
→
0
+
h
2
2
m
F
i
→
0
\overrightarrow{r_{i}}^{1}=-\overrightarrow{r_{i}}^{0}+h\overrightarrow{v_{i}}^{0}+\frac{h^{2}}{2m}\overrightarrow{F_{i}}^{0}
ri1=−ri0+hvi0+2mh2Fi0
(迭代计算):
a
i
→
n
=
F
i
→
n
m
i
=
−
1
m
i
∂
U
(
r
→
)
∂
r
i
→
n
\overrightarrow{a_{i}}^{n}=\frac{\overrightarrow{F_{i}}^{n}}{m_{i}}=-\frac{1}{m_{i}}\frac{\partial U\left(\overrightarrow{r}\right)}{\partial \overrightarrow{r_{i}}^{n}}
ain=miFin=−mi1∂rin∂U(r)
r
i
→
n
+
1
=
−
r
i
→
n
+
1
+
2
r
i
→
n
+
h
2
a
i
→
n
\overrightarrow{r_{i}}^{n+1}=-\overrightarrow{r_{i}}^{n+1}+2\overrightarrow{r_{i}}^{n}+h^{2}\overrightarrow{a_{i}}^{n}
rin+1=−rin+1+2rin+h2ain
v
i
→
n
=
1
2
h
(
r
i
→
n
+
1
−
r
i
→
n
−
1
)
\overrightarrow{v_{i}}^{n}=\frac{1}{2h}\left(\overrightarrow{r_{i}}^{n+1}-\overrightarrow{r_{i}}^{n-1}\right)
vin=2h1(rin+1−rin−1)
2. 蛙跳算法
将粒子的位置展开为Taylor级数
r
i
→
(
t
+
△
t
)
=
r
i
→
(
t
)
+
v
i
→
(
t
)
△
t
+
1
2
a
i
→
(
t
)
△
t
2
+
⋯
\overrightarrow{r_{i}}\left(t+\bigtriangleup t\right)=\overrightarrow{r_{i}}\left(t\right)+\overrightarrow{v_{i}}\left(t\right)\bigtriangleup t+\frac{1}{2}\overrightarrow{a_{i}}\left(t\right) \triangle t^{2}+\cdots
ri(t+△t)=ri(t)+vi(t)△t+21ai(t)△t2+⋯
r
i
→
(
t
−
△
t
)
=
r
i
→
(
t
)
−
v
i
→
(
t
)
△
t
+
1
2
a
i
→
(
t
)
△
t
2
−
⋯
\overrightarrow{r_{i}}\left(t-\bigtriangleup t\right)=\overrightarrow{r_{i}}\left(t\right)-\overrightarrow{v_{i}}\left(t\right)\bigtriangleup t+\frac{1}{2}\overrightarrow{a_{i}}\left(t\right) \triangle t^{2}-\cdots
ri(t−△t)=ri(t)−vi(t)△t+21ai(t)△t2−⋯
对两式分别求导得到
v
i
→
(
t
+
△
t
)
=
v
i
→
(
t
)
+
a
i
→
(
t
)
△
t
+
1
2
a
i
→
(
t
)
△
t
2
+
⋯
\overrightarrow{v_{i}}\left(t+\bigtriangleup t\right)=\overrightarrow{v_{i}}\left(t\right)+\overrightarrow{a_{i}}\left(t\right)\bigtriangleup t+\frac{1}{2}\overrightarrow{a_{i}}\left(t\right) \triangle t^{2}+\cdots
vi(t+△t)=vi(t)+ai(t)△t+21ai(t)△t2+⋯
v
i
→
(
t
−
△
t
)
=
v
i
→
(
t
)
−
a
i
→
(
t
)
△
t
+
1
2
a
i
→
(
t
)
△
t
2
−
⋯
\overrightarrow{v_{i}}\left(t-\bigtriangleup t\right)=\overrightarrow{v_{i}}\left(t\right)-\overrightarrow{a_{i}}\left(t\right)\bigtriangleup t+\frac{1}{2}\overrightarrow{a_{i}}\left(t\right) \triangle t^{2}-\cdots
vi(t−△t)=vi(t)−ai(t)△t+21ai(t)△t2−⋯
上述两式相减,得
v
i
→
(
t
+
△
t
)
=
v
i
→
(
t
−
△
t
)
+
2
a
i
→
(
t
)
△
t
\overrightarrow{v_{i} } \left ( t+\triangle t\right ) =\overrightarrow{v_{i} } \left ( t-\triangle t\right )+2\overrightarrow{a_{i} } \left ( t \right ) \triangle t
vi(t+△t)=vi(t−△t)+2ai(t)△t
令
△
t
→
△
t
2
,
则
v
i
→
n
+
1
2
=
v
i
→
n
−
1
2
+
h
a
i
→
n
令\triangle t\to\frac{\triangle t}{2},则\overrightarrow{v_{i}}^{n+\frac{1}{2}}=\overrightarrow{v_{i}}^{n-\frac{1}{2}}+h\overrightarrow{a_{i}}^{n}
令△t→2△t,则vin+21=vin−21+hain
对位置矢量重新整理,得
r
i
→
(
t
+
△
t
)
=
r
i
→
(
t
)
+
△
t
(
v
i
→
(
t
)
+
1
2
a
i
→
(
t
)
△
t
)
=
r
i
→
(
t
)
+
△
t
⋅
v
i
→
(
t
+
△
t
2
)
\overrightarrow{r_{i}}\left(t+\triangle t\right)=\overrightarrow{r_{i}}\left(t\right)+\triangle t\left(\overrightarrow{v_{i}}\left(t\right)+\frac{1}{2}\overrightarrow{a_{i}}\left(t\right)\triangle t\right)=\overrightarrow{r_{i}}\left(t\right)+\triangle t\cdot\overrightarrow{v_{i}}\left(t+\frac{\triangle t}{2}\right)
ri(t+△t)=ri(t)+△t(vi(t)+21ai(t)△t)=ri(t)+△t⋅vi(t+2△t)
即
r
i
→
n
+
1
=
r
i
→
n
+
h
v
i
→
n
+
1
2
\overrightarrow{r_{i}}^{n+1}=\overrightarrow{r_{i}}^{n}+h\overrightarrow{v_{i}}^{n+\frac{1}{2}}
rin+1=rin+hvin+21
(初始化):
v
i
→
−
1
2
=
v
i
→
0
−
h
2
m
i
F
i
→
0
\overrightarrow{v_{i}}^{-\frac{1}{2}}=\overrightarrow{v_{i}}^{0}-\frac{h}{2m_{i}}\overrightarrow{F_{i}}^{0}
vi−21=vi0−2mihFi0
(迭代计算):
a
i
→
n
=
F
i
→
n
m
i
=
−
1
m
i
∂
U
(
r
→
)
∂
r
i
→
n
\overrightarrow{a_{i}}^{n}=\frac{\overrightarrow{F_{i}}^{n}}{m_{i}}=-\frac{1}{m_{i}}\frac{\partial U\left(\overrightarrow{r}\right)}{\partial \overrightarrow{r_{i}}^{n}}
ain=miFin=−mi1∂rin∂U(r)
v
i
→
n
+
1
2
=
v
i
→
n
−
1
2
+
h
a
i
→
n
\overrightarrow{v_{i}}^{n+\frac{1}{2}}=\overrightarrow{v_{i}}^{n-\frac{1}{2}}+h\overrightarrow{a_{i}}^{n}
vin+21=vin−21+hain
r
i
→
n
+
1
=
r
i
→
n
+
h
v
i
→
n
+
1
2
\overrightarrow{r_{i}}^{n+1}=\overrightarrow{r_{i}}^{n}+h\overrightarrow{v_{i}}^{n+\frac{1}{2}}
rin+1=rin+hvin+21
v
i
→
n
=
1
2
(
v
i
→
n
+
1
2
+
v
i
→
n
−
1
2
)
\overrightarrow{v_{i}}^{n}=\frac{1}{2}\left(\overrightarrow{v_{i}}^{n+\frac{1}{2}}+\overrightarrow{v_{i}}^{n-\frac{1}{2}}\right)
vin=21(vin+21+vin−21)
3. 速度Verlet算法(lammps中默认使用)
将粒子的位置展开为Taylor级数
r
i
→
(
t
+
△
t
)
=
r
i
→
(
t
)
+
v
i
→
(
t
)
△
t
+
1
2
a
i
→
(
t
)
△
t
2
+
⋯
\overrightarrow{r_{i}}\left(t+\bigtriangleup t\right)=\overrightarrow{r_{i}}\left(t\right)+\overrightarrow{v_{i}}\left(t\right)\bigtriangleup t+\frac{1}{2}\overrightarrow{a_{i}}\left(t\right) \triangle t^{2}+\cdots
ri(t+△t)=ri(t)+vi(t)△t+21ai(t)△t2+⋯
对上式求导可得
v
i
→
(
t
+
△
t
)
=
v
i
→
(
t
)
+
a
i
→
(
t
)
△
t
+
⋯
=
v
i
→
(
t
)
+
1
2
[
a
i
→
(
t
)
+
a
i
→
(
t
+
△
t
)
]
△
t
+
⋯
\overrightarrow{v_{i}}\left(t+\triangle t\right) =\overrightarrow{v_{i}}\left(t\right)+\overrightarrow{a_{i}}\left(t\right)\triangle t+\cdots\\ \quad\quad\quad\quad\quad\quad\quad=\overrightarrow{v_{i}}\left(t\right)+\frac{1}{2}\left[\overrightarrow{a_{i}}\left(t\right)+\overrightarrow{a_{i}}\left(t+\triangle t\right)\right]\triangle t+\cdots
vi(t+△t)=vi(t)+ai(t)△t+⋯=vi(t)+21[ai(t)+ai(t+△t)]△t+⋯
上述两式等价于
v
i
→
(
t
+
△
t
2
)
=
v
i
→
(
t
)
+
1
2
a
i
→
(
t
)
△
t
\overrightarrow{v_{i}}\left(t+\frac{\triangle t}{2}\right)=\overrightarrow{v_{i}}\left(t\right)+\frac{1}{2}\overrightarrow{a_{i}}\left(t\right)\triangle t
vi(t+2△t)=vi(t)+21ai(t)△t
r
i
→
(
t
+
△
t
)
=
r
i
→
(
t
)
+
v
i
→
(
t
+
△
t
2
)
△
t
\overrightarrow{r_{i}}\left(t+\triangle t\right)=\overrightarrow{r_{i}}\left(t\right)+\overrightarrow{v_{i}}\left(t+\frac{\triangle t}{2}\right)\triangle t
ri(t+△t)=ri(t)+vi(t+2△t)△t
v
i
→
(
t
+
△
t
)
=
v
i
→
(
t
+
△
t
2
)
+
1
2
a
i
→
(
t
+
△
t
)
△
t
\overrightarrow{v_{i}}\left(t+\triangle t\right)=\overrightarrow{v_{i}}\left(t+\frac{\triangle t}{2}\right)+\frac{1}{2}\overrightarrow{a_{i}}\left(t+\triangle t\right)\triangle t
vi(t+△t)=vi(t+2△t)+21ai(t+△t)△t
即
v
i
→
n
+
1
2
=
v
i
→
n
+
h
2
a
i
→
n
\overrightarrow{v_{i}}^{n+\frac{1}{2}}=\overrightarrow{v_{i}}^{n}+\frac{h}{2}\overrightarrow{a_{i}}^{n}
vin+21=vin+2hain
r
i
→
n
+
1
=
r
i
→
n
+
h
v
i
→
n
+
1
2
\overrightarrow{r_{i}}^{n+1}=\overrightarrow{r_{i}}^{n}+h\overrightarrow{v_{i}}^{n+\frac{1}{2}}
rin+1=rin+hvin+21
v
i
→
n
+
1
=
v
i
→
n
+
1
2
+
h
2
a
i
→
n
+
1
\overrightarrow{v_{i}}^{n+1}=\overrightarrow{v_{i}}^{n+\frac{1}{2}}+\frac{h}{2}\overrightarrow{a_{i}}^{n+1}
vin+1=vin+21+2hain+1
(迭代计算):
v
i
→
n
+
1
2
=
v
i
→
n
+
h
2
a
i
→
n
\overrightarrow{v_{i}}^{n+\frac{1}{2}}=\overrightarrow{v_{i}}^{n}+\frac{h}{2}\overrightarrow{a_{i}}^{n}
vin+21=vin+2hain
r
i
→
n
+
1
=
r
i
→
n
+
h
v
i
→
n
+
1
2
\overrightarrow{r_{i}}^{n+1}=\overrightarrow{r_{i}}^{n}+h\overrightarrow{v_{i}}^{n+\frac{1}{2}}
rin+1=rin+hvin+21
a
i
→
n
+
1
=
F
i
→
n
m
i
=
−
1
m
i
∂
U
(
r
→
)
∂
r
i
→
n
+
1
\overrightarrow{a_{i}}^{n+1}=\frac{\overrightarrow{F_{i}}^{n}}{m_{i}}=-\frac{1}{m_{i}}\frac{\partial U\left(\overrightarrow{r}\right)}{\partial \overrightarrow{r_{i}}^{n+1}}
ain+1=miFin=−mi1∂rin+1∂U(r)
v
i
→
n
+
1
=
v
i
→
n
+
1
2
+
h
2
a
i
→
n
+
1
\overrightarrow{v_{i}}^{n+1}=\overrightarrow{v_{i}}^{n+\frac{1}{2}}+\frac{h}{2}\overrightarrow{a_{i}}^{n+1}
vin+1=vin+21+2hain+1
积分步长的选择
过长的步长会造成分子间的剧烈碰撞,导致数据溢出;过小的步长会降低模拟过程中搜索相空间的能力。积分步长需小于系统中最快周期的十分之一。
约化单位制
采用国际单位制会导致物理量的数值非常小,在计算中会引起严重的舍入误差,因此,实际模拟时通常采用约化单位制。
常见物理量的约化表示:
能量
:
E
∗
=
E
ε
数密度
:
ρ
∗
=
ρ
σ
3
速率
:
v
∗
=
v
m
ε
黏度
:
η
∗
=
η
σ
2
m
ε
长度
:
r
∗
=
r
σ
温度
:
T
∗
=
k
B
T
ε
力
:
f
∗
=
f
σ
ε
电荷
:
q
∗
=
q
4
π
ϵ
0
σ
ε
时间
:
t
∗
=
t
ε
m
σ
2
压强
:
p
∗
=
p
σ
3
ε
力矩
:
τ
∗
=
τ
ε
偶极矩
:
μ
∗
=
μ
4
π
ϵ
0
σ
3
ε
能量:E^{*}=\frac{E}{\varepsilon}\quad\quad数密度:\rho ^{*} =\rho \sigma ^{3}\quad\quad速率:v^{*}=v\sqrt{\frac{m}{\varepsilon}}\quad\quad黏度:\eta^{*}=\eta\frac{\sigma^{2}}{\sqrt{m\varepsilon}} \\ 长度:r^{*}=\frac{r}{\sigma}\quad\quad\quad温度:T^{*}=\frac{k_{B}T}{\varepsilon}\quad\quad力:f^{*}=f\frac{\sigma}{\varepsilon}\quad\quad\quad\quad电荷:q^{*}=\frac{q}{\sqrt{4\pi\epsilon^{0}\sigma\varepsilon }}\\ 时间:t^{*}=t\sqrt{\frac{\varepsilon }{m\sigma^{2}}}\quad压强:p^{*}=p\frac{\sigma^{3}}{\varepsilon}\quad\quad力矩:\tau ^{*} =\frac{\tau }{\varepsilon } \quad\quad\quad\quad偶极矩:\mu ^{*}=\frac{\mu}{\sqrt{4\pi\epsilon^{0}\sigma^{3}\varepsilon }}
能量:E∗=εE数密度:ρ∗=ρσ3速率:v∗=vεm黏度:η∗=ηmεσ2长度:r∗=σr温度:T∗=εkBT力:f∗=fεσ电荷:q∗=4πϵ0σεq时间:t∗=tmσ2ε压强:p∗=pεσ3力矩:τ∗=ετ偶极矩:μ∗=4πϵ0σ3εμ
由于文章篇幅有限,剩下内容请点击[分子随记]分子动力学模拟基础(二)