樊哲勇《分子动力学模拟》python实例 |5.压力控制算法mtkk的编程实现

樊哲勇老师最近出版的书籍《分子动力学模拟》中，使用C++语言进行分子模拟程序的实现。未了更好地理解分子动力学模拟，AI4S1O1-材料小队学习志愿者决定针对樊老师书籍中的内容，进行python语言的程序实现以及注释说明。

本节内容对应樊老师书中的第四章

樊哲勇《分子动力学模拟》python实例 | 概述

本章目标

1. 了解分子动力学模拟中的压力控制常用算法
2. 进行python编码实现
3. 进行不同算法之间的对比
4. 控压知识扩展

一、压力控制算法

1.1定义

压力控制算法，也称为压力耦合算法或恒压算法，是在分子动力学（MD）模拟中用来调节系统压力的方法。这些算法通过调节系统的体积或者施加外部压力，使得系统在模拟过程中保持在预设的压力状态下。

在MD模拟中，压力可以通过以下几种方式进行控制：

1. 体积调节（Volume Adjustment）：通过改变模拟盒的体积来调节系统的压力。
2. 应力张量（Stress Tensor）：通过施加外部应力，调节系统内部的力平衡。 常用的压力控制算法

1.2 常用的压力控制算法

Berendsen 控压算法：该算法通过逐步调整系统的体积，使得系统压力逐渐趋向目标压力。简单，高效，适合初步平衡；不能保证正确的热力学统计，越接近目标压力调节力度减弱。

SCR（stochastic cell rescaling) 控压算法：在模拟过程中随机地对模拟盒的体积进行重新缩放，从而调整系统的压力。快速平衡，保持了动力学一致，易于实现；随机性引入了不稳定性。

MTTK（Martyna-Tuckerman-Tobias-Klein）控压算法：算法基于扩展系综方法，通过引入额外的动态变量（如体积的变化率和应力张量）来调节系统的压力。精确控制压力，适用于复杂系统；计算复杂度高，趋于平衡的过程表现较上两个差。

PRB（Parrinello-Rahman Barostat）恒压算法:算法引入一个变换矩阵（通常表示为 ℎ）来描述模拟盒的形状和大小变化，通过将模拟盒的形状和体积作为动态变量，使系统在恒定外界压力下进行演化。能处理各向异性应变，精确热力学统计，适用于复杂系统；实现复杂，计算开销大，参数调节困难。

Nose-Hoover 链恒压算法：算法通过引入虚拟变量（称为热浴变量和体积变量），使系统能够在恒温恒压（NPT）系综下演化，适用于高精度模拟和复杂系统。缺点是实现复杂，调节参数困难。

因为篇幅有限，本章节主要示范Berendsen 控压算法，*SCR控压算法，MTTK控压算法

二、控温算法代码实现

Berendsen 控温算法

Berendsen 控压算法通过将系统体积按指数方式调节，使得系统的压力逐渐逼近目标压力。具体来说，算法根据当前压力和目标压力之间的差异，调整系统的体积。

考虑各向同性时，体系的压强变化率正比于实际压强 p 与目标压强p0 的差值： $$\frac{\mathrm{d}p}{\mathrm{d}t}=\frac{1}{{\tau }_p}(p_0-p).$$ 当体系的实际压强小于目标压强时， dp/dt 为正，体系压强有增加的趋势；当体系的 实际压强大于目标压强时， dp/dt 为负，体系压强有减小的趋势。参数 τp 具有时间的量纲， 它反映了压强变化的快慢，可被称为控压算法的弛豫时间。引入等温压缩率 $$\beta =-\frac{1}{V}\frac{\partial V}{\partial p}$$ 可进一步的表达，如同对体积做如下标度变换： $$V\left(1-\frac{\beta \Delta t}{{\tau }_p}(p_0-p)\right)$$ 相当于改变坐标： $${\mathbf{r}}_i{\left(1-\frac{\beta \Delta t}{{\tau }_p}(p_0-p)\right)}^{1/3}$$ 由于每一步变化不大，坐标变换等同如下： $${\mathbf{r}}_i\left(1-\frac{\beta \Delta t}{3{\tau }_p}(p_0-p)\right)$$ 同时如果盒子为立方体，盒子边长L改变： $$L\left(1-\frac{\beta \Delta t}{3{\tau }_p}(p_0-p)\right)$$ 推广到一般的情况，体系的盒子矩阵 hαβ 和原子坐标 ri 都按照一个形变矩阵 µαβ 变换： $${\mu }_{\alpha \beta }={\delta }_{\alpha \beta }-\frac{{\beta }_{\alpha \beta }\Delta t}{3{\tau }_p}(p_{0\alpha \beta }-p_{\alpha \beta }).$$ 其中， ${\mu }_{\alpha \beta }$是等温压缩率张量。这代表最一般的 6 自由度控压。

SCR算法

SCR 算法的核心思想是在模拟过程中随机地对模拟盒的体积进行重新缩放，从而调整系统的压力。与传统的压力控制算法相比，SCR 采用了一种随机过程来控制压力，这种随机性有助于系统更快地达到平衡。与 Berendsen 控压算法相比， SRC 控压算法给形变矩阵增加了一个随机项 $${\mu }_{\alpha \beta }^{\text{stochastic}}=\sqrt{\frac{1}{D_{\mathrm{couple}}}}\sqrt{\frac{{\beta }_{\alpha \beta }\Delta t}{3{\tau }_p}\frac{2k_{\mathrm{B}}{T}_0}{V}}{R}_{\alpha \beta }$$ 其中， ${\mu }_{\alpha \beta }$是均值为零、方差为一的正态分布的随机数， $D_{\mathrm{couple}}$ 是一个与盒子自由度耦合 有关的因子。当采用各向同性控压时，该因子为 3，其他情况下为 1。

