GPAW局部自旋约束实操

局部电荷约束控制一个区域中的总电子数；局部自旋约束控制该区域中上、下自旋电子数之差。两者不是同一个参数。本例让H2的总磁矩保持为零，同时让第一个H附近的高斯分区具有+0.03或−0.03的局部自旋差，观察两个区域相反的极化。

1 下载与计算#

下载局部自旋约束实操包，保存为gpaw-spin.zip。从登录节点提交可在解压后使用sbatch job.slurm；已经分配到CPU资源时，直接执行：

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

例程使用GPAW 26.7.0 CPU构建，1 MPI进程、4线程。results/保存实际计算、网格与数值结果，plots/保存本次图片。setup.sh会展开并核对原图使用的参考网格；下载包中的压缩是字节无损的。

请使用完整附件，保留spin_compat.py。 它适配本例所用26.7.0 FD自旋势接口，只作用于当前示例计算器，不修改公共安装目录。不要仅复制run.py到其他版本环境中。本例不计算约束态力、跃迁率或自旋激发谱。

2 结构、总自旋与局部目标#

两个H分别位于(4,4,4.63)和(4,4,5.37) Å，键长0.740 Å，盒为8×8×10 Å，三个方向不周期。几何保持固定。

run.py中的主要设置为：

atoms = Atoms('H2',
              positions=[[4, 4, 4.63], [4, 4, 5.37]],
              cell=[8, 8, 10], pbc=False,
              magmoms=[0.03, -0.03])
base = GPAW(legacy_gpaw=True, mode='fd', h=0.25,
            xc='PBE', spinpol=True,
            occupations=FermiDirac(0.01, fixmagmom=True),
            convergence={'energy': 1e-7, 'density': 1e-6,
                         'eigenstates': 1e-9},
            maxiter=200, txt='SCF.txt')
compat.attach(base)  # run.py已按文件加载随包的兼容层
c = CDFT(base, atoms,
         spin_regions=[[0]], spins=[0.03],
         spin_coefs=[-1.5],
         mu={'H': 0.5}, Rc={'H': 20.0},
         method='BFGS',
         minimizer_options={'gtol': 1e-4, 'maxiter': 50},
         tol=1e-4, compute_forces=False, txt='CDFT.txt')
atoms.calc = c
energy = atoms.get_potential_energy()

fixmagmom=True配合初始磁矩之和为零，固定全体系的自旋差；spin_regions=[[0]]指定只约束第0号原子的分区。spins=[0.03]是分区电子数差，不是“这个原子的量子自旋S=0.03”。在本例约定下，它对应约0.03 μB的部分自旋磁矩。

spin_coefs是外层优化器的初始势系数，单位eV；它不是最终自旋值。mu为高斯宽度、Rc为高斯支持半径，单位Å。本小盒使用20 Å使高斯支持覆盖整个计算盒；这个数值不是其他材料的通用推荐值。

目标与宽度在PLAN.json的conditions中，几何和GPAW设置在run.py中。更换体系时，要一起检查边界、网格、自旋状态和区域定义，并在新示例目录中运行。

3 本集群实际结果#

全体系始终有2个电子，总磁矩接近零。相同宽度的正、负目标给出相同能量，实际差远小于1 meV；改变分区宽度后，即便目标相同，定义的约束态和所需势系数也会改变。

前两幅图是实际赝自旋密度n↑−n↓的切片，红、蓝代表相反符号；黑点是实际原子位置。右图是第一个H的高斯分区权重。

本例计算的局部量为高斯加权的赝自旋密度积分，再按GPAW cDFT的原子分区约定添加PAW增广修正。它不是Bader自旋、形式氧化态，也不能仅凭该局部数值认定实验自旋多重度。

表中能量使用c.Edft，表示约束后密度的DFT部分。c.Ecdft还带有所施加势的贡献；不要把这两个字段混用，也不要直接把不同分区的能量差当作完整实验激发能。

4 终端导出表格#

source env.sh
"$YS_PY" - <<'PYCODE'
import json, csv
from pathlib import Path
v = json.loads(Path('results/VALIDATION.json').read_text())
with open('local-spin.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['target_e', 'width_A', 'region_spin_e',
                'DFT_energy_eV', 'coefficient_eV', 'total_spin_e'])
    for r in v['results']:
        w.writerow([r['spin_target_e'], r['gaussian_width_A'],
                    r['region_spin_e'], r['KS_constrained_energy_eV'],
                    r['lagrange_eV'], r['total_magnetization_e']])
        print(r['name'], r['region_spin_e'],
              r['total_magnetization_e'])
PYCODE

5 自己读取网格积分和重画原图#

下面直接读取某一条件的完整网格，算一次区域积分：

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
from pathlib import Path
f = sorted(Path('results').glob('*-FIELDS.npz'))[0]
d = np.load(f)
pseudo = (d['weight'] * d['pseudo_spin_density_bohr3']).sum()
pseudo *= float(d['grid_volume_bohr3'])
aug = d['augmentation_e'][0] - d['augmentation_e'][1]
print('pseudo spin integral =', pseudo)
print('PAW augmentation =', aug)
print('region spin =', pseudo + aug)
PYCODE
bash plot.sh --reference
ls reference-plots/*.png

归档图在reference-plots/，本次计算图在plots/。再次执行相同完整输入会读取现有结果；更换目标或模型后请保留旧结果，使用新目录，便于对照。

构建环境：GCC 13.3、HPC-X/OpenMPI 5、MKL 2023.2、Libxc 7.0.0，GPAW 26.7.0、PAW数据26.8.0，随包兼容层仅适配所注明CPU FD接口。区域约定见GPAW约束DFT说明；局部电荷的操作见约束电荷实操。