分子动力学入门 | DPMD+LAMMPS运算实例 | Bohrium-玻尔科研空间站

分子动力学（Molecular Dynamics，MD）是一种通过模拟分子系统的运动来研究物质性质的计算方法，在材料、化学、生物医药等领域有着广泛的应用。

在模拟分子运动时，我们需要考虑其运动方程。目前一般仍然是用经典的分子动力学，即在相空间内描述，认为分子运动满足哈密顿正则方程。我们应用verlet算法，分别对 $t + Δ t$ 和 $t - Δ t$ 时刻的位置泰勒展开，可以得到 $r (t + Δ t) + r (t - Δ t) = 2 r (t) + \frac{f ( t )}{m} Δ t^{2} + O (Δ t^{4})$ 这时 $r (t + Δ t)$ 的误差仅为 $Δ t^{4}$ 量级，计算精度高。但该算法对描述粒子速度不够准确。粒子速度没有真正进入粒子的Verlet运动方程积分公式之中。算法所需初始数据（位置、速度、力）需额外算法求出。作为一个改进，Velocity-Verlet算法给出了更好的粒子速度描述，算法所需初始数据可全部得到，因此得到较广泛的应用。Velocity Verlet速度公式为 $v (t + δ t) = v (t) + \frac{f ( t + δ t ) + f ( t )}{2 m} ∆ t$ 算出速度等微观量之后，我们就可以得到系统的宏观量。比如温度可以通过 $E_{k} = \frac{3}{2} (N - 1) k_{B} T$ 来计算（N-1）是因为受到质心系约束。我们进行分子动力学模拟时也需要选择系综，比如NVT，NVE，NPT系综，在模拟过程中会进行系综调节，系综调节主要是指在进行分子动力学计算过程中，对温度和压力参数的调节，包括调温技术和调压技术。

为了算各个分子的运动，我们需要描述分子间的相互作用,这通过势能面来描述.描述一个系统在给定参数下的势能，这个系统通常是原子的集合，而这些参数通常是原子的坐标.如果获得了精确的势能面,我们就可以很好的描述材料的性质.

描述势能面的方法一是直接自己给定一个经验的解析的表达式,计算效率高但结果不精确.方法二是用第一性原理,即密度泛函理论来计算,这需要解薛定谔方程,需要很大的运算量,也比较精确.

现在比较火热的新方法是用深度学习的方法,我们先用第一性原理精确建模,再用深度学习方法拟合高维的势函数.这种方法可以实现计算速度和精度的统一.

关于深度学习方法是如何构造势函数的，我们给出如下说明。考虑体系的原子坐标为 $R = {r_{1}, r_{2}, \dots, r_{N}}$ ，我们一般可以通过第一性原理精确的计算得到体系能量 $E = E (R)$ ，从能量进而可以算出第i个原子的受力 $F_{i}$ ，和体系的virial张量 $Ξ_{α β} = - \frac{\partial E}{\partial h _{γ α}} h_{γ β}$ 。利用这三个量作为标签可以得到损失函数 $L = \frac{1}{∣ B ∣ N ^{2}} {k \in B \sum (p_{e} ∥ E (R^{(k)}) - E^{w} (R^{(k)}) ∥^{2} + p_{f} ∥ F (R^{(k)}) - F^{w} (R^{(k)}) ∥^{2} + p_{v} ∥ Ξ (R^{(k)}) - Ξ^{w} (R^{(k)}) ∥^{2})}$ 其中 $B$ 为批次大小，三项前各自有个系数 $p$ 调整对损失函数的贡献。我们只要最小化损失函数就完成了训练。

