MACE能量与偶极联用：预测和短时MD

希望沿轨迹同时读取能量、原子力和偶极，可以使用 EnergyDipolesMACE。本例以真实H₂O DFT标签训练一个小模型，再在P100上完成20步、总2 fs的Velocity Verlet演示，保存21帧完整坐标、动量、能量、力和偶极。

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

bash run.sh evaluate
bash run.sh trajectory
head -n 6 results/trajectory/MD-observables.csv

前者用固定模型做3个独立测试构型的推断，后者读取已有21帧轨迹，不再跑一次MD。

2 训练设置与三种模型的区别#

完整configs/energy-dipole.yml：

train_file: ../../data/train.xyz
valid_file: ../../data/valid.xyz
energy_key: REF_energy
forces_key: REF_forces
dipole_key: REF_dipole
polarizability_key: REF_polarizability
charges_key: REF_charges
total_charge_key: REF_total_charge
E0s: average
batch_size: 5
valid_batch_size: 3
max_num_epochs: 200
default_dtype: float64
seed: 20261108
ema: true
ema_decay: 0.99
swa: false
num_workers: 0
device: cuda
log_dir: logs
model_dir: models
checkpoints_dir: checkpoints
results_dir: results
lr: 0.005
weight_decay: 5.0e-10
config_type_weights: '{''Water'':1.0}'
r_max: 4.0
hidden_irreps: 16x0e+16x1o+16x2e
MLP_irreps: 16x0e+16x1o
num_interactions: 2
correlation: 2
max_ell: 2
max_L: 2
clip_grad: 10.0
patience: 200
scheduler_patience: 30
eval_interval: 1
save_cpu: true
compute_stress: false
name: water-energy-dipole
model: EnergyDipolesMACE
loss: energy_forces_dipole
error_table: EnergyDipoleRMSE
dipole_weight: 100.0
polarizability_weight: 100.0
energy_weight: 1.0
forces_weight: 10.0
num_radial_basis: 8
num_cutoff_basis: 5

能量单位eV/整个分子、力eV/Å、偶极e·Å。能量、力、偶极权重分别1、10、100；权重用于平衡数值损失，不是将偶极转成其他单位。测试不参与训练。

sbatch train.slurm energy-dipole
# 完成后
ls work/train/energy-dipole/models

实际训练200轮，附件模型是已训练、可CPU读取的自有教学权重。这个模型具有零场能量头和偶极头，不是随外加电场变化的能量函数；不能仅凭同时输出偶极，就认定mu=−dE/dField或可以模拟强电场驱动。

3 用ASE同时读取三个量#

python - <<'PYCODE'
import numpy as np
from ase.io import read
from mace.calculators import MACECalculator
a = read('qualified/data/test.xyz', index=0)
for key in list(a.info):
    if key.startswith('REF_'): del a.info[key]
for key in list(a.arrays):
    if key.startswith('REF_'): del a.arrays[key]
a.set_initial_charges(np.zeros(len(a)))
a.calc = MACECalculator(model_paths='models/energy-dipole.model',
    model_type='EnergyDipoleMACE', device='cpu', default_dtype='float64')
print('E [eV/molecule]:', a.get_potential_energy())
print('F [eV/A]:', a.get_forces())
print('mu [e A]:', a.get_dipole_moment())
PYCODE

独立测试3构型：能量RMSE 4.2453493 meV/原子、全部原子笛卡尔力分量RMSE 0.02368355362 eV/Å、整分子偶极向量RMSE 0.00316133341 Debye。能量相减得到力的有限差分、旋转协变及平移不变性也由配套脚本独立检查。

4 查看实际2 fs轨迹#

本例从构型5出发，Maxwell初速度目标25 K，固定随机种子20261108；去掉质心动量后不重新缩放温度，因此不能把后续轨迹称为恒定25 K。采用NVE、时间步0.1 fs、20步，没有温控器。全轨迹总能量极差约 6.53837×10⁻⁶ eV/分子。

直接读取全精度ASE轨迹：

python - <<'PYCODE'
from ase.io import read
frames = read('molecular-MD/MD.traj', ':')
for i, a in enumerate(frames):
    print(i * 0.1, a.get_potential_energy(), a.get_kinetic_energy(),
          a.get_dipole_moment())
PYCODE

MD.traj保留所有原子坐标和动量，以及同一步计算的能量、力、偶极；MD-FIELDS.json是完整数值副本。计算开始前已移除DFT参考标签，轨迹中的属性来自MACE预测。

要按相同方法做一次新的动力学演示，可提交：

sbatch new-md.slurm

脚本在新目录results/new-MD/生成轨迹，核心操作如下：

from ase import units
from ase.md.verlet import VelocityVerlet
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary
MaxwellBoltzmannDistribution(atoms, temperature_K=25,
                            rng=np.random.default_rng(20261108))
Stationary(atoms, preserve_temperature=False)
dyn = VelocityVerlet(atoms, timestep=0.1 * units.fs)
dyn.run(20)

本片段接在模型计算器绑定和参考标签清理之后；附件tools/run_new_md.py包含可直接执行的完整程序和逐步保存逻辑。

5 能量稳定之外，还应看什么#

本训练集只含两根O–H等长结构；实际短轨迹已经出现不等长运动，两键最大相差约0.00377 Å，键长范围0.9600–0.97294 Å，夹角约103.716–104°。这说明即使总能量数值稳定，也不等于力场已在这些运动方向获得DFT精度验证。

因此本例用于学习多属性接口、积分与轨迹读取，不能当作通用水势，也不从2 fs轨迹计算或宣传IR谱。若目标是液态水、Raman、IR或更长动力学，应补充非对称伸缩、热采样、分子间相互作用等相应DFT标签，再独立检查力、偶极、极化率和时间步/采样收敛。