3.2.1 做init需要准备的文件2024.04.04

DPGEN init_bulk模块运行过程大致解释，DPGEN软件主要是起到调度ABACUS+dpdata+DPDiscapter进行相关的计算，运行过程中包含四部分的内容。1.优化结构获得稳定构型、2.对稳定构型进行微扰获得更多有用数据构型、3.对构型进行DFT的MD计算获得力和能量的结果、4.收集计算结果中的力和能量信息并转换成DeePMD能识别的数据集。

从四个过程可以看出需要的文件包括1.ABACUS做结构优化需要的INPUT，做MD需要的INPUT、STRU、*.upf、*.orb文件，以及DPGEN作业安排init.json和machine.json

1. 做结构优化获得平衡结构的稳定构型需要的ABACUS的INPUT文件。

做基态优化就需要调用ABACUS来进行计算，因此需要准备STRU文件，INPUT文件，势函数文件，由于前面已经讲过相关文件的准备。此处只给出INPUT文件，需要注意，做结构优化和后续微扰最好去掉晶体对称性，进而获得足够多的构型。我的INPUT文件以INPUT_cell_realx来命名，内容如下（不一定是合适的，因为我DFT计算经验不足，只供参考，欢迎给出改进意见）

另外也需要DPGEN运行的输入文件init.json和机器配置文件machine.json，见下一节设置。

1. ABACUS的INPUT文件示例

INPUT_PARAMETERS
#Parameters (1.General)
calculation         cell-relax

#Parameters (2.Iteration)
ecutwfc             100
scf_nmax            100
#Parameters (3.Basis)
basis_type          lcao
symmetry            0
kspacing            0.08

#Parameters (4.Smearing)
smearing_method     gauss
smearing_sigma      0.002

#Parameters (5.Mixing)
mixing_type         pulay
mixing_beta         0.7

cal_stress          1
cal_force           1

2. 需要注意的地方；

calculation         cell-relax#对晶胞进行弛豫优化，获得稳定构型。

关于稳定构性，我碰到过一些自认为不合理的事儿，因此需要补充一下我的看法。第一次准备初试训练集时，以为不用平衡结构来做md也是可以的，因此，跳过了计算耗时的结构优化部分。现在仔细考虑，感觉不应该跳过初始稳定结构的优化。原因如下，后续的探索过程基本都是在平衡态附近进行，如果初始构型包含的平衡态比较少的话在探索计时阶段容易获得较大的力得偏差，此处可能不利于我们快速获得一个可迭代的势函数（当然经过多次的迭代，这个问题应该是可以抵消的，只是我们可能需要额外处理一些情况，例如LAMMPS刚开始跑就崩溃的情况，需要我们降低步数来进行前期迭代）。

symmetry            0#关闭对称性

关于对称性分析的思考，刚开始使用软件时，由于没有经验，开启了对称性功能，以为这样可以加快计算速度。在对DPGEN流程大致了解之后，感觉应该关闭对称性。主要是因为DPGEN对平衡结构进行后续的微扰时会破坏晶体的对称性，这时开启对称性反而是不对的做法。

kspacing            0.08#使用kspacing指定K点密度，不用KPT文件指定

关于K点的思考，刚开始计算时使用的是KPT文件，对所有构型都使用5 5 5 0 0 0 的K点，后续计算结果出现能量偏差恒定常数可能跟这里相关（但目前没确定到真正原因）。另外，在后续的微扰过程中晶格的大小也发生变化，如果使用特定的K点，则会出现计算精度分布不一致的情况，此时单原子的能量可能出现一定的偏移。

cal_stress          1#计算压力 cal_force           1#计算力

添加了计算压力和计算力的选项，因为DP需要使用单个原子的受力（默认值我没关注，也许默认就开着）。

2. DPGEN init.json和machine.json设置。对稳定构性进行微扰，获得足够多的空间非稳定构型

DPGEN是作业调度软件，在做计算时需要准备他的相关文件。一个是控制程序指定功能的文件init.json。另外一个是调度服务器进行计算涉及的输入文件machine.json。

