计算物理（三）分子动力学方法

本篇文章是西安交大《计算物理与程序设计》课程分子动力学部分的笔记和作业。

分子动力学

Verlet

Velocity Verlet 算法是一个二阶算法，能够根据前一时刻的位置、速度和作用力，计算出下一时刻的位置、速度和作用力。即对于方程

$$\begin{cases} \dot{\bm{r}} = \bm{v} \\ m \ddot{\bm{r}} = \bm{F} \end{cases}$$

Velocity Verlet 算法给出

$$\begin{cases} \boldsymbol{r}(t+\Delta t)=\boldsymbol{r}(t)+ \boldsymbol{v}(t) \Delta t+ \frac{1}{2} \boldsymbol{a}(t) (\Delta t)^2 \\ \boldsymbol{v}(t+\Delta t)=\boldsymbol{v}(t)+ \frac{1}{2} [\boldsymbol{a}(t)+\boldsymbol{a}(t+\Delta t)]\Delta t \end{cases}$$

式中 $\bm{a}(t)=\bm{F}[\bm{r}(t)] / m$ 是根据 $t$ 时刻系统位形确定的粒子加速度。Velocity Verlet 算法的精度只有 $O(\mathrm{d}t^3)$ ，但相比四阶RK算法其优点在于 Verlet 算法保辛形式 (symplectic integrator)，满足能量守恒，故分子动力学模拟中常常使用 Velocity Verlet 算法。

Velocity Verlet 算法的 MATLAB 代码为

function [r,v] = verlet(force,mass,r0,v0,t)
% verlet.m
dim = length(r0);
N = length(t);
dt = t(2) - t(1);
r = zeros(dim,N); v = zeros(dim,N); % initialization
r(:,1) = r0; v(:,1) = v0; % initial value

for i = 1:N-1
    r(:,i+1) = r(:,i) + v(:,i)*dt + 0.5*force(r(:,i))/mass*dt*dt;
    v(:,i+1) = v(:,i) + 0.5*(force(r(:,i))+force(r(:,i+1)))/mass*dt;
end
end

质心速度调整

由于速度是围绕零点随机分布的，系统的总动量接近零，但不会精确地为零。为了消除系统的这种空间漂移，必须调整质心速度为零，一般把原子的速度变换到质心坐标系中

$$\boldsymbol{v}_C=\frac{\sum{m_i\boldsymbol{v}_i}}{\sum{m_i}}\,\,, \boldsymbol{v}_i \to \boldsymbol{v}_i-\boldsymbol{v}_C$$

周期性边界条件

为了保证模拟区域的粒子数不变要采用周期性的边界条件。当粒子运动时，如果出了模拟区域，施加周期性边界条件使粒子回到系统中。

最小镜像法则

最小镜像法则保证每个原子计算作用力时，仅计算离它最近的原子，而非多次计算周期性边界带来的无限原子。

• 若 $r_i-r_j>{L/2}$，则取 $r_{ij}=r_i-\left(r_j+L \right)$；
• 若 $r_i-r_j<{-L/2}$，则取 $r_{ij}=r_i-\left(r_j-L \right)$；
• 否则取 $r_{ij}=r_i-r_j$ 。

速度标度法

经典分子动力学模拟中常常使用微正则系综(microcanonical ensemble)，微正则系综具有确定的粒子数 $N$、体积 $V$ 和能量 $E$ 。然而现实实验环境中，相比系统的总能量 $E$ ，温度 $T$ 通常更容易维持恒定，只需设置一大热源与系统接触即可，这一过程称为控温(thermostat)。因此统计物理中更常用的是粒子数 $N$ 、体积 $V$ 和温度 $T$ 都相同的热力学体系组成的正则系综(canonical ensemble)。

由于统计物理系统的特性，精确地控制初态是困难的。模拟过程中人为设置的系统初态不会具有正则系综所要求的温度 $T_\text{bath}$ ，并且系统的状态不是一个平衡态。在系统趋衡过程中需要不断调节速度，使系统的温度达到给定系统的温度，即所谓速度标度法。

对速度的重新标度通过标度因子 $\beta$ 实现^[1]

$$\beta=\left(\frac{\overline{E_\text{k}}}{E_\text{k}}\right)^{\frac{1}{2}}=\left[\frac{3(N-1) k_{\text{B}} T}{ \sum_{i=1}^N m v_i^2}\right]^{\frac{1}{2}}$$

无量纲化后的标度因子为^[2]

$$\beta^* = \left[\frac{(N-1) T^*}{16 \sum_{i=1}^N v_i^{* 2}}\right]^{\frac{1}{2}}$$

Nosé-Hoover控温法

速度标度法并不具有物理意义，速度的突变无法描述热源与体系之间的相互作用及其带来的物理化学效应。Nosé和Hoover提出^[3][4]可以通过扩展系统哈密顿量，引入粒子与热源的相互作用参数 $\zeta$ 作为系统额外的自由度，从而达到描述和调控体系温度的目的，即所谓Nosé-Hoover控温法。

Nosé-Hoover控温法给出的扩展系统哈密顿量为

$$\mathcal{H} \left( P,R,p_s,s \right) =\sum_i{\frac{\boldsymbol{p}_{i}^{2}}{2ms^2}}+\frac{1}{2}\sum_{i\ne j}{U}\left( \boldsymbol{r}_i-\boldsymbol{r}_j \right) +\frac{p_{s}^{2}}{2Q}+gk_{\mathrm{B}}T_{\text{bath}}\ln s$$

