分子动力学模拟基础(一)

雪薇儿

已于 2022-07-24 17:33:45 修改

阅读量2.6k

点赞数 4

文章标签：科技

于 2022-04-28 21:31:26 首次发布

本文链接：https://blog.csdn.net/qq_44963682/article/details/124434098

版权

文章目录

1.分子模拟技术

分子模拟的分类

Monte Carlo 模拟:1949年，利用随机行走产生微观态。通过系综平均，得到体系的平衡性质。
分子动力学模拟：1957年，利用运动轨迹产生微观态。通过时间平均，得到体系的动力学性质。
混合方法模拟：QM/MM方法，QM+MD方法，MC+MD方法等。

2.分子动力学的基本概念

分子动力学

初始化条件:势函数，相互作用，系综，外界条件，初始构型
动力学方程：Newton力学方程
一次信息：原子的坐标，三维结构，原子的速度和加速度，原子运动
二次信息:热力学性质，动力学性质，光学性质

初始化条件——模拟盒子

边界条件:
在小体系中，边界效应一般不能忽略。
如果研究的体系受环境影响较小，可采用自然边界条件；如果研究的体系会受到环境的影响，则最好采用周期边界条件。
周期边界条件：
中心盒子在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 为下取整$
利用周期边界条件可以基于少量粒子模拟的同时考虑周围环境分子多研究体系的影响。
势函数(力场)：
描述分子内相互作用和分子间相互作用的数学模型。
分子内相互作用：
(1)键相互作用，例如bond_style命令
$U_{bond}\left ( r \right )=k_{r} \left ( r-r_{0} \right )^{2}$
(2)角相互作用，例如angle_style命令
$U_{angle}\left ( \theta \right )=k_{\theta } \left ( \theta -\theta _{0} \right )^{2}$
(3)二面角相互作用，例如dihedral_style命令
$U_{dihedral}\left ( \varphi \right )=k_{\varphi}\left [ 1\pm cos\left ( n\varphi \right ) \right ]$
(4)离面角相互作用，例如improper_style命令
$U_{improper}\left (\chi \right )=k_{\chi }\left ( \chi -\chi _{0} \right )^{2}$
分子间相互作用
(1)van der Waals作用，例如pair_style lj命令
$U\left (r_{ij}\right)=4\varepsilon \left [ \left ( \frac{\sigma }{r_{ij} } \right ) ^{12} -\left ( \frac{\sigma }{r_{ij} } \right )^{6} \right ]$
(2)静电相互作用，例如pair_style coul命令
$U\left (r_{ij}\right)=k\frac{q_{i}q_{j} }{\epsilon r_{ij} }$
(3)三体相互作用，例如pair_style tersoff命令
$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 ]$
其中
$\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} }$
(4)多体相互作用，例如pair_style eam命令
$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 )$

初始化条件——短程势函数的处理

当分子间相互作用使积分
$\int_{r_{cut}}^{\infty} U\left ( r \right ) 4\pi r^{2}dr$ 收敛，则为短程相互作用,否则为长程相互作用。

最小映像原理：
中心盒子中的一个粒子只与此盒子中的其他N-1个粒子，或它们的最近映像发生相互作用。
$\mu_{ij} = \mu_{ij} -\left [ \frac{\mu_{ij} }{L_{\mu } } \right ]\quad\left ( \mu =x,y,z \right ) \quad\left [\quad \right ] 为就近取整$
势能截断：
不可能将粒子与它们的全部映像间的相互作用都进行计算，势函数可以在一定距离r_{cut}处进行截断，一般势函数的截断半径不大于模拟盒子的最短边长的一半。
(1)简单截断
势能在截断处不连续(分子跨过截断半径时总能量不守恒)，力在截断处为无穷大(MD过程不稳定)
$\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.$
(2)位移截断
势能在截断处连续，力在截断处为有限值但不连续。
$\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.$
(3)位移—力截断
势能在截断处连续，力在截断处连续
$\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.$

初始化条件——长程势函数的处理

Ewald加和：
长程相互作用(如静电相互作用)需要将粒子与它们的全部映像间的相互作用都进行计算。
$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为中心盒子，否则为映像盒子$
直接计算静电相互作用收敛慢，故采用Ewald加和算法：在每个点电荷周围加一个各向同性、数值相同、符号相反的Gauss型电荷(屏蔽电荷) $\rho _{i}\left ( r \right )=\frac{q_{i}\alpha ^{3}}{\pi ^{\frac{3}{2} }} e^{-\alpha ^{2}r^{2}}$
将静电相互作用分为一个短程力和一个长程力。短程力在实空间做截断，长程力在倒易空间做截断。
1. 实空间部分
  $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 越大收敛越快$
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}}$
3. 自能项
  $U_{3}\left(r\right)=-\frac{\alpha}{\sqrt{\pi}}\sum_{i}\frac{q_{i}^{2}}{4\pi \epsilon _{0}}$
