0. 写在前边

本教程首先介绍了DeePKS训练数据的格式，然后介绍了一种通过DeePKS init功能进行训练数据生产的方法，该方法只适合于使用ABACUS进行高精度的训练数据的生产，需要用户对DeePKS流程有基本的了解，新手可以先阅读后续的教程熟悉DeePKS流程，在需要生产自己的训练数据时再查看本教程。

1. DeePKS训练数据介绍

DeePKS的训练数据主要包括两部分，一部分是算例体系的结构文件，主要包括原子种类及坐标信息以及晶胞的信息；另一部分则是算例的标签数据，例如能量、力、应力和轨道等。结构文件的格式及各种属性标签的对应格式总结如下：

1.1 算例体系的结构文件

1.1.1 原子种类及坐标信息atom.npy（推荐格式）

atom.npy 文件的形状为 [nframes, natoms, 4]：

• nframes 指的是算例体系的帧（结构）数量；
• natoms 指的是算例体系的原子数量，例如，对于单一水分子系统，natoms = 3；
• 最后一维 4 对应于给定原子的核电荷及其 xyz 坐标（可以是笛卡尔坐标形式或分数坐标形式，需要在 DeePKS的配置文件scf_abacus.yaml (DeePKS中配置ABACUS scf计算的参数文件）中用关键字 coord_type 指定）。

其中一帧的结构文件示例如下：

[[55, 3.302000045776367, 10.319330215454102, 14.370149612426758], # 每一行表示一个原子的种类和坐标信息，第一个整数表示核电荷数（原子序数），后边三个数字表示原子坐标
  [55, 1.9708499908447266, 9.589559555053711, 6.998380184173584],
  [55, 3.530519962310791, 5.52239990234375, 15.28225040435791],
  [55, 2.157900094985962, 5.562459945678711, 3.440109968185425],
  [82, 1.8799699544906616, 2.1771700382232666, 18.001850128173828],
  [82, 4.564519882202148, 10.913680076599121, 1.2713600397109985],
  [82, 0.9406599998474121, 6.323229789733887, 11.14264965057373],
  [82, 4.157649993896484, 2.7756199836730957, 7.601940155029297],
  [53, 4.622320175170898, 4.023340225219727, 18.982349395751953],
  [53, 1.986649990081787, 9.189929962158203, 1.4250600337982178],
  [53, 3.3898699283599854, 8.489720344543457, 10.775850296020508],
  [53, 2.120500087738037, 1.3380500078201294, 9.347249984741211],
  [53, 1.1738500595092773, 1.945930004119873, 1.6067700386047363],
  [53, 4.532810211181641, -0.024819999933242798, 17.70655059814453],
  [53, 2.0408899784088135, 5.284709930419922, 7.748879909515381],
  [53, 3.8335800170898438, 4.460169792175293, 11.49275016784668],
  [53, 1.1351399421691895, 2.606300115585327, 15.121549606323242],
  [53, 4.954619884490967, 8.381369590759277, 4.626279830932617],
  [53, 0.8476099967956543, 7.216929912567139, 14.072750091552734],
  [53, 4.062069892883301, 2.1832399368286133, 4.9308600425720215]],

注意

如果 atom.npy 中保存的坐标单位为 Bohr，则在 scf_abacus.yaml 中应将 lattice_constant 设置为 1。如果 atom.npy 中保存的坐标单位为 Angstrom，则 lattice_constant 应设置为 1.8897259886。如果使用分数坐标，lattice_constant 也需要相应设置。有关更多详细信息，请参见 ABACUS 文档。

1.1.2 晶胞信息box.npy

box.npy 指定每个帧的晶格矢量，形状为 [nframe, 9]，是将晶胞参数的9个数字横向拼接形成。 示例如下：

[5.092175483703613, 0, 0, 0.2556837499141693, 11.27470874786377, 0, 0.3732239007949829, -0.425565630197525, 19.532777786254883]

其单位与atom.npy同样需要考虑单位的转换问题，有关更多详细信息，请参见 ABACUS 文档。

1.2 算例体系的标签

DeePKS目前支持训练的标签包括能量、力、应力和轨道（能带），展开介绍如下：

为了训练 DeePKS 模型，需要感兴趣系统的目标能量，其格式应遵循结构文件的格式。还可以训练额外的属性，包括力、应力和带隙。请注意，带隙标签对应于每个 k 点的价带与导带之间的能量差，目前此标签仅适用于具有偶数电子的半导体。结构文件的格式（以 atom.npy 为例）及各种属性标签的对应格式总结如下：

1.2.1 能量标签energy.npy

energy.npy的形状为[nframes, 1]，即每个算例对应一个energy值，其单位为Hartree，注意ABACUS输出的能量单位通常为Ry，需要进行转换，1Ry = 0.5 Hartree。

1.2.2 力标签force.npy

force.npy的形状为[nframes, natoms, 3]，即每个算例中每一个原子的三维坐标上的受力值，其单位为Hartree/Bohr，注意ABACUS输出的能量单位通常为Ry/Bohr，需要进行转换，1Ry = 0.5 Hartree。

1.2.3 应力标签stress.npy

stress.npy的形状为[nframes, 9]，即每个算例的应力矩阵，其单位为Hartree，注意ABACUS输出的能量单位通常为Ry，需要进行转换，1Ry = 0.5 Hartree。

1.2.4 轨道（能带）标签orbital.npy

orbital.npy的形状为[nframes,nkpt,1]，即每个算例的每个k点对应的数值，其单位为Hartree，注意ABACUS输出的能量单位通常为Ry，需要进行转换，1Ry = 0.5 Hartree。

2. DeePKS训练数据生产流程

目前DeePKS最重要的特点是可以在不显著增加时间成本的情况下大幅度提高精度，考虑到DeePKS自身训练的算力成本，目前DeePKS主要有两大应用：

1. 用于生产高精度的DeePMD的训练数据（研究体系中的原子数、元素不发生变化）； 2. 大规模生产特定材料体系的高精度计算数据（研究体系中的原子数、元素均可能变化）；

目前针对用于生产高精度的DeePMD的训练数据这一应用已有Notebook从 DFT 先去 DeePKS 再到 DeePMD | DeePKS基础篇和从 DFT 先去 DeePKS 再到 DeePMD | DeePKS案例篇 + 增强采样可以参考，本示例将针对第二种应用的训练数据生产进行介绍。

2.1 算例结构文件的准备（atom.npy,box.py)