在物理上考虑，能量应该是广延的，即可以写成对每个原子的分解式。 $E^{w} = i = 1 \sum N E_{i}^{w} = i = 1 \sum N E_{i}^{w} (R_{i})$ 从而我们只需要对每个 $E_{i}$ 分别建模。其次，体系的能量有对称性要求，包括平移，旋转，交换对称性，这体现为，能量有约束 $E (R) = E (U R)$ U代表某种变换。但一般的神经网络是不会满足这样的对称性的，所以我们可以引入描述子： $D (R_{i}) = D (U R_{i}), \forall U \in U$ 其为对原子坐标的映射，满足所有的对称性要求，再将建模的能量函数写成关于描述子的函数，这样能量就自然的满足了对称性。描述子和神经网络的具体构造可以参考王涵.机器学习原子间相互作用建模[J].计算数学, 2021, 43(3):18.DOI:10.12286/jssx.j2021-0833.

我们下面将给出一个使用deepMD的运行实例.deepMD的各个包的作用如下 alt

我们将以ZrO2材料为例,计算其在不同温度下的密度从而获得物质的不同相.并且我们将直接使用AIS Square中的现成模型：models/215/ZrO2-PBEsol-model (aissquare.com)，经测试在dpmd:2.2.8-cuda12.0这个镜像下可以兼容运行。直接用现成模型相当于就不需要自己从头训练深度势能模型,而是可以直接应用于下游的分子动力学模拟.我们将使用dpmd+lammps来进行分子动力学模拟

设置ZrO2的结构如下,此文件为ZrO2.data

# LAMMPS data file 
36 atoms
2 atom types
0.0 4.957159 xlo xhi
0.0 5.916337 ylo yhi
0.0 16.316733 zlo zhi

Masses

1 91.224  # Zr
2 15.999  # O

Atoms  # charge

1 1 0.0 0.012204525	2.046159235	1.371470359
2 1 0.0 0	0.100962291	4.07918325
3 1 0.0 4.944954475	2.046159235	6.786896141
4 1 0.0 0	5.815374709	12.23754975
5 1 0.0 0.012204525	3.870177765	9.529836859
6 1 0.0 4.944954475	3.870177765	14.94526264
7 1 0.0 2.490784025	0.912009265	14.94526264
8 1 0.0 2.466374975	0.912009265	9.529836859
9 1 0.0 2.4785795	2.857206209	12.23754975
10 1 0.0 2.490784025	5.004327735	6.786896141
11 1 0.0 2.466374975	5.004327735	1.371470359
12 1 0.0 2.4785795	3.059130791	4.07918325
13 2 0.0 1.005227573	3.720222287	0.493059038
14 2 0.0 1.005227573	2.196114713	8.651425538
15 2 0.0 1.391970247	1.144870373	11.35536926
16 2 0.0 1.391970247	4.771466627	3.197002764
17 2 0.0 1.037518507	0.399950298	5.917246906
18 2 0.0 1.037518507	5.516386702	14.07561341
19 2 0.0 1.441060993	3.358118798	5.917246906
20 2 0.0 1.441060993	2.558218202	14.07561341
21 2 0.0 1.086609253	1.813298127	3.197002764
22 2 0.0 1.086609253	4.103038873	11.35536926
23 2 0.0 1.473351927	0.762053787	0.493059038
24 2 0.0 1.473351927	5.154283213	8.651425538
25 2 0.0 3.483807073	0.762053787	7.665307462
26 2 0.0 3.483807073	5.154283213	15.82367396
27 2 0.0 3.870549747	1.813298127	4.961363736
28 2 0.0 3.870549747	4.103038873	13.11973024
29 2 0.0 3.516098007	2.558218202	10.39948609
30 2 0.0 3.516098007	3.358118798	2.241119594
31 2 0.0 3.919640493	0.399950298	2.241119594
32 2 0.0 3.919640493	5.516386702	10.39948609
33 2 0.0 3.565188753	1.144870373	13.11973024
34 2 0.0 3.565188753	4.771466627	4.961363736
35 2 0.0 3.951931427	2.196114713	15.82367396
36 2 0.0 3.951931427	3.720222287	7.665307462

