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

雪薇儿

已于 2022-07-24 17:34:35 修改

阅读量2.3k

点赞数 6

文章标签：科技

于 2022-04-28 21:32:04 首次发布

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

版权

没有看过上一篇的小伙伴，可以点击上篇[分子随记]分子动力学模拟基础(一)

文章目录

- 4.限制性分子动力学(NVT、NpT、NpH...)

4.限制性分子动力学(NVT、NpT、NpH…)

系综

在一定宏观条件下，大量宏观状态相同而微观状态不同的独立系统的集合，称为统计系综。

	守恒量与特征函数	宏观条件	典型适用场所
微正则系综	$S (N, V, E)$	孤立系统	输运过程
正则系综	$F (N, V, T)$	封闭系统	凝聚系统过程
巨正则系统	$J(\mu,V,T)$	开放系统	吸附过程
等温等压系综	$G (N, p, T)$		相变过程
等压等焓系综	$\wp(N,p,H)$		节流过程

NVT系综分子动力学

温度控制
在NVE系综的分子动力学基础上，通过温度控制使系统的温度维持在给定值，或根据环境温度使系统温度发生涨落，即可实现NVT系综分子动力学。
(1)Woodcock速度标度法(fix temp/rescale命令):
根据能量均分定理进行标度
$\frac{1}{2}\sum_{i=1}^{N}m_{i}\overrightarrow{v_{i}}^{2}=\frac{3}{2}Nk_{B}T$
$\lambda=\sqrt{\frac{\left(3N-3\right)k_{B}\tilde{T} }{\sum_{i=1}^{N}m_{i}\overrightarrow{v_{i}}^{2}\left(t\right)}}=\sqrt{\frac{\tilde{T} }{T\left(t\right)}}\quad\quad\quad\overrightarrow{v_{i}}\left(t\right)\to\lambda\overrightarrow{v_{i}}\left(t\right)$
(2)Berendsen热浴标度法(fix temp/berendsen命令)：
根据Fourier定律进行标度
$\frac{dT\left(t\right)}{dt}=\frac{1}{\tau}\left[\tilde{T}-T\left(t\right)\right]$
$\lambda=\sqrt{1+\left(\frac{\tilde{T}}{T\left(t\right)}-1\right)\frac{\triangle t}{\tau}}\quad\quad\quad\overrightarrow{v_{i}}\left(t+\triangle t\right)\to\lambda\overrightarrow{v_{i}}\left(t\right)$
改进版本(fix temp/csvr命令):在Fourier定律基础上引入随机项
$\frac{dT\left(t\right)}{dt}=\frac{1}{\tau}\left[\tilde{T}-T\left(t\right)\right]+\frac{2}{\sqrt{\tau}}\sqrt{\frac{\tilde{T}T\left(t\right)}{3N}}\frac{dW\left(t\right)}{dt}$
(3)Andersen热浴法:
系统中的粒子以一定的频率v与温度为T的热浴进行碰撞，碰撞后粒子从温度为T的Maxwell-Boltzmann分布中随机获得一个速度。
$v=\frac{\kappa}{\sqrt[3]{\rho N^{2}}}$
(4)限制控温法/阻尼力法/等动能法(fix nvk命令):
根据恒温动力学方程 $\left\{\begin{array}{l}\frac{d\overrightarrow{r_{i}}}{dt}=\overrightarrow{v_{i}}\left(t\right) \\ \frac{d\overrightarrow{v_{i}}}{dt}=\frac{\overrightarrow{F_{i}}\left(t\right)}{m_{i}}-\xi\left(t\right)\overrightarrow{v_{i}}\left(t\right) \end{array}\right.$
引入温度限制 $\frac{dT}{dt}=0$ (即恒温)，则待定系数 $\xi(t)$ 的选择应使动力学方程轨迹与温度限制轨迹的差异最小化。
由Gauss最小约束原理，有
$\xi\left(t\right)=\frac{\sum_{i}m_{i}\overrightarrow{v_{i}}\left(t\right)\cdot\overrightarrow{F_{i}}\left(t\right)}{\sum_{i}m_{i}^{2}\overrightarrow{v_{i}}^{2}\left(t\right)}$
(5)Nose-Hoover热浴法(fix nvt命令):
根据恒温动力学方程，增加额外的自由度，在扩展的微正则系综中完成体系的正则系综模拟，则待定系数 $\xi(t)$ 满足如下方程。
$\frac{d\xi\left(t\right)}{dt}=\frac{gk_{B}}{Q}\left[T\left(t\right)-\tilde{T}\right]=\frac{1}{Q}\left[\sum_{i}m_{i}\overrightarrow{v_{i}}^{2}\left(t\right)-gk_{B}\tilde{T}\right]\quad\quad\quad其中Q=gk_{B}\tilde{T}\tau^{2}$
(6)Langevin动力学法(fix langevin命令):
执行Brown动力学，每个粒子受力
$\overrightarrow{F_{i}}=\sum_{j}\left[-\bigtriangledown_{i}U\left(r_{ij}\right)-\frac{m_{i}}{\tau}\left(\overrightarrow{v_{ij}}\cdot\hat{\overrightarrow{r_{ij}} } \right)\hat{\overrightarrow{r_{ij}} } +\sqrt{2k_{B}\tilde{T}\frac{m_{i}}{\tau}}\frac{dW\left(t\right)}{dt}\hat{\overrightarrow{r_{ij}} } \right]$
改进版本(fix temp/csld命令)：对粒子速度的更新是在原速度上线性叠加一个温度为 $\overrightarrow{T}$ 的Gauss分布的随机速度，来实现控温。
(7) 耗散粒子动力学法(pair_style dpd/tstat命令):
与Langvin动力学法原理相同，不同的是耗散力和随机力在Langvin动力学法中通过fix命令施加的，而在耗散粒子动力学中是通过pair_style命令由力场施加的。
附加到粒子上的耗散力和随机力为
$\overrightarrow{F_{i}}=\begin{cases}\sum_{j}\left[-\gamma\left(1-\frac{r_{ij}}{r_{cut}}\right)^{2}\left(\overrightarrow{v_{ij}}\cdot\hat{\overrightarrow{r_{ij}} } \right)\hat{\overrightarrow{r_{ij}} } +\sqrt{2k_{B}\tilde{T}_{\gamma }}\left(1-\frac{r_{ij}}{r_{cut}}\right)\frac{dW\left(t\right)}{dt}\hat{\overrightarrow{r_{ij}} }\right]\quad&\left(r_{ij}\le r_{cut}\right) \\0&\left(r_{ij}> r_{cut}\right) \end{cases}$

NpT系综分子动力学

压强控制
在NVT系综的分子动力学基础上，通过体积控制使系统的压力维持在给定值，或根据环境压力使系统压力发生涨落，即可实现NpT系综分子动力学。
(1)Berendsen体积标度法(fix press/berendsen命令):
根据移动方程进行标度
$\frac{dp\left(t\right)}{dt}=\frac{1}{\tau}\left[\tilde{p}-p\left(t\right)\right]$
$\lambda_{\mu v}=\delta_{\mu v}+\frac{1}{V}\left(\frac{\partial V}{\partial p_{\mu v}} \right)_{T}\left[\tilde{p}_{\mu v}-p_{\mu v}\left(t\right)\right]\frac{\triangle t}{\tau}$
$\begin{pmatrix}x_{i}\left(t+\triangle t\right) \\y_{i}\left(t+\triangle t\right) \\z_{i}\left(t+\triangle t\right) \end{pmatrix}\to \begin{pmatrix} \lambda_{xx}&\lambda_{xy}+\lambda_{yx}&\lambda_{xz}+\lambda_{zx}\\ 0& \lambda_{yy}&\lambda_{yz}+\lambda_{zy}\\ 0& 0 &\lambda_{zz} \end{pmatrix}^{\frac{1}{3}}\begin{pmatrix}x_{i}\left(t\right) \\y_{i}\left(t\right) \\z_{i}\left(t\right) \end{pmatrix}$
(2)Nose-Hoover热浴压浴法(fix npt命令):
根据恒温恒压动力学方程
$\begin{cases}\frac{d\overrightarrow{v_{i}}}{dt}=\frac{\overrightarrow{F_{i}}\left(t\right)}{m_{i}}-\xi\left(t\right)\overrightarrow{v_{i}}\left(t\right)-\zeta\left(t\right)\overrightarrow{v_{i}}\left(t\right)-\left(1+\frac{d}{g}\right)\eta\left(t\right)\overrightarrow{v_{i}}\left(t\right) \\ \frac{d\overrightarrow{r_{i}}}{dt}=\overrightarrow{v_{i}}\left(t\right)+\zeta\left(t\right)\overrightarrow{r_{i}}\left(t\right)+\eta\left(t\right)\overrightarrow{r_{i}}\left(t\right) \end{cases}$
增加额外的温度限制 $\xi (t)$ 、压强限制 $\zeta (t)$ 和 $\eta (t)$ ，则待定系数满足如下方程。
$\frac{d\xi\left(t\right)}{dt}=\frac{gk_{B}}{Q_{T}}\left[T\left(t\right)-\overrightarrow{T}\right]+\frac{Q_{p}}{Q_{T}}\eta^{2}\left(t\right)+\frac{Q_{0}}{Q_{T}}tr\left(\zeta^{\top}\left(t\right)\zeta\left(t\right)\right)-\frac{d^{2}k_{B}}{Q_{T}}\overrightarrow{T}\quad\quad Q_{T}=gk_{B}\overrightarrow{T}\tau_{T}^{2} \\ \frac{d\zeta\left(t\right)}{dt}=\frac{V\left(t\right)}{Q_{0}}\left[p\left(t\right)-\overrightarrow{p}I\right]-\frac{V\left(t\right)}{dQ_{0}}tr\left(\left[p\left(t\right)-\overrightarrow{p}I\right]I\right)-\xi\left(t\right)\zeta\left(t\right)\quad\quad\quad Q_{p}=gk_{B}\overrightarrow{T}\tau_{p}^{2} \\ \frac{d\eta\left(t\right)}{dt}=\frac{dV\left(t\right)}{Q_{p}}\left[p\left(t\right)-\overrightarrow{p}\right]+\frac{dk_{B}}{Q_{p}}T\left(t\right)-\xi\left(t\right)\eta\left(t\right) \\ \frac{dV\left(t\right)}{dt}=d\eta\left(t\right)V\left(t\right)$

NpH系综分子动力学

等焓压强控制
(1)Andersen压浴法
根据恒压动力学方程
$\begin{cases} \frac{d\overrightarrow{v_{i}}}{dt}=\frac{\overrightarrow{F_{i}}\left(t\right)}{m_{i}}-\eta\left(t\right)\overrightarrow{v_{i}}\left(t\right) \\ \frac{d\overrightarrow{r_{i}}}{dt}=\overrightarrow{v_{i}}\left(t\right)+\eta\left(t\right)\overrightarrow{r_{i}}\left(t\right) \end{cases}$
将一个虚拟质量为W，势能为 $\overrightarrow{p}V$ ,动能为 $\frac{1}{2}M(\frac{dV}{dt})^2$ 的活塞与体系相关联，则待定系数 $\eta (t)$ 满足如下方程。
$\eta\left(t\right)=\frac{1}{3}\frac{dlnV\left(t\right)}{dt}\quad\quad其中M\frac{d^{2}V}{dt^{2}}=\frac{1}{3V\left(t\right)}\left(\sum_{i}m_{i}\overrightarrow{v_{i}}^{2}\left(t\right)-\sum_{i}\overrightarrow{r_{i}}\left(t\right)\frac{\partial U}{\partial \overrightarrow{r_{i}}}\right) -\tilde{p}$
(2)Parrinello-Rahman压浴法(fix nph命令):
$记\overrightarrow{s_{i}}= \begin{pmatrix} prj_{\overrightarrow{a_{1}}}\overrightarrow{r_{i}} \\prj_{\overrightarrow{a_{2}}}\overrightarrow{r_{i}} \\prj_{\overrightarrow{a_{3}}}\overrightarrow{r_{i}} \end{pmatrix},引入晶胞度规矩阵G=\overrightarrow{h}^{\top}\overrightarrow{h}=\begin{pmatrix}\overrightarrow{a_{1}} \\ \overrightarrow{a_{2}} \\ \overrightarrow{a_{3}} \end{pmatrix}\left(\overrightarrow{a_{1}}\quad\overrightarrow{a_{2}}\quad\overrightarrow{a_{3}}\right)$
则动力学方程转化为:
$\begin{cases}M\frac{d^{2}\overrightarrow{h}}{dt^{2}}=\frac{V^{2}}{2\pi}\left[\frac{1}{V}\left(\sum_{i}m_{i}\overrightarrow{v_{i}}\overrightarrow{v_{i}}-\sum_{i}\sum_{j\neq i}\frac{1}{r_{ij}}\frac{dU\left(r_{ij}\right)}{dr_{ij}}\overrightarrow{r_{ij}}\overrightarrow{r_{ij}}\right)-p\right]\left(\overrightarrow{h}^{\top}\right)^{-1}\quad &\left(晶胞运动方程\right) \\ \frac{d^{2}\overrightarrow{s_{i}}}{dt^{2}}=-\frac{1}{2}\sum_{j\neq i}\frac{1}{m_{i}r_{ij}}\frac{dU\left(r_{ij}\right)}{dr_{ij}}\left(\overrightarrow{s_{i}}-\overrightarrow{s_{j}}\right)-G^{-1}\frac{dG}{dt}\frac{d\overrightarrow{s_{i}}}{dt} & \left(晶胞内粒子运动方程\right) \end{cases}$