由于大规模生产特定材料体系的高精度计算数据时算例的研究体系中的原子数、元素均可能变化，不适合使用aimd的方法生成轨迹来生产数据。这里需要用户自己收集好相关的算例结构文件，编写脚本将其整理成atom.npy和box.npy格式。用户也可以在开源平台aissquare.com下载相关数据集，本篇教程中使用的算例文件为直接从科学智能广场下载的开源钙钛矿数据集。

用户可以根据自己的需要自定义结构文件数据集，训练数据的选取应尽可能覆盖所有研究体系。

2.2 算例体系标签的生产

算例体系标签的生产简而言之就是运行对应的高精度泛函的DFT计算，例如通过ABACUS的HSE计算，收集计算结果并将其整理成energy.npy,force.npy,stress.npy以及orbital.npy等。用户可以通过批处理库例如dpdata来实现这一过程。

DeePKS集成了批量提交计算任务的功能，在本教程中，我们提供一种简单的利用DeePKS init过程通过ABACUS生产算例标签的方法：

1. 修改scf_abacus.yaml文件中的init_scf_abacus部分（初始轮次对应的参数）：（主要修改泛函种类、输出力、应力、能带标签的参数）

示例如下：

init_scf_abacus:
  orb_files: ["Cs_gga_10au_100Ry_4s2p1d.orb","Pb_gga_7au_100Ry_2s2p2d1f.orb","I_gga_7au_100Ry_2s2p2d1f.orb","C_gga_7au_100Ry_2s2p1d.orb","N_gga_7au_100Ry_2s2p1d.orb","H_gga_6au_100Ry_2s1p.orb"]
  pp_files: ["Cs_ONCV_PBE-1.0.upf","Pb_ONCV_PBE-1.0.upf","I_ONCV_PBE-1.0.upf","C_ONCV_PBE-1.0.upf","N_ONCV_PBE-1.0.upf","H_ONCV_PBE-1.0.upf"]
  proj_file: ["jle.orb"]
  # ntype: 3
  ecutwfc: 100
  scf_thr: 1e-7
  scf_nmax: 50
  dft_functional: "hse" #高精度泛函种类
  gamma_only: 0
  kspacing: 0.1
  cal_force: 1 #生产力的标签
  cal_stress: 1 #生产应力的标签
  deepks_bandgap: 1  #生产轨道（能带）的标签
  lattice_constant: 1.88972613
  lattice_vector: [[28, 0, 0], [0, 28, 0], [0, 0, 28]]
  coord_type: "Cartesian"
  #cmd args
  run_cmd : "mpirun"
  abacus_path: "abacus"

2. 正常提交DeePKS任务

bash run.sh（本地运行） 或者 run_dpdispatcher.sh （提交服务器运行）

取消注释可进行hse高精度标签数据生产，大约需要150元机时费用。

3. init轮次任务结束后手动停止任务，在systems目录下运行stats.py脚本

import os
import numpy as np

# 根目录
systems_dir = "."

