GPAW超精细耦合与氢谱线

本篇演示从自旋极化计算提取超精细耦合张量，再换算氢原子的谱线频率。完整操作包括：建立H原子模型、计算电子态、调用GPAW原生工具、乘核g因子、按磁矩归一化、导出张量并绘图。

这是小模型的使用参考，不能把一次LDA或PBE计算当作精密预测实验21 cm谱线的方法。

1 下载、提交和输出#

下载超精细耦合实操包，保存为gpaw-hyperfine.zip。从登录节点提交时，解压后使用sbatch job.slurm；已进入分配好的CPU资源时，依次执行：

python3 -m zipfile -e gpaw-hyperfine.zip .
cd gpaw-hyperfine
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

任务使用1个MPI进程、1线程，运行6个条件。每个条件的原生输出与完整数值数组在results/条件名/；汇总为results/VALIDATION.json，两张图在plots/。

2 自旋极化H原子输入#

模型是中性H：1个电子，初始磁矩+1，位于立方晶胞中心。使用周期平面波模型，以真空间隔近似孤立原子，不能把周期镜像误差当作自动消失。

from ase import Atoms
from gpaw import GPAW, PW, FermiDirac

L = 6.0
atoms = Atoms('H', positions=[[L/2, L/2, L/2]],
              cell=[L, L, L], pbc=True, magmoms=[1])
atoms.calc = GPAW(
    legacy_gpaw=True, mode=PW(600), xc='PBE',
    spinpol=True, nbands=2, symmetry='off',
    occupations=FermiDirac(0.01, fixmagmom=True),
    convergence={'energy': 1e-8,
                 'density': 1e-6,
                 'eigenstates': 1e-9}, maxiter=200,
    txt='native.txt')
energy = atoms.get_potential_energy()
atoms.calc.write('state.gpw', mode='all')

附件比较PBE的400、600、800 eV截止能；另算600 eV、8 Å晶胞，600 eV、磁矩−1，以及LDA、400 eV。固定磁矩占据在这里用于维持同一个H原子的自旋状态，不是给材料随意指定磁矩的通用做法。

3 提取张量并完成单位换算#

原生hyperfine_parameters对每个原子返回一个3×3张量，单位为eV。该结果尚需乘对应核的g因子。对于本例总磁矩为±1的单H原子，再除以总磁矩，比较同一个正号约定下的耦合常数：

import numpy as np
from ase import units
from gpaw.hyperfine import hyperfine_parameters

raw_A = hyperfine_parameters(atoms.calc)[0]  # eV，尚未乘核g因子
g_proton = 5.586                            # 本例使用的1H核g因子
M = atoms.calc.get_magnetic_moment()
A = raw_A * g_proton / M
a_iso = np.trace(A) / 3
A_aniso = A - np.eye(3) * a_iso
frequency_Hz = a_iso * units._e / units._hplanck
wavelength_cm = units._c / frequency_Hz * 100

a_iso是各向同性部分，A_aniso是各向异性部分。对近似球对称的H基态，三个对角分量应接近相等，非对角项应很小。不能直接拿一个未经核g因子处理的数值当作频率，也不能忘记eV到焦耳的转换。

多原子或反铁磁体系不能机械照搬本例的总磁矩除法：净磁矩可能为零，需要按目标局域磁性区域定义归一化，先明确所研究的物理量。

4 鸢算科技集群实际结果#

600→800 eV的频率变化约1.07%，6→8 Å的变化约0.78%。400→600 eV的变化更明显，说明超精细耦合对数值参数的要求可能比总能直观变化更高。附件分别按2%和5%的预设比较标准核对这两个变化，它们不是实验频率误差保证。

磁矩翻转后，未经归一化的张量会变号；乘相同核g因子并除以M后，两种自旋方向给出相同的正频率。LDA条件得到23.196 cm，与GPAW官方小例中约23.2 cm的输出一致；它仍不同于实验约21 cm的谱线。泛函、PAW设置与数值参数的误差需要分开理解。

5 看实际自旋密度和张量#

左图是n↑−n↓的全电子自旋密度切片；右图为本例600 eV、6 Å条件的归一化质子耦合张量，显示单位为μeV。H没有冻结芯层，所以对本例使用exclude_core=True和默认调用得到相同结果；有芯层的元素不能据此推断芯层极化贡献总是为零。

附件保留原生平滑密度贡献、PAW修正贡献以及它们的和，完整两通道三维密度也保留在fields.npz中，没有只保存图上的切片。

6 导出CSV并重新画图#

source env.sh
"$YS_PY" - <<'PYCODE'
from pathlib import Path
import csv, json
v = json.loads(Path('results/VALIDATION.json').read_text())
with open('hyperfine.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['condition', 'frequency_MHz', 'wavelength_cm',
                'isotropic_A_eV', 'magnetic_moment'])
    for r in v['results']:
        w.writerow([r['name'], r['frequency_MHz'], r['wavelength_cm'],
                    r['isotropic_A_eV'], r['magnetic_moment']])
        print(r['name'], 'A(eV) =', r['normalized_proton_A_eV'])
PYCODE

只查看附件参考数据，不重新计算时：

bash setup.sh
bash plot.sh --reference
ls reference-plots

图片在reference-plots/。自己的完整结果在results/，附件原始数据在reference/；下载包不包含PAW数据集、GPW检查点或程序二进制。

7 环境与其他元素#

使用GPAW 26.7.0、Python 3.11.16，CPU构建采用GCC 13.3、HPC-X Open MPI 5.0.10rc2、MKL 2023.2、Libxc 7.0.0。setup.sh加载gpaw/26.7.0-8173m；本例通过legacy接口调用原生超精细工具，没有修改安装文件。

换成缺陷、自由基或其他核素时，需要明确同位素和核g因子，检查自旋态、局域磁矩、芯层处理、晶胞和基组。结果对核附近密度敏感，应针对目标张量重新检查截止能或网格，不能只看总能收敛。

方法、核g因子约定与氢小例：GPAW超精细耦合官方说明。