4.1.3.3 构型探索相关分析2024.04.19, 2024.04.25

DP势函数拟合之后，由于初始构型准备时可能存在不完备的情况，因此需要对训练集进行进一步的补充。DPGEN会调用LAMMPS利用前期拟合的四个势函数进行相关结构动力学演化计算，并计算四个势函数在同一个构型下的原子受力的不同。如果原子受力偏差大于某个设定值则认为训练的势函数不能准确描述该结构的构型。涉及到的关键字是：

"model_devi_f_trust_lo": 0.15, 当然在探索过程可能会出现非物理的情况，比如原子距离太近等，这种情况可能也不是我们后期模拟所需要的，因此需要设定一个力的偏差值上限，即大于某个力之后认为物理结构超出研究范围，不再进行补充（注：此设定更多的是为了节省计算资源，放弃没必要训练的场景，降低计算量）。

"model_devi_f_trust_hi": 0.35,

其他解释参见上边章节 4.1.1.1.3 模型探索相关设置

在进行LAMMPS动力学模拟时会生成对应的文件，在inter*中的01.model_devi包括task.*.*。在task下是具体执行一次LAMMPS模拟的输入和输出文件，其中model_devi.out包含了四个势函数计算的结构在不同步数下的力偏差值。

├── task.000.000001
│   ├── 75283215e7e2226df31d3c4113d90ee5c6598934_task_tag_finished
│   ├── all.lammpstrj
│   ├── conf.lmp -> ../confs/000.0000.lmp
│   ├── dpgen.restart.10000
│   ├── input.lammps
│   ├── job.json
│   ├── log.lammps
│   ├── model_devi.log
│   └── model_devi.out

由于DPGEN中设定了不同的温度，压力，和探索构型，因此每次迭代有很多次LAMMPS的计算数据需要处理，此时就需要自行修改脚本针对某次迭代进行画图。这里先给一个单次迭代中某个探索构型的最大偏差力分析图。python脚本不再解释，可自行查阅chatgpt或智普清言

#!/usr/bin/env python
import os
import numpy as np
import matplotlib.pyplot as plt

for i in range(8,9):    # Specify the number of iterations
    for j in range(0, 4):    # Specify the STRU for LAMMPS model_devi
        max_devi_f_values=[]
        for k in range(0,20): # 20 LAMMPS tasks munber for before stru
            directory = "/home/zxg/BeCu/dpgen/run/iter.{:06d}/01.model_devi/task.{:03d}.{:06d}".format(i, j, k%20)
            file_path = os.path.join(directory, "model_devi.out")
            data = np.genfromtxt(file_path, skip_header=1, usecols=4)
            max_devi_f_values.append(data)
    
        # Convert the list to a numpy array
        max_devi_f_values = np.concatenate(max_devi_f_values)
    
        # Use numpy.histogram() to calculate the frequency of each value
        hist, bin_edges = np.histogram(max_devi_f_values, range=(0, 0.9), bins=100)
    
        # Normalize the frequency to percentage
        hist = hist / len(max_devi_f_values) * 100
    
        # Use matplotlib.pyplot.plot() to plot the frequency of each value
        plt.plot(bin_edges[:-1], hist, label = 'iter{:02d}__dev{:02d}'.format(i,j))
        plt.legend()
        plt.xlabel("max_devi_f eV/Å")
        plt.ylabel("Distribution %")
    
        with open(f'./iter{i:02d}_dev{j:02d}dist-max-devi-f.txt', 'w') as f:
            f.write("{:>12} {:>12}\n".format("bin_edges", "hist"))
            for x, y in zip(bin_edges[:-1], hist):
                f.write('{:>12.3f} {:>12.3f}\n'.format(x, y))
        plt.savefig(f'iter{i:02d}_dev{j:02d}.png',dpi=300)
        plt.clf()

由于我迭代时使用了多次势函数的调整，包括探索步长，训练步数等，进而导致训练力规律并未展现出来，一般情况，迭代是从右向左迭代的，例如inter00的蓝线向inter06的浅紫线必经收敛。

以上脚本最后两行如果缩进至最内层循环，且控制变量在最终的迭代可分析最后的力的偏差，结果如下

Be108_fcc

Be128_hcp

Be16_hcp

Be32_fcc

结果图是迭代8次之后的力的偏差，由于初始集准备的比较充足，因此在后续迭代过程中很快就会得到一个准确度在制定数据之下的DP势函数，降低了迭代次数。当然，可以修改脚本把不同的迭代次数画到一个图中。

关于迭代中力的偏差，刚开始偏差值一般会很大，随着迭代次数增多，偏差值会一直降低，等降低至某个收敛值后不再降低，这时候基本就训练好DP势函数了。

4. 针对特定问题增加DPGEN迭代次数

通过势函数质量分析，可大致判断势函数的准确性，进而可针对特定问题进行势函数准确性的验证。