1. 其中init.json包含了微扰的相关设置，以及优化和md的相关设置。

{
    "init_fp_style": "ABACUS",
    "stages": [1, 2, 3, 4],
    "cell_type": "fcc",
    "super_cell": [ 1, 1, 1],
    "elements": ["Be", "Cu"],
    "from_poscar": true,
    "from_poscar_path": "/home/zxg/BeCu/STRU/Be/STRU_file/Be15_Cu1_hcp.STRU",
    "potcars": [
        "/home/zxg/BeCu/SG15_v1.0_ONCVPP_and_Orbitals/Be_ONCV_PBE-1.0.upf", "/home/zxg/BeCu/SG15_v1.0_ONCVPP_and_Orbitals/Cu_ONCV_PBE-1.0.upf"
    ],
    "orb_files": [
        "/home/zxg/BeCu/SG15_v1.0_ONCVPP_and_Orbitals/Be_gga_10au_100Ry_4s1p.orb", "/home/zxg/BeCu/SG15_v1.0_ONCVPP_and_Orbitals/Cu_gga_9au_100Ry_4s2p2d1f.orb"
    ],
    "relax_incar": "/home/zxg/BeCu/abacus_SET/INPUT_rlx_cell",
    "md_incar": "/home/zxg/BeCu/abacus_SET/INPUT_md",
    "relax_kpt": "/home/zxg/BeCu/abacus_SET/KPT_rlx",
    "md_kpt": "/home/zxg/BeCu/abacus_SET/KPT_md",
    "md_nstep": 10,
    "scale": [0.94, 0.96, 0.98, 1.00, 1.02, 1.04, 1.06],
    "skip_relax": false,
    "pert_numb": 1,
    "pert_box": 0.01,
    "pert_atom": 0.01,
    "coll_ndata": 5000,
    "type_map": ["Be", "Cu"],
    "_comment1": "that's all"
}

"init_fp_style": "ABACUS",#指定调用ABACUS进行计算

"stages": [1, 2, 3, 4],#对应软件控制过程，1--结构优化，2--稳定结构微扰，3--微扰结构计算，4--收集数据并转换格式。

"cell_type": "fcc",#指定晶体构型，由于后续没有让DPGEN自己创建文件，因此该命令在以上例子中没起到作用。

"super_cell": [ 1, 1, 1],#对结构进行扩胞，主要是调用dpdata。根据我的计算结果，利此处扩胞在后续结构优化时会增多计算时间（貌似atomkit也是），因此我是通过materials stuido扩的包，此处设置表示不扩胞。

"elements": ["Be", "Cu"],#元素种类，我是按先后顺序来写元素种类。

"from_poscar": true,#通过文件获取晶体结构信息。由于我的测试中扩胞会增大计算，因此我都是导入扩胞后的文件，这里true表示要导入相关文件。

"from_poscar_path": "/home/zxg/BeCu/STRU/Be/STRU_file/Be15_Cu1_hcp.STRU",
 "potcars": ["/home/zxg/BeCu/SG15_v1.0_ONCVPP_and_Orbitals/Be_ONCV_PBE-1.0.upf", "/home/zxg/BeCu/SG15_v1.0_ONCVPP_and_Orbitals/Cu_ONCV_PBE-1.0.upf"],
 "orb_files": ["/home/zxg/BeCu/SG15_v1.0_ONCVPP_and_Orbitals/Be_gga_10au_100Ry_4s1p.orb", "/home/zxg/BeCu/SG15_v1.0_ONCVPP_and_Orbitals/Cu_gga_9au_100Ry_4s2p2d1f.orb"],
 "relax_incar": "/home/zxg/BeCu/abacus_SET/INPUT_rlx_cell",
 "md_incar": "/home/zxg/BeCu/abacus_SET/INPUT_md",
 "relax_kpt": "/home/zxg/BeCu/abacus_SET/KPT_rlx",
 "md_kpt": "/home/zxg/BeCu/abacus_SET/KPT_md",
 "md_nstep": 10,

