一、 分子动力学简介
分子动力学（Molecular Dynamics，MD）是一套分子模拟方法，该方法主要是通过牛顿力学来模拟分子体系的运动，之后再通过统计方法得到系统的宏观性质。与密度泛函理论（DFT）相比，分子动力学的优点在于计算效率很高，但是缺点就是计算精度不咋地。增加求解系统的分子数量可以有效地提高计算精度，但是由于MD要考虑两个粒子之间的对能，因此计算成本也在指数级上升。周期性边界条件（Periodic Boundary Conditions，PBCs）可以解决这一难题。通过将模拟区域周期性地扩大，实际上一个系统只需要计算64个分子，它的一些宏观性质（如密度、温度）就能够计算得很好，但是像扩散这种大尺度的物理量依然无法计算。
本文通过DeePMD-kit:2.2.8-cuda12.0镜像中自带的LAMMPS软件，对LiCl的熔化进行深度势能分子动力学（Deep Potential Molecular Dynamics，DPMD）计算，从而通过密度变化的相变来预测LiCl的熔点。
二、分子动力学的计算原理
2.1 运动方程
在MD计算中，所有分子在相空间内都必须满足哈密顿正则方程，即：
在非相对论情况下可以对上式泰勒展开，由时刻和时刻的位置坐标可以求出时刻的位置如下：
可以看到这里的可以精确到量级。然而，随着模拟的推进，的算法由于引进了之前累积下来的误差，其误差会越来越大。可以证明，如果模拟的时长，那么的累积误差可以达到量级。
同样，我们可以写出速度的递推关系式：
可以精确到量级（速度误差不累计）。
以上两式被合称为Verlet算法，被广泛应用于MD模拟中。对于Verlet算法的改进有很多，在此处就不一一列举了，感兴趣的朋友可以前往分子动力学的积分算法深入了解。
2.2 势能面
现在我们了解了MD的运动方程，相信大家都发现了一个现实的问题：怎么算？一般情况下，我们需要计算出系统的势能面，根据势能面求导得到的值。这个势能面是以所有分子的位置为参数的函数，给定的分子的空间分布，也就可以得到唯一的势能面。
那么如何计算势能面呢？一种方法是依靠一些简单的物理图像和经验参数，自己给出一个经验公式，比如最经典的L-J势：
就是根据短程排斥和长程吸引的物理图像拟合出来的一个经验公式，常用来描述弱相互作用力，例如范德化力。这种依靠经验公式计算势能面的MD方法被称为经典分子动力学（Classical Molecular Dynamics，CMD）。CMD的计算速度很快，但是结果的精度不高。
另一种办法是依靠DFT等方法进行第一性原理计算，被称为从头算分子动力学（ab initio molecular dynamics，AIMD）。AIMD的计算量极大，耗时耗力，但是结果的精度很高。
最近的AI浪潮催生了用深度学习实现势能面的计算，即先用AIMD精确建模，之后利用深度学习拟合出一个高维的势函数。这种方法被称为深度势能分子动力学（Deep Potential Molecular Dynamics，DPMD），即本文将要使用的MD计算方法。DPMD综合了CMD和AIMD的优点，既能算得快，又能算得好，因此被广为使用。网上有人整理了关于DPMD的文章，并建立了一个网站：科学智能广场。我们可以在这里搜索已经训练好的深度势能模型，下载下来为我们使用。
三、 计算实例（LiCl的熔化）
在数据集中，我们准备了所需的所有文件。之后的教程会依次解释各类文件的内容，但是在此之前，由于我们无法修改系统文件夹中的数据，我们建议先把数据集复制到data文件夹下：
./data └── LiCl-molten-DPMD-2i5r └── v1 └── LiCl-molten ├── LiCl-molten-model.pb ├── LiCl-molten-model_new.pb ├── licl.data ├── licl.data1000 ├── licl.data700 ├── licl.data750 ├── licl.data800 ├── licl.data850 ├── licl.data900 ├── licl.data950 ├── licl.in ├── log.lammps ├── log1000.lammps ├── log700.lammps ├── log750.lammps ├── log800.lammps ├── log850.lammps ├── log900.lammps ├── log950.lammps └── log_lammps.py 3 directories, 20 files
进入LiCl-molten文件夹，开始我们的教程：
当前路径为：/data/LiCl-molten-DPMD-2i5r/v1/LiCl-molten
3.1 准备深度势能模型
首先先在科学智能广场上下载LiCl的熔化模型。我们已经为你事先准备好了，命名为LiCl-molten-model.pb。由于该模型的版本与LAMMPS软件的版本不一样，所以我们用以下命令更新模型的版本：
2024-05-31 20:33:34.053747: E external/local_xla/xla/stream_executor/cuda/cuda_dnn.cc:9261] Unable to register cuDNN factory: Attempting to register factory for plugin cuDNN when one has already been registered 2024-05-31 20:33:34.053812: E external/local_xla/xla/stream_executor/cuda/cuda_fft.cc:607] Unable to register cuFFT factory: Attempting to register factory for plugin cuFFT when one has already been registered 2024-05-31 20:33:34.119886: E external/local_xla/xla/stream_executor/cuda/cuda_blas.cc:1515] Unable to register cuBLAS factory: Attempting to register factory for plugin cuBLAS when one has already been registered 2024-05-31 20:33:34.214334: I tensorflow/core/platform/cpu_feature_guard.cc:182] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations. To enable the following instructions: SSE4.1 SSE4.2 AVX AVX2 AVX512F FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags. WARNING:tensorflow:From /opt/deepmd-kit-2.2.8/lib/python3.11/site-packages/tensorflow/python/compat/v2_compat.py:108: disable_resource_variables (from tensorflow.python.ops.variable_scope) is deprecated and will be removed in a future version. Instructions for updating: non-resource variables are not supported in the long term WARNING:root:To get the best performance, it is recommended to adjust the number of threads by setting the environment variables OMP_NUM_THREADS, TF_INTRA_OP_PARALLELISM_THREADS, and TF_INTER_OP_PARALLELISM_THREADS. See https://deepmd.rtfd.io/parallelism/ for more information. DEEPMD INFO the converted output model (2.2.8 support) is saved in LiCl-molten-model_new.pb
更新版本后的模型被命名为LiCl-molten-model_new.pb。之后我们就用这个文件作为势能面的计算输入模型。
3.2 设置输入文件licl.in
准备好模型后，我们可以用cat命令查看输入文件licl.in：
variable a loop 7 pad variable b equal $a-1 variable f equal $b*50 variable t equal 700+$f log log$t.lammps units metal boundary p p p atom_style atomic read_data licl.data replicate 2 2 2 mass 1 6.94 mass 2 35.45 pair_style deepmd ./LiCl-molten-model_new.pb out_freq 10 out_file model_devi$t.out pair_coeff * * velocity all create $t 23456789 fix 1 all npt temp $t $t 0.1 iso 0 0 0.5 timestep 0.01 thermo_style custom step temp pe ke etotal press density lx ly lz vol thermo 10 run 800 write_data licl.data$t clear next a jump licl.in
让我们来一段一段地解释这些参数的含义：
第一段的variable命令用于设置一个新的变量，这在循环执行文件的时候很有用。
a loop 7 pad定义了一个1~7的整数循环数组a，即a=[1,2,3,4,5,6,7]；
b equal $a-1定义了b=a-1，即b=[0,1,2,3,4,5,6]；
f equal $b*50定义了f=b*50，即f=[0,50,100,150,200,250,300]；
t equal 700+$f定义了t=700+f，即t=[700,750,800,850,900,950,1000]；
总而言之，第一段定义了一个数组t，之后我们将其视为系统的温度，单位为K。
第二段的log用于指定日志输出文件，这里我们把不同温度的日志放到各自对应的文件中，文件名中的$t代表温度。
第三段是初始化设置，指定了模拟中所用的单位系统为metal，边界是周期性的，原子类型为atomic。
第四段从licl.data文件中读取初始构型，其中replicte 2 2 2命令使得模拟区域扩大了2x2x2=8倍，从原来的64个原子增加到512个原子，可以模拟得更准（友情提醒：如果你想自己试着运行一遍，建议把这条命令删除，这样计算时长不会长到难以接受）。
另外，mass命令会将初始构型文件licl.data的原子质量强制更改成这里设置的值，其中第一个参数代表原子的种类（在licl.data中我们定义1代表Li原子，2代表Cl原子），第二个参数代表原子质量。
第五段的pair_style代表对能项（即势能面）从LiCl-molten-model_new.pb中读取，后面的out_freq设置了输出频率（每10步输出一次）和输出文件（文件名按t分类）。
pair_coeff代表对能项的输入参数，这里我们直接用模型自己的参数。
第六段的velocity all create $t 23456789给了所有原子一个随机速度（满足麦克斯韦分布），平均温度由t确定，后面是随机数种子，随便写。
第七段fix命令采用Nose-Hoover方法保证了npt系统的等温-等压性，具体来说，temp $t $t 0.1表示将温度控制在t附近，阻尼系数为0.1；iso 0 0 0.5表示将压强控制在0附近，阻尼系数为0.5。
timestep设置了模拟一步的时间步长为0.01 ps。
第八段的thermo_style用于设置热力学输出的格式。这里我们主要关心密度（density）。
thermo 10代表每10步输出一次热力学量。
第九段run命令设置了模拟总步长为800步，即总时间步长为8 ps。
write_data命令指定了原子位置、速度的初始构型输出按t分类，保存在相应的文件内。
第十段用于循环执行代码，即执行完一个温度下的模拟后用clear命令清除缓存，重新开始下一个温度的模拟。
以上的参量足以用于本例的模拟，更多的参量解析详见LAMMPS Commands 官方文档。
3.3 设置初始构型文件licl.data
利用cat命令查看初始构型文件licl.data：
# LAMMPS data file written by OVITO Basic 3.7.5 64 atoms 2 atom types 0.0 11.858 xlo xhi 0.0 11.858 ylo yhi 0.0 11.858 zlo zhi Masses 1 6.941 # Li 2 35.453 # Cl Atoms # atomic 1 2 4.74804 2.06958 1.06046 2 2 0.764981 2.41914 5.91261 3 2 2.61221 6.84017 12.2559 4 1 4.97269 4.0693 2.75017 5 2 4.23422 1.04616 7.12379 6 1 2.96187 6.81192 2.4587 7 1 0.927452 4.04304 4.46175 8 2 2.72948 4.78368 3.7049 9 2 6.44989 3.2905 4.20552 10 2 11.038 5.05657 5.4604 11 1 9.14871 3.8713 4.52127 12 1 7.58785 1.59778 2.79984 13 2 9.16295 1.07271 4.06425 14 1 11.3064 0.961372 3.34963 15 2 12.0987 2.50655 2.11096 16 2 0.986399 8.12845 2.68796 17 1 6.97154 8.16389 2.35447 18 2 5.32921 6.67095 3.51104 19 1 4.8602 12.0704 1.54688 20 2 2.57294 11.452 2.16311 21 1 4.72099 11.4296 7.26722 22 1 1.69269 10.0613 3.45313 23 1 5.80875 6.85847 6.22263 24 2 9.14309 8.6146 2.47043 25 2 7.07655 11.4539 1.36294 26 1 10.5019 9.42075 4.31031 27 2 11.8339 11.0282 3.89201 28 1 11.3553 6.61012 6.67926 29 2 10.9247 8.76171 6.78613 30 1 8.64752 12.7639 6.11301 31 1 12.8399 9.79091 7.4341 32 2 7.33926 11.033 6.91932 33 1 4.72085 2.75675 5.93236 34 2 7.18844 3.95799 10.3186 35 1 6.99681 2.65627 11.9433 36 1 11.3598 3.21072 7.09302 37 2 11.2323 3.32534 10.0614 38 2 9.45638 2.73118 7.16077 39 1 12.3998 5.56773 10.2205 40 1 10.3965 3.60358 12.4057 41 2 11.184 6.15175 11.7663 42 1 8.47607 2.78672 9.22135 43 2 8.45651 3.10368 1.63533 44 1 12.4799 6.44673 1.20093 45 2 1.11325 11.5997 11.3394 46 2 5.1174 4.60921 7.04953 47 1 4.60374 7.2949 11.347 48 2 2.99858 10.1566 5.99636 49 2 6.95215 7.55699 12.0361 50 1 3.13114 5.95639 8.14884 51 2 4.60896 7.18053 9.05281 52 1 12.1804 10.3914 12.5538 53 1 1.77627 12.4441 6.40526 54 2 1.11479 6.28549 8.14817 55 1 8.54528 9.39882 7.31403 56 1 2.04173 1.22117 1.40882 57 2 7.82716 7.24856 7.28153 58 1 6.75125 5.64433 8.8055 59 1 9.10628 7.66698 11.5961 60 1 12.0207 12.6303 9.86709 61 2 10.0854 9.69115 11.0475 62 2 12.3009 11.8275 7.96787 63 2 9.04836 0.874194 10.3233 64 1 8.73489 11.3448 11.4987
第一段用于初始化设置，设置了64个原子（但是由于licl.in文件内的设置，实际上要计算256个原子），2种原子类型，以及边长为11.858 的盒子。
第二段设置了原子的编号与质量，如前所述，编号1代表Li原子，编号2代表Cl原子。
第三段一共64行，每一行设置了对应编号原子的类型与初始位置。
3.4 运行代码
准备好.pb、.in和.data文件后，就可以用下面这行代码正式进行DPMD计算了：
lmp -i licl.in
由于计算时间较长（5h以上），数据集内直接给出了算好了的输出文件。（再次提醒，如果想要在本教程提供的节点内计算的话，最好把.in文件内的replicate 2 2 2语句去掉）。
我们可以从.lammps文件中提取想要的系统信息。数据集内准备了批量提取数据的.py文件，可以用cat命令查看：
import numpy as np
def calculate_mean(file_name, column_name):
with open(file_name, 'r') as file:
lines = file.readlines()
step_line = next(i for i, line in enumerate(lines) if 'Step' in line)
loop_line = next(i for i, line in enumerate(lines) if 'Loop' in line)
column_index = lines[step_line].split().index(column_name)
data = [float(line.split()[column_index]) for line in lines[step_line+1:loop_line]]
mean = sum(data[21:]) / len(data[21:]) #skip the first 200 md steps
return mean
temps=[700,750,800,850,900,950,1000]
for temp in temps:
mean_density = calculate_mean(f'log{temp}.lammps', 'Density')
mean_Lx = calculate_mean(f'log{temp}.lammps', 'Lx')
mean_Eng = calculate_mean(f'log{temp}.lammps','TotEng')/512
print(round(mean_density,3), round(mean_Lx,15), round(mean_Eng,15))
其中我们忽略了每次模拟的钱200步数据，因为那时候系统还没有平衡。运行这个.py文件就可以看到模拟结果了：
1.624 22.305886216666668 -3.610760804036458 1.6 22.417404916666673 -3.590539934895835 1.577 22.523178783333332 -3.572294033203126 1.555 22.62926185 -3.554371552734374 1.54 22.70282385 -3.539264615885416 1.513 22.835978933333326 -3.517140921223959 1.498 22.9154875 -3.501806022135416
画成折线图长这样：
可以发现两幅图在900K-950K的地方都有突变，且在突变两侧的斜率不同，说明在这个温度范围内存在一级相变点，即熔点。实验测量得LiCl的熔点在870K左右，符合得较好。
系统大概在2ps左右趋于平衡，因此总模拟时长至少要5ps；而由于模拟需要的现实时间较长（每一步约3s），所以为了让总模拟时间可以接受（8h以内），可以将每一步的计算时间dt设置为0.01ps。 如果改变dt到0.1ps，会发生粒子损失的错误，不予采用。