4. 校正项
  $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}$

3.初等分子动力学(NVE)

运动方程的有限差分算法

1. Verlet算法:
将粒子的位置展开为Taylor级数
$\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$
$\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$
上述两式相加，得
$\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}$
即 $\overrightarrow{r_{i}}^{n+1}=-\overrightarrow{r_{i}}^{n+1}+2\overrightarrow{r_{i}}^{n}+h^{2}\overrightarrow{a_{i}}^{n}$
上述两式相减，得
$\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}$
即 $\overrightarrow{v_{i}}^{n}=\frac{1}{2h}\left(\overrightarrow{r_{i}}^{n+1}-\overrightarrow{r_{i}}^{n-1}\right)$
(初始化): $\overrightarrow{r_{i}}^{1}=-\overrightarrow{r_{i}}^{0}+h\overrightarrow{v_{i}}^{0}+\frac{h^{2}}{2m}\overrightarrow{F_{i}}^{0}$
(迭代计算): $\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}}$
$\overrightarrow{r_{i}}^{n+1}=-\overrightarrow{r_{i}}^{n+1}+2\overrightarrow{r_{i}}^{n}+h^{2}\overrightarrow{a_{i}}^{n}$
$\overrightarrow{v_{i}}^{n}=\frac{1}{2h}\left(\overrightarrow{r_{i}}^{n+1}-\overrightarrow{r_{i}}^{n-1}\right)$
2. 蛙跳算法
将粒子的位置展开为Taylor级数
$\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$
$\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$
对两式分别求导得到
$\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$
$\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$
上述两式相减，得
$\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$
$令\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}$
对位置矢量重新整理，得
$\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)$
即 $\overrightarrow{r_{i}}^{n+1}=\overrightarrow{r_{i}}^{n}+h\overrightarrow{v_{i}}^{n+\frac{1}{2}}$
(初始化): $\overrightarrow{v_{i}}^{-\frac{1}{2}}=\overrightarrow{v_{i}}^{0}-\frac{h}{2m_{i}}\overrightarrow{F_{i}}^{0}$
(迭代计算):
$\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}}$
$\overrightarrow{v_{i}}^{n+\frac{1}{2}}=\overrightarrow{v_{i}}^{n-\frac{1}{2}}+h\overrightarrow{a_{i}}^{n}$
$\overrightarrow{r_{i}}^{n+1}=\overrightarrow{r_{i}}^{n}+h\overrightarrow{v_{i}}^{n+\frac{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)$
3. 速度Verlet算法(lammps中默认使用)
将粒子的位置展开为Taylor级数
$\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$
对上式求导可得
$\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$
上述两式等价于
$\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$
$\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$
$\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$
即
$\overrightarrow{v_{i}}^{n+\frac{1}{2}}=\overrightarrow{v_{i}}^{n}+\frac{h}{2}\overrightarrow{a_{i}}^{n}$
$\overrightarrow{r_{i}}^{n+1}=\overrightarrow{r_{i}}^{n}+h\overrightarrow{v_{i}}^{n+\frac{1}{2}}$
$\overrightarrow{v_{i}}^{n+1}=\overrightarrow{v_{i}}^{n+\frac{1}{2}}+\frac{h}{2}\overrightarrow{a_{i}}^{n+1}$
(迭代计算):
$\overrightarrow{v_{i}}^{n+\frac{1}{2}}=\overrightarrow{v_{i}}^{n}+\frac{h}{2}\overrightarrow{a_{i}}^{n}$
$\overrightarrow{r_{i}}^{n+1}=\overrightarrow{r_{i}}^{n}+h\overrightarrow{v_{i}}^{n+\frac{1}{2}}$
$\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}}$
$\overrightarrow{v_{i}}^{n+1}=\overrightarrow{v_{i}}^{n+\frac{1}{2}}+\frac{h}{2}\overrightarrow{a_{i}}^{n+1}$

积分步长的选择

过长的步长会造成分子间的剧烈碰撞，导致数据溢出；过小的步长会降低模拟过程中搜索相空间的能力。积分步长需小于系统中最快周期的十分之一。

约化单位制

采用国际单位制会导致物理量的数值非常小，在计算中会引起严重的舍入误差，因此，实际模拟时通常采用约化单位制。
常见物理量的约化表示:
$能量: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 }}$