Alder, Berni Julian, and Thomas Everett Wainwright. “Phase transition for a hard sphere system.” The Journal of chemical physics 27.5 (1957): 1208-1209.
2.第一个连续势场的分子动力学模拟
在1960年，连续排斥的Born-Mayer相互作用势能首次用于MD在布鲁克海文国家实验室进行的铜靶辐射损伤模拟(J.B. Gibson，Goland，M.Milgram和G.H.Vineyard，Dynamics of radiation damage，Phys. Rev. 120，1229-1253，1960)。这可能是材料科学中MD方法的第一个应用。
计算晶胞由446至998个铜原子组成。一个时间积分步骤在IBM704电脑上大约需要一分钟。
Gibson, J. B., et al. “Dynamics of radiation damage.” Physical Review 120.4 (1960): 1229.
Klaus Schulten, University of Illinois at Urbana-Champaign
DL_POLY
商业
W. Smith and T.R. Forester,Daresbury Laboratory
Materials Studio
商业
Biovia
3.分子动力学的步骤
1.输入信息
1.1 初始化-位置和速度
初始速度： 初始速度一般按玻尔兹曼分布或者高斯分布取得，通常计算前检查粒子总动量为零，否则计算的系统本身产生移动而导致总能量的不稳定。 Maxwell-Boltzmann速率分布
f
(
v
)
=
4
π
(
m
2
k
T
)
3
2
v
2
e
−
m
v
2
2
k
T
f(v)=\frac{4}{\sqrt{\pi}}(\frac{m}{2kT})^{\frac{3}{2}}v^{2}e^{-\frac{mv^{2}}{2kT}}
f(v)=π4(2kTm)23v2e−2kTmv2
动力学演化： N个原子组成的分子体系，第i个原子运动的牛顿方程是：
f
i
=
m
i
d
2
r
i
d
t
2
=
m
i
r
ˉ
i
…
(
i
=
1
,
2
,
3
,
.
.
.
,
N
)
f_{i}=m_{i}\frac{d^{2}r_{i}}{dt^{2}}=m_{i}\bar{r} _{i}\dots(i=1,2,3,...,N)
fi=midt2d2ri=mirˉi…(i=1,2,3,...,N)对其速度进行积分即为粒子的位置轨迹：
r
(
t
)
=
∫
0
t
v
(
t
)
d
t
r(t)=\int_{0}^{t}v(t)dt
r(t)=∫0tv(t)dt
r
(
t
+
Δ
t
)
=
r
(
t
)
+
v
(
t
)
Δ
t
+
O
(
Δ
t
)
2
r(t+\Delta t)=r(t)+v(t)\Delta t+O(\Delta t)^{2}
r(t+Δt)=r(t)+v(t)Δt+O(Δt)2
v
(
t
+
Δ
t
)
=
v
(
t
)
+
F
(
t
)
m
Δ
t
+
O
(
Δ
t
)
2
v(t+\Delta t)=v(t)+\frac{F(t)}{m}\Delta t+O(\Delta t)^{2}
v(t+Δt)=v(t)+mF(t)Δt+O(Δt)2 通过算法对下一时刻原子速度和位置进行求解，算法主要包括：Euler算法、Velert算法、Leap-frog算法、Velocity-velert算法
T
=
1
3
1
N
k
B
⟨
∑
i
=
1
N
m
i
v
⃗
i
2
⟩
v
⃗
i
2
=
v
⃗
i
⋅
v
⃗
i
T=\frac{1}{3}\frac{1}{Nk_{B}}\left \langle \sum_{i=1}^{N} m_{i} \vec{v}_{i}^{2}\right \rangle \quad \vec{v}_{i}^{2}=\vec{v}_{i}\cdot \vec{v}_{i}
T=31NkB1⟨∑i=1Nmivi2⟩vi2=vi⋅vi
Direct
Pressure
P
=
1
3
V
⟨
∑
i
=
1
N
(
m
i
v
⃗
i
2
+
r
⃗
i
⋅
f
⃗
i
)
⟩
P=\frac{1}{3V}\left \langle \sum_{i=1}^{N}(m_{i}\vec{v}_{i}^{2}+\vec{r}_{i}\cdot \vec{f}_{i}) \right \rangle
P=3V1⟨∑i=1N(mivi2+ri⋅fi)⟩
Direct
Stress
σ
i
j
=
1
Ω
⟨
(
−
∑
α
m
α
u
α
,
i
u
α
,
j
+
1
2
∑
α
,
β
,
α
≠
β
∂
ϕ
(
r
)
∂
r
r
i
r
⋅
r
j
∣
r
=
r
α
β
)
⟩
\sigma _{ij}=\frac{1}{\Omega}\left\langle (-\sum_{\alpha}m_{\alpha}u_{\alpha,i}u_{\alpha,j}+\frac{1}{2}\sum_{\alpha,\beta,\alpha\ne\beta}\frac{\partial \phi (r) }{\partial r}\frac{r_{i}}{r}\cdot r_{j}\mid_{r=r_{\alpha\beta}}) \right\rangle
σij=Ω1⟨(−∑αmαuα,iuα,j+21∑α,β,α=β∂r∂ϕ(r)rri⋅rj∣r=rαβ)⟩
Direct
MSD
⟨
Δ
r
2
(
t
)
⟩
=
1
N
∑
i
(
r
i
(
t
)
−
r
i
(
t
=
0
)
)
2
\left\langle \Delta r^{2}(t) \right\rangle=\frac{1}{N}\sum_{i}(r_{i}(t)-r_{i}(t=0))^{2}
⟨Δr2(t)⟩=N1∑i(ri(t)−ri(t=0))2
Diffusivity
RDF
g
(
r
)
=
⟨
N
(
r
±
Δ
r
2
)
Ω
(
r
±
Δ
r
2
)
ρ
⟩
g(r)=\left\langle \frac{N(r\pm \frac{\Delta r}{2})}{\Omega(r\pm \frac{\Delta r}{2})\rho} \right\rangle
g(r)=⟨Ω(r±2Δr)ρN(r±2Δr)⟩
Atomic structure(signature)
VAF
⟨
v
(
0
)
v
(
t
)
⟩
=
1
N
∑
i
=
1
N
1
N
i
∑
k
=
1
N
i
v
i
(
t
k
)
v
i
(
t
k
+
t
)
\left\langle v(0)v(t) \right\rangle=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{N_{i}}\sum_{k=1}^{N_i}v_{i}(t_{k})v_{i}(t_{k}+t)
⟨v(0)v(t)⟩=N1∑i=1NNi1∑k=1Nivi(tk)vi(tk+t)
此近似有效性的一个重要判断依据是德布罗意波长：当德布罗意波长
Λ
\Lambda
Λ远大于粒子间距
d
d
d的时候，量子效应将会非常的显著
对于热运动，可求得德布罗意波长，T=300K平衡态时
Λ
t
h
=
1
2
π
m
k
B
T
\Lambda_{th}=\frac{1}{\sqrt{2\pi mk_{B}T}}
Λth=2πmkBT1
Λ
t
h
=
1
A
˚
f
o
r
a
H
a
t
o
m
(
m
H
=
1
a
m
u
)
\Lambda_{th}=1\AA \; \mathrm{for \; a \; H \; atom \; (m_{H}=1amu)}
Λth=1A˚foraHatom(mH=1amu)
Λ
t
h
=
0.19
A
˚
f
o
r
a
S
i
a
t
o
m
(
m
S
i
=
28
a
m
u
)
\Lambda_{th}=0.19\AA \; \mathrm{for \; a \; Si \; atom \; (m_{Si}=28amu)}
Λth=0.19A˚foraSiatom(mSi=28amu)
Λ
t
h
=
0.07
A
˚
f
o
r
a
A
u
a
t
o
m
(
m
A
u
=
197
a
m
u
)
\Lambda_{th}=0.07\AA \; \mathrm{for \; a \; Au \; atom \; (m_{Au}=197amu)}
Λth=0.07A˚foraAuatom(mAu=197amu)
所以在绝大多数足够高温度的情况下，除了H、He、Ne等轻原子外，原子可以用质点来描述(
d
>
>
Λ
d>>\Lambda
d>>Λ)，可以用牛顿方程来描述原子的运动
扫一扫
1008
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