这三行时设定结构文件位置，势函数位置，原子轨道基组位置和名称。结构优化INPUT_rlx_cell文件位置，以及后续做AIMD需要的INPUT_md文件位置。还有K点设置文件，由于INPUT文件中有kspacing，因此此处K点文件不起作用（由于起始阶段使用错误，此处没删）。

"md_nstep": 10,#表示做10不AIMD的模拟，关于AIMD模拟步长和步数在下一节讨论不行。

"scale": [0.94, 0.96, 0.98, 1.00, 1.02, 1.04, 1.06],#

模拟盒子缩放尺度。

关于缩放尺度，我想进行一下说明。在缩放中尽可能涉及到我们研究项目的尺度上，而且各个制度中间的间隔不易过小。DeePMD会根据我们提供的数据进行势函数拟合，按理来说给出的准确信息越多越有利于拟合出准确的势函数。自主的说法制度可以参考前面的EOS计算和弹性模量的计算。例如，我想模拟5Gpa的压力范围，前边有计算的铍晶体样式模量信息（不一定准确，未做多次计算确认准确性），杨氏模量和单轴拉伸模量在一二百数量级，因此压缩和拉伸6%，理论上能到5Gpa（我这个你可以吧？预估的，到底能不能到我不确定）。

152.49  141.66  141.66   -0.00    0.00   -0.00 
 141.66  152.49  141.66   -0.00    0.00    0.00 
 141.66  141.66  152.49   -0.00    0.00    0.00 
  -0.00   -0.00   -0.00  123.96    0.00    0.00 
  -0.00   -0.00   -0.00    0.00  123.96   -0.00 
   0.00    0.00    0.00    0.00   -0.00  123.96 
# Bulk   Modulus BV = 145.27 GPa
# Shear  Modulus GV = 76.54 GPa
# Youngs Modulus EV = 195.32 GPa
# Poission Ratio uV = 0.28

如果要模拟熔化，建议此处范围设置大一些，主要是较长原子距离和能量跟熔化相关。

例如，如果后续我想研究磁控溅射的高能入射，那这里我可能需要考虑较大的压缩情况，让单原子的能量与最低点差值达到实验上的几十eV对。

具体模拟盒子的缩放尺度应该根据研究问题来确定。因为计算范围越大，需要计算的AIMD越多也越耗时，不过我自以为尽可能选的范围大一些好，可能前期准备数据集没用，浪费部分时间，但如果后续有用可能会降低后续的DPGEN探索时间。

"skip_relax": false,#不要跳过结构优化，我也不建议跳过，索然结构优化需要时间，但我们应该是需要他的。

"pert_numb": 1,#在上边缩放之后，对于每个缩放构型设置几个微扰，从0开始计数，1表示有两个微扰结构（可能包含原结构）（微扰结构需要多少个我没有经验，对于我项目可能涉及的结构微扰例如hcp的Be晶体，我会设置的微扰个数比较多，比如说微扰6次。而对只是可能相关的不太重要的结构，例如cubic的Be，我的围绕结构设置会比较少例如设定为1。此处经验也不足，有好建议可提出）。

"pert_box": 0.01,#晶胞盒子微扰尺度，由于有上边的缩放，因此这里使用的缩放尺度并不大（这个我没经验，可以认为是随意使用了一个较小的值） "pert_atom": 0.01,#原子坐标的微扰尺度，同样由于上边有盒子缩放因此此处的微扰设定也比较小

"coll_ndata": 5000,对于一次计算最大搜集的构型数量，一般我就几百个结构，达不到5000，如果本次计算达到这个数则随机挑选一些作为训练集。 "type_map": ["Be", "Cu"],#元素类型匹配。这里需要提示一下，如果即拟合单质也拟合掺杂或合金需要把元素匹配类型加入。刚开始模拟我没加入匹配类型，因此后续手动修改了type_map.raw文件，但手动设置容易出错，因此还是在此设定上好一些。 "_comment1": "that's all"#最后做个评论下划线开头表示这是注释没有作用，只是告诉一些相关的信息。