其中 $s$ 和 $p_s$ 为引入的额外坐标和共轭动量，$Q$ 为控温质量，$T_\text{bath}$ 为热源温度，$g$ 为系统独立自由度，对于三维 $N$ 粒子系统质心静止故有 $g = 3(N-1)$。由哈密顿正则方程，可得运动方程为

$$\begin{gathered} \frac{\mathrm{d} \boldsymbol{r}_i}{\mathrm{d} t}=\frac{\partial \mathcal{H}}{\partial \boldsymbol{p}_i}=\frac{\boldsymbol{p}_i}{m s^2} \\ \frac{\mathrm{d} s}{\mathrm{d} t}=\frac{\partial \mathcal{H}}{\partial p_s}=\frac{p_s}{Q} \\ \frac{\mathrm{d} \boldsymbol{p}_i}{\mathrm{d} t}=-\frac{\partial \mathcal{H}}{\partial \boldsymbol{r}_i}=-\nabla_i U(R)=-\sum_{i<j} \nabla_i U\left(\boldsymbol{r}_i-\boldsymbol{r}_j\right) \\ \frac{\mathrm{d} p_s}{\mathrm{d} t}=-\frac{\partial \mathcal{H}}{\partial s}=\frac{1}{s}\left(\sum m v_i^2 / s^2-g k_{\mathrm{B}} T_{\text{bath}}\right) \end{gathered}$$

引入耦合参量 $\zeta =s{p_s / Q}$ 描述系统与热源之间的相互作用，运动方程化简为^[5]

$$\begin{gathered} m\ddot{\boldsymbol{r}}=\boldsymbol{F}-\zeta m\boldsymbol{v}\\ \dot{\zeta}=\frac{1}{Q}\left[ \sum_{i=1}^N{\frac{1}{2}mv_{i}^{2}}-\frac{3\left( N-1 \right)}{2}k_{\mathrm{B}}T_{\text{bath}} \right] \end{gathered}$$

我们希望当系统温度 $T_{\text{sys}}=\dfrac{\sum_{i=1}^N{\frac{1}{2}m_iv_{i}^{2}}}{\frac{3}{2}\left( N-1 \right) k_{\text{B}}}\implies T_{\text{sys}}^{*}=\dfrac{32K^*}{N-1}$ 趋于热源温度 $T_\text{bath}$ 时，运动方程中的热源导致的耦合项 $\zeta m\boldsymbol{v}$ 能够消失，即 $T_{\text{sys}}\to T_{\text{bath}}$ 时有 $\zeta \to 0$，粒子的运动完全由粒子内部的分子势决定。Berendsen 提出 $\zeta$ 的初值可取^[6]

$$\zeta _0=1-\frac{T_{\text{bath}}}{T}$$

尔后按照运动方程演化即可。根据 Verlet 算法可将运动方程差分为

$$\begin{gathered} \boldsymbol{r}_i(t+\delta t)=\boldsymbol{r}_i(t)+\boldsymbol{v}_i(t)\delta t+\left[ \frac{\boldsymbol{F}_i(t)}{m_i}-\zeta (t)\boldsymbol{v}_i(t) \right] \frac{\delta t^2}{2}\\ \boldsymbol{v}_i(t+\delta t/2)=\boldsymbol{v}_i(t)+\frac{\delta t}{2}\left[ \frac{\boldsymbol{F}_i(t)}{m_i}-\zeta (t)\boldsymbol{v}_i(t) \right] \\ \boldsymbol{F}_i(t+\delta t)=\boldsymbol{F}_i\left( \boldsymbol{r}_i(t+\delta t) \right) \\ \zeta (t+\delta t)=\zeta \left( t \right) +\frac{\delta t}{2Q}\left[ \sum_{i=1}^N{m_i}\frac{\boldsymbol{v}_i(t)^2}{2}+\sum_{i=1}^N{m_i}\frac{\boldsymbol{v}_i(t+\delta t/2)^2}{2}-3\left( N-1 \right) k_BT \right] \\ \boldsymbol{v}_i(t+\delta t)=\frac{\left[ \boldsymbol{v}_i(t+\delta t/2)+\frac{\delta t}{2}\frac{\boldsymbol{F}_i(t+\delta t)}{m_i} \right]}{1+\frac{\delta t}{2}\zeta (t+\delta t)} \end{gathered}$$

运用分子动力学方法即可计算。

---

1. Woodcock, L. V. Isothermal Molecular Dynamics Calculations for Liquid Salts. Chem. Phys. Lett. 1971, 10, 257−261. ↩︎

2. 齐楠, 王圣军. (2022). 分子动力学方法中速度标度因子的推导. 大学物理, 41(4), 61. doi:10.16854/j.cnki.1000-0712.210204 ↩︎

3. Nosé, S. (1984). A unified formulation of the constant temperature molecular dynamics methods. Journal of Chemical Physics, 81(1), 511–519. ↩︎

4. Hoover, W. G. (1985). Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A, 31, 1695–1697. ↩︎

5. Thijssen, J. (2007). Computational Physics (2nd ed.). Cambridge: Cambridge University Press. ↩︎

6. Berendsen, H. J., Postma, J. V., Van Gunsteren, W. F., DiNola, A. R. H. J., & Haak, J. R. (1984). Molecular dynamics with coupling to an external bath. The Journal of chemical physics, 81(8), 3684-3690. ↩︎