# 处理每个group并输出energy.npy、force.npy、orbital.npy和stress.npy文件
def process_group(group_dir):
    energy_list = []
    force_list = []
    orbital_list = []
    stress_list = []
    print(f"正在处理组：{group_dir}")
    
    # ABACUS目录路径
    abacus_dir = os.path.join(group_dir, "ABACUS")
    
    # 确保ABACUS目录存在
    if not os.path.isdir(abacus_dir):
        print(f"在以下目录中未找到ABACUS目录：{group_dir}")
        return
    
    # 遍历ABACUS目录下的子目录，并按照数字顺序排序
    subdirs = [subdir for subdir in os.listdir(abacus_dir) if os.path.isdir(os.path.join(abacus_dir, subdir))]
    subdirs_sorted = sorted(subdirs, key=lambda x: int(x))  # 将子目录按数字顺序排序
    
    for subdir in subdirs_sorted:
        subdir_path = os.path.join(abacus_dir, subdir)
        
        out_abacus_dir = os.path.join(subdir_path, "OUT.ABACUS")
        
        # 检查并处理deepks_etot.npy文件（能量）
        deepks_etot_path = os.path.join(out_abacus_dir, "deepks_etot.npy")
        if os.path.exists(deepks_etot_path):
            etot = np.load(deepks_etot_path)  # 示例：[-1517.8723503798371]
            energy_hartree = etot / 2  # Ry 转换为 Hartree
            energy_list.append(energy_hartree)
        else:
            print(f"  在以下路径未找到deepks_etot.npy：{deepks_etot_path}")
        
        # 检查并处理deepks_ftot.npy文件（力）
        deepks_ftot_path = os.path.join(out_abacus_dir, "deepks_ftot.npy")
        if os.path.exists(deepks_ftot_path):
            force = np.load(deepks_ftot_path)
            force_hartree = force / 2  # Ry 转换为 Hartree
            force_list.append(force_hartree)
        else:
            print(f"  在以下路径未找到deepks_ftot.npy：{deepks_ftot_path}")
        
        # 检查并处理deepks_otot.npy文件（轨道）
        deepks_otot_path = os.path.join(out_abacus_dir, "deepks_otot.npy")
        if os.path.exists(deepks_otot_path):
            otot = np.load(deepks_otot_path)  # 假设为三维数组
            otot_hartree = otot / 2.0  # Ry 转换为 Hartree
            orbital_list.append(otot_hartree)
        else:
            print(f"  在以下路径未找到deepks_otot.npy：{deepks_otot_path}")
        
        # 检查并处理deepks_stot.npy文件（应力）
        deepks_stot_path = os.path.join(out_abacus_dir, "deepks_stot.npy")
        if os.path.exists(deepks_stot_path):
            stot_upper = np.load(deepks_stot_path)  # 示例：[sxx, sxy, sxz, syy, syz, szz]
            
            # 检查数据形状，如果是一维的，则需要扩展为二维
            if stot_upper.ndim == 1:
                stot_upper = stot_upper.reshape(1, -1)
            
            # 转换单位，从Ry到Hartree
            stot_upper = stot_upper / 2.0
            
            # 对于每个样本，构建完整的应力张量扁平列表
            for upper in stot_upper:
                # 上三角部分的元素
                sxx = upper[0]
                sxy = upper[1]
                sxz = upper[2]
                syy = upper[3]
                syz = upper[4]
                szz = upper[5]
                # 对称性
                syx = sxy
                szx = sxz
                szy = syz

                # 构建扁平列表
                stress_flat = [sxx, sxy, sxz, syx, syy, syz, szx, szy, szz]
                # 将扁平列表加入 stress_list
                stress_list.append(stress_flat)
        else:
            print(f"  在以下路径未找到deepks_stot.npy：{deepks_stot_path}")
    
    # 保存能量数据
    if energy_list:
        energy_output_path = os.path.join(group_dir, "energy.npy")
        np.save(energy_output_path, np.array(energy_list))
        print(f"已将energy.npy保存到 {energy_output_path}")
    else:
        print(f"在{group_dir}中未找到任何能量数据")
    
    # 保存力数据
    if force_list:
        force_output_path = os.path.join(group_dir, "force.npy")
        np.save(force_output_path, np.array(force_list))
        print(f"已将force.npy保存到 {force_output_path}")
    else:
        print(f"在{group_dir}中未找到任何力数据")
    
    # 保存轨道数据
    if orbital_list:
        orbital_array = np.array(orbital_list)
        orbital_output_path = os.path.join(group_dir, "orbital.npy")
        np.save(orbital_output_path, orbital_array)
        print(f"已将orbital.npy保存到 {orbital_output_path}")
    else:
        print(f"在{group_dir}中未找到任何轨道数据")
    
    # 保存应力数据
    if stress_list:
        stress_array = np.array(stress_list)  # 形状：(n_samples, 9)
        stress_output_path = os.path.join(group_dir, "stress.npy")
        np.save(stress_output_path, stress_array)
        print(f"已将stress.npy保存到 {stress_output_path}")
    else:
        print(f"在{group_dir}中未找到任何应力数据")

# 遍历systems目录下的所有group
for group in os.listdir(systems_dir):
    group_dir = os.path.join(systems_dir, group)
    
    # 只处理以"group"开头的目录
    if os.path.isdir(group_dir) and group.startswith('group'):
        process_group(group_dir)
    else:
        print(f"跳过 {group_dir}，因为它不是一个group目录")

这样即可在systems下的子目录中构建所有标签数据。