2. 运行DPGEN的machine.json文件

同上边autotest的计算，这里只贴出一样的machine.json，不再解释，一个文件可以用到所有的计算。据图内容目前还是建议自行理解一下

{
    "api_version": "1.0",
    "deepmd_version": "2.1.5",
     "train" :[
      {
        "command": "dp",
        "machine": {
          "batch_type": "Shell",
          "context_type": "LazyLocalContext",
          "local_root" : "./",
          "remote_root" : "./"
        },
        "resources": {
          "number_node": 1,
          "cpu_per_node": 12,
          "gpu_per_node": 1,
          "strategy"  : {"if_cuda_multi_devices": true },
          "queue_name": "zw",
          "group_size": 4,
          "source_list": ["/home/zxg/BeCu/conda_env/deepmd"]
        }
      }
    ],
    "model_devi":
      [{
        "command": "lmp",
        "machine": {
          "batch_type": "Shell",
          "context_type": "LazyLocalContext",
          "local_root" : "./",
          "remote_root" : "./"
        },
        "resources": {
          "number_node": 1,
          "cpu_per_node": 20,
          "gpu_per_node": 1,
          "strategy"  : {"if_cuda_multi_devices": true },
          "queue_name": "zw",
          "group_size": 3500,
          "source_list": ["/home/zxg/BeCu/conda_env/deepmd"]
        }
      }
    ],
    "fp":
    [  
      {
        "command": "OMP_NUM_THREADS=1 mpirun -n 12 abacus",
        "machine": {
          "batch_type": "Shell",
          "context_type": "LazyLocalContext",
          "local_root" : "./",
          "remote_root" : "./",
          "clean_asynchronously": false
        },
        "resources": 
        {
          "number_node": 1,
          "cpu_per_node": 12,
          "gpu_per_node": 1,
          "group_size": 500, 
          "para_deg": 1,
          "source_list": ["/home/zxg/BeCu/conda_env/abacus"]
        }
      }
    ]
  }  

与conda环境配套的DeePMD文件如下，即我的DeePMD是通过conda安装的

#!/bin/bash
source activate deepmd

与conda环境配套的ABACUS文件如下，即我的ABACUS也是通过conda安装的（后来conda安装出问题，我再base下安装了，就需要把machine.json文件中的"source_list"给注释掉，"_source_list"）

#!/bin/bash
source activate abacus_env

3. 进行DFT的AIMD模拟需要的INPUT文件，增多采样构型。2024.04.04

由于上边做了结构微扰，但微扰的数据可能量还不够，因此这里进行一下AIMD的模拟以补充数据，进而丰富模拟初始训练集中的构型。计算AIMD是在ABACUS中进行，需要准备结构文件STRU、计算控制文件INPUT_md、势函数文件等，并把相关文件的位置准确加入到DPGEN运行的init.json文件中。其他文件跟前期准备的文件基本一致，这里只示例ABACUS做AIMD的INPUT文件

1. ABACUS做AIMD的INPUT文件

INPUT_PARAMETERS 
#Parameters (1.General) 
calculation md 

#Parameters (2.Iteration) 
ecutwfc 100 
scf_nmax 100 

#Parameters (3.Basis) 
basis_type lcao
symmetry            0
kspacing            0.08

#Parameters (4.Smearing) 
smearing_method gauss 
smearing_sigma 0.002

#Parameters (5.Mixing) 
mixing_type pulay 
mixing_beta 0.7

cal_stress 1
cal_force  1

#MD
chg_extrap second-order
md_type nvt
md_tfirst 300
md_tlast  300 
md_nstep 10 
md_tfreq 0.5

可以看出，除了md设置相关以外，其他地方跟前边的INPUT设置完全一致，进而确保计算精度一致。（也许md和scf的默认收敛标准不一样，我用的默认值可能会出错，后续再查，自认为软件设置的默认精度是一样的）

