MACE分子属性：准备偶极与极化率标签

MACE除了能量和力，还可以学习分子偶极和极化率。本组教程从标签开始，再分别训练偶极模型、偶极–极化率模型及能量–偶极模型。标签来自鸢算科技集群实际完成的 CP2K2026.2 PBE/GTH/DZVP计算，模型使用 MACE0.3.16，并通过ASE3.29读取。

1 下载并查看结构与标签#

下载本组完整输入、真实DFT标签、三个教学模型、轨迹和绘图脚本。

python3 -m zipfile -e mace-response-practical-wave315-*.zip .
cd mace-response-practical
bash run.sh

bash run.sh从保存的数值生成results/figures/下五张图，不重复DFT、训练或MD。以下Python操作在mace/0.3.16环境中执行；推断、训练与动力学请在已分配的计算资源内运行，登录节点用于编辑和提交。

module purge
module load mace/0.3.16

本例16种中性孤立H₂O构型：两根O–H键等长，键长取0.92、0.96、1.00、1.04 Å，夹角取98、104、110、116°。每个构型都得到能量、三个原子的全力向量、总偶极和完整3×3极化率张量。另有12次有限电场、4次原子位移计算用于独立差分，不进入训练集。

2 第一步：准备一份完整DFT输入#

inputs/geom-05-zero/input.inp的完整内容如下。它是固定结构单点，不会自动优化原子，也不优化晶胞：

&GLOBAL
 PROJECT response
 RUN_TYPE ENERGY_FORCE
 PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
 METHOD QS
 &DFT
  &PRINT
   &MOMENTS ON
    PERIODIC .FALSE.
    REFERENCE COM
   &END MOMENTS
  &END PRINT
  BASIS_SET_FILE_NAME ../../DATA/BASIS_MOLOPT
  POTENTIAL_FILE_NAME ../../DATA/GTH_POTENTIALS
  &QS
   EPS_DEFAULT 1.0E-12
  &END QS
  &MGRID
   CUTOFF 600
   REL_CUTOFF 50
  &END MGRID
  &SCF
   EPS_SCF 1.0E-8
   MAX_SCF 100
   SCF_GUESS ATOMIC
   &OT
    MINIMIZER DIIS
    PRECONDITIONER FULL_SINGLE_INVERSE
   &END OT
  &END SCF
  &XC
   &XC_FUNCTIONAL PBE
   &END XC_FUNCTIONAL
  &END XC
  &POISSON
   PERIODIC NONE
   POISSON_SOLVER MT
  &END POISSON
 &END DFT
 &SUBSYS
  &CELL
   ABC 12 12 12
   PERIODIC NONE
  &END CELL
  &COORD
  O 6.000000000000 6.000000000000 6.000000000000
  H 6.756490323462 6.000000000000 6.591035016313
  H 5.243509676538 6.000000000000 6.591035016313
  &END COORD
  &KIND O
   BASIS_SET DZVP-MOLOPT-SR-GTH
   POTENTIAL GTH-PBE-q6
  &END KIND
  &KIND H
   BASIS_SET DZVP-MOLOPT-SR-GTH
   POTENTIAL GTH-PBE-q1
  &END KIND
 &END SUBSYS
 &PRINT
  &FORCES ON
  &END FORCES
 &END PRINT
 &PROPERTIES
  &LINRES
   EPS 1.0E-8
   MAX_ITER 100
   PRECONDITIONER FULL_SINGLE_INVERSE
   &POLAR ON
    DO_RAMAN .TRUE.
    PERIODIC_DIPOLE_OPERATOR .FALSE.
   &END POLAR
  &END LINRES
 &END PROPERTIES
&END FORCE_EVAL

12 Å是计算盒边长，分子不周期；POISSON PERIODIC NONE与MT边界共同处理孤立体系。CUTOFF 600、REL_CUTOFF 50单位为Ry，EPS_SCF 1e-8控制电子自洽。LINRES/POLAR计算静态电子极化率，这不是有限温度介电常数，也不含核振动响应。基组和赝势参考数据已放在DATA/，输入从实际计算目录用../../DATA读取。

先查看附件已有结果，再按需要提交新的单点：

cat inputs/geom-05-zero/input.inp
sbatch new-dft.slurm geom-05-zero
# 作业结束后
rg 'SCF run converged|Polarizability tensor|PROGRAM ENDED' cases/geom-05-zero/output.log

new-dft.slurm加载cp2k/2026.2-p100，申请一张P100及两个CPU线程。实际执行命令是cp2k.psmp -i input.inp -o output.log。脚本拒绝覆盖已有同名计算目录；需要再做不同设置时，另存新目录和新输入。

3 第二步：把原生单位转成训练单位#

采用同版本CP2K原生CODATA2006常数：1 Hartree=27.2113838565563 eV，1 Bohr=0.52917720859 Å，1原子偶极单位=2.5417462288785617 Debye。因此本数据1 e·Å=4.803204271875351 Debye。极化率由原子单位乘Bohr³转成ų。不要把Debye数直接写入这里的e·Å标签，也不要期待模型自动推断单位。

查看标签和拆分，不需重新算DFT：

head -n 5 qualified/data/train.xyz
cat qualified/data/SPLIT.json
python - <<'PYCODE'
from ase.io import read
frames = read('qualified/data/all.xyz', ':')
a = frames[5]
print(a.info['REF_energy'])
print(a.arrays['REF_forces'])
print(a.info['REF_dipole'])
print(a.info['REF_polarizability'])
PYCODE

构型5的参考能量为−468.56645399746344 eV/分子，偶极主要沿z、大小约1.793129 Debye；极化率对角为1.1893882635、0.6920594778、0.9767916923 ų。分子在xz平面，y是面外方向。微小非对角量与张量反对称数值残差均保留，不用事后对称化修改原始标签。

4 第三步：用有限电场检查标签方向与符号#

附件inputs/field-g05-*给出完整电场输入。CP2K外势写为A*(X-X0)等形式，物理电场E=−A。本例分别沿x、y、z，使用A=±0.0005、±0.001原子单位；同一分子保持原位。因而：

mu_i = [Energy(+h_i) - Energy(-h_i)] / (2h)
alpha_ii = -[Energy(+h_i) - 2Energy(0) + Energy(-h_i)] / h^2
alpha_ij = -[mu_i(+h_j) - mu_i(-h_j)] / (2h)

这里能量和偶极先转换为原子单位。构型5原生αxx=8.026389951259 a.u.，h=0.0005的能量曲率给8.026392180227；αzz=6.591717157951，h=0.001的能量曲率给6.591727135685。qualified/FD-CONTROLS.json保留两步长、三方向及全部张量分量，原子力另用±0.0005、±0.001 Å坐标位移检查F=−dE/dr。

5 第四步：按构型拆分，再进入模型训练#

训练10个构型[0,3,12,15,5,6,9,10,1,14]，验证3个[2,7,8]，测试3个[4,11,13]。三种模型共用同一拆分，测试文件不传给训练程序。extXYZ保留17位有效数字，避免力比较被坐标输出精度影响。

这些都是同一分子附近、两键等长的窄范围结构，不覆盖断键、溶剂、其他电荷态或任意不对称伸缩。需要实际生产模型时，应先扩展目标结构分布，再做独立测试。

接下来阅读偶极训练、极化率训练和能量–偶极联用。