之后我们考虑设计lammps的输入文件.设计如下.各行的意思在注释中给出.考虑到会有相变,相变时系统的压强和温度是不变的.那么应该是选npt系综。

variable t equal 2000

log             log_$t.lammps #指定结果的输出文件

units           metal #选系统的单位.对于metal，时间的单位是皮秒（ps），长度的单位是埃（Å），质量的单位是原子质量单位（amu），能量的单位是电子伏特（eV），温度的单位是开尔文（K），压力的单位是巴（bar），速度的单位是埃/皮秒（Å/ps）。
boundary        p p p #用于设置模拟的边界条件。在本教程情况下，我们在 x, y, z 三个方向上均使用周期性（periodic）边界条件。p 表示周期性，f (fixed) 表示固定边界条件。
atom_style      charge #用于设置原子类型和属性。在这个示例中，我们使用带电原子模型，因此使用了 charge 类型。

read_data	ZrO2.data #用于读取数据文件ZrO2.data，同时也是该案例中 MD 模拟所使用的初始构型。
mass 		1 91.224
mass		2 15.999 #设置两种原子的质量
group		Zr  type 1
group		O  type 2 #用于根据原子类型创建两个组，分别包含元素类型为 1（Zr）和元素类型为 2（O）的所有原子。

pair_style	deepmd ./ZrO2-PBEsol.pb #设置模型为ZrO2-PBEsol.pb
pair_coeff  	* *	#自动根据模型设置

velocity        all create $t 23456789 #所有原子初始速度随机

fix             1 all npt temp $t _$t 0.1 iso 0 0 0.5 #设置为npt系综,温度为t,并使用Nose-Hoover 算法进行温度控制。0.5 是温度阻尼系数。iso代表对压强的控制
timestep        0.001 #时间步长为0.001ps

thermo_style    custom step temp pe ke etotal press density lx ly lz vol #设置输出格式
thermo          100 #每隔100步输出信息

dump		1 all custom 100 ZrO2_$t.dump id type x y z #用于在模拟过程中输出原子的坐标信息。在这个示例中，使用 `dump 1 all` 表示输出模拟系统中的所有原子，使用 `custom` 表示使用自定义的输出格式，使用 `100` 表示输出频率，使用 `ZrO2_$t.dump` 表示输出文件的名称，使用 `id type x y z` 表示输出的原子信息，包括原子的 ID、类型以及在 x、y、z 方向上的坐标。

run             5000 #运行5000步,即5ps

由于bohrium平台的问题，上面fix一行需要把_$t前的_删掉才是正确的输入文件

我们用脚本来获得不同温度下的密度

#!/bin/bash

for ((t=2400; t<=3400; t+=100))
do
    sed -i "1s/.*/variable t equal $t/" "ZrO2.in"
    lmp -i ZrO2.in

    echo $t >> res.txt
    awk 'NR==176 {print $7}' log_$t.lammps  >> res.txt
done

最终画出来的图像如下

alt

这跟模型的原文献里给出的相图有些类似

J.-Y. Zhang, G. Huynh, F.-Z. Dai, A. Tristan, S. H. Zhang, S. Ogata, D. Rodney, A deep-neural network potential to study transformation-induced plasticity in zirconia, J. Eur. Ceram. Soc. 44 (6) (2024) 4243-4254. (doi: 10.1016/j.jeurceramsoc.2024.01.007)

alt

网上搜到的熔点的标准值约为3000K，在画的图像中，3000K时也看到了明显的密度降低，这证明了我们的计算的合理性。

参考:

【张与之：DeePMD原理介绍与上机实践】https://www.bilibili.com/video/BV1Nh4y1t73L/share_source=copy_web&vd_source=c2bb13397c89aa4ca527531eb8f1dd8c

从 DFT 到 MD｜超详细「深度势能」材料计算上手指南 - Bohrium (dp.tech)