calculation md #计算方式要设置成md，

chg_extrap second-order#这个参数我不熟，刚开始学习时在案例中拷贝（也许需要我后续确认）md_type nvt#模拟奚综选择了nvt，也是在案例中拷贝的，也许需要npt（后续再确认，希望这里影响不大）但是反过来想，如果用npt则会导致模拟盒子弛豫平衡，可能会导致起始状态做的微扰效果丢失，还是用nvt吧。 md_tfirst 300#AIMD设置起始温度为300K，根据温度给定原子速度。 md_tlast  300 #AIMD设置最终温度为300K，根据起始和终止温度调节原子速度。 md_nstep 10 #md模拟步数，这里跟DPGEN的init.json里边设置一样，他俩有一个是控制步数的，如果想做更多的模拟可以修改步数

md_tfreq 0.5#模拟过程中速度更新对应的频率0.5。

2. 关于AIMD的讨论

AIMD模拟中我感觉需要讨论的是md模拟的步数，以及速度跟新md模拟的步长，模拟步长我没经验暂不讨论。

下边讨论一下模拟的总步数。

关于AIMD的步数问题，我认为是想通过AIMD来获与初始构型密切关联的构型，进而在DP拟合中给出势函数能量的准确导数信息，并获取比较准确的受力。而经过微扰的结构一般不会在稳定构型附近，因此在做AIMD时，我们选择的模拟步长最好能够让一些AIMD模拟到平衡结构，这样应该有利于后续DPGEN的探索时间降低。关于多少步模拟能够达到平衡结构，我觉得可以通过输出最大力来验证。下边给一个画最终计算结果最大力的python脚本，该脚本是我把所有init计算的结果都复制到deepmd_file下的对应的文件夹之后再做的处理脚本（脚本参考于贵州师范大学张为老师）。

即在使用脚本之前需要先执行，然后再deepmd_file下会生成很多个init_bulk获得的文件夹。

cp -r /you/path/to/deepmd /order/path/to/savedataset/deepmd_file/

#!/usr/bin/env python
import os
import dpdata
import numpy as np
import matplotlib.pylab as plt
folder_path = '/home/zxg/BeCu/dpgennew/deepmd_file/'
file_list = os.listdir(folder_path)
for file in file_list:
    print(file)
    folder_list = os.path.join(folder_path, file, 'deepmd/')
    print(folder_list)
    data=dpdata.LabeledSystem(folder_list, fmt="deepmd/raw")
    forces = np.absolute(data["forces"])
    ff = []
    for i in range(len(data)):
        ff.append(np.max(forces[i]))
    plt.plot(ff,'-o')
    plt.savefig(file+'max.png')
    plt.clf()

脚本解释请自行咨询chatgpt，大致为读取对应文件夹，并在对应文件夹下利用dpdata读取数据，然后获取数据中的最大力并画图。然后删除数据，再整下一个文件夹（如果不删除数据，则下次循环所有的数据都会输出）

当然脚本修改一下就可以读取单次init_bulk数据的值并画图

#!/usr/bin/env python
import os
import dpdata
import numpy as np
import matplotlib.pylab as plt
folder_path = '/home/zxg/BeCu/dpgennew/deepmd_file/deepmd/'
print(folder_path)
data=dpdata.LabeledSystem(folder_path, fmt="deepmd/raw")
forces = np.absolute(data["forces"])
ff = []
for i in range(len(data)):
    ff.append(np.max(forces[i]))
plt.plot(ff,'-o')
plt.savefig('max_f_Be16.png')
plt.clf()

4. 收集计算结果，提取不同构型的力能量信息，并准备初始训练集。

这个过程是自动化的，目前无需人为参与，是DPGEN的流程操作过程，使用dpdata搜集数据并进行格式转化，此处目前没什么想法。

但是对于初学者来说，最好知道dpdata如何转换相关的数据，在 《DFT 到 MD｜超详细「深度势能」材料计算上手指南》中的3.1节有介绍。

同样DPGEN的input.json和machine.json同3.12节。