MTTK控压算法

其算法公式较为复杂，请参考樊老师的书籍第七章P142

下面的代码实现了一个分子动力学模拟，通过Nosé-Hoover链（NHC）热浴控制粒子和盒子的温度，并计算粒子的位置、速度、势能和压力变化。具体的计算思路如下：

1. 初始化参数和变量：

   • 初始化系统常数、盒子参数、粒子参数和NHC热浴参数。
   • 定义时间步长 dt 和半时间步长 dt2。
   • 初始化粒子位置 x 和速度 v，并根据初始位置和盒子长度计算力 f。

2. 定义辅助函数：

   • sinh_factor(x2)：计算双曲正弦函数的近似值。
   • find_force(x, h)：根据粒子的位置和盒子长度计算粒子所受的力。
   • find_potential(x, h)：计算粒子在给定位置处的势能。
   • find_pressure(ek2, x, h)：根据粒子的动能、位置和盒子长度计算压力。
   • nhc(...)：实现Nosé-Hoover链热浴的动力学计算，包括更新热浴的速度、位置和计算缩放因子。

3. 时间演化循环：

   • 进行多步时间演化模拟，每步更新系统的状态（位置、速度、盒子长度等）。
   • 每隔一定步数输出系统的状态数据到文件 data.txt，包括粒子的位置、盒子长度、总能量和速度。

4. 时间演化步骤：

   • 在每个时间步内，分别执行多种演化操作：
     • 更新盒子的动量。
     • 更新粒子的速度。
     • 更新盒子的动量。
     • 更新粒子的速度和位置。
     • 更新盒子长度。
     • 更新力。
     • 更新粒子的速度。
     • 更新盒子的动量。

5. 能量守恒检查：

   • 每个采样间隔计算系统的总能量，并输出到文件中，用于检查能量是否守恒。

三、不同算法之间的对比

使用了樊哲勇老师等人编写的GPUMD软件，可以实现不同算法之间的对比。具体的输入脚本，请参考：

`potential   Si_Tersoff_1989.txt
velocity    300
#ensemble    npt_scr 300 300 100 0 50 1000
#ensemble    npt_ber 300 300 100 0 50 1000
ensemble    npt_mttk temp 300 300 iso 0 0
time_step   1

dump_thermo 1

run         10000`

运行结果可视化后如下：

在 Berendsen 控 压算法中，压强和体积都以指数函数的方式趋于平衡值，达到平衡后便几乎没有涨落。 SCR控压算法在趋于平衡的过程与 Berendsen 控压算法有类似的表现，但在体系达到平衡后压强和体积都有合理的涨落，遵循 NpT 系综的统计分布。在 MTTK 控压算法中，压强和体积在体系趋于平衡的过程中都有较大的振荡。虽然 MTTK 和 SCR 算法给出的压强和盒子长度涨落的局部细节不一样，但是它们确实都给出正确的 NpT 系综。

四、控压知识扩展全自由度控压。

各向异性控压

上述方法是各向同性控压，也可以进行三轴独立控压和全自由度控压。 使用GPUMD软件实现这一过程： #三轴独立控压 potential Si_Tersoff_1989.txt velocity 300

#ensemble    npt_ber 300 300 100 10 0 0 50 50 50 1000
ensemble    npt_scr 300 300 100 10 0 0 50 50 50 1000
#ensemble    npt_mttk temp 300 300 x 10 10 y 0 0 z 0 0
time_step   1
dump_thermo 2
run         20000

#全自由度控压
potential   Si_Tersoff_1989.txt
velocity    300

#ensemble    npt_ber 300 300 100 0 0 0 0 0 10 50 50 50 50 50 50 1000
ensemble    npt_scr 300 300 100 0 0 0 0 0 10 50 50 50 50 50 50 1000
#ensemble    npt_mttk temp 300 300 x 0 0 y 0 0 z 0 0 yz 0 0 xz 0 0 xy 10 10
time_step   1
dump_thermo 2
dump_exyz   20000
run         20000

结果如下：

选择合适的控压算法

具体算法选择建议：

1. Berendsen Barostat：适用于快速达到平衡、计算效率高、适用于一般的分子动力学模拟，不需要严格的统计特性。
2. Parrinello-Rahman Barostat：适用于需要精确控制压力和处理各向异性系统的模拟，特别是在应力和应变关系的研究中。
3. MTTK控压算法：适用于高精度和长时间模拟，适合复杂生物分子和材料系统，能够提供高稳定性和准确性。
4. Stochastic Cell Rescaling (SCR)：适用于一般的分子动力学模拟，实现简单，但不适用于高精度模拟。
5. Andersen Barostat：适用于需要较好压力控制的中等精度模拟，实现复杂，计算效率较低。

实际操作建议：

1. 测试和调试：在选择控压算法时，可以先进行短时间的测试模拟，观察系统在不同算法下的平衡速度和压力波动情况。
2. 结合使用：在某些情况下，可以结合使用不同的控压算法。例如，可以先使用Berendsen Barostat快速达到平衡，然后切换到MTTK控压算法进行精确的长时间模拟。
3. 参考文献：查阅相关领域的文献，了解在类似研究中常用的控压算法，以便选择适合的算法。