GPAW最大重叠法激发态计算

本篇演示如何用最大重叠法（MOM）保持一组非基态轨道占据：先求CO基态，再把一个电子从占据轨道移到空轨道，让程序沿目标电子态收敛，最后检查能量、力和自旋密度。

普通自洽计算通常倾向于回到更低能的电子排布。MOM根据轨道重叠跟踪指定电子态，帮助避免这一过程。本例实际比较子空间投影和逐轨道重叠两种原生方式。

1 下载和运行#

下载CO最大重叠法算例包，保存为gpaw-mom-co.zip。从登录节点提交时，解压后使用sbatch job.slurm；已进入分配好的CPU资源时，可以执行：

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

任务使用1 MPI进程、1线程，依次计算7个条件：基态、两种MOM方法的中心结构，以及每种方法的两个真实键长位移。完整结果在results/VALIDATION.json，图片在plots/，各条件的原生输出和自己生成的GPW文件也保留在results/。

2 建立CO与基态输入#

模型为中性CO，C—O距离1.13 Å，6×6×6 Å非周期盒子。分子绕x方向旋转1°，便于固定这一小例的轨道方向；后面的位移计算始终保持同一晶胞和分子方向。

from ase import Atoms
from gpaw import GPAW, FermiDirac

atoms = Atoms('CO',
              positions=[[3, 3, 3-0.565], [3, 3, 3+0.565]],
              cell=[6, 6, 6], pbc=False)
atoms.rotate(1, 'x', center=[3, 3, 3])
atoms.calc = GPAW(
    legacy_gpaw=True, mode='lcao', basis='dzp',
    nbands=7, h=0.18, xc='PBE', spinpol=True,
    occupations=FermiDirac(0.01), symmetry='off',
    convergence={'energy': 1e-8, 'density': 1e-10,
                 'eigenstates': 1e-9, 'bands': 'all'},
    maxiter=300, txt='ground.txt')
ground_energy = atoms.get_potential_energy()
atoms.calc.write('ground.gpw', mode='all')

CO有14个物理电子；这里C、O的1s²各作为冻结芯层，显式价电子共10个。基态两个自旋通道各占据5条轨道。7条带包含本例要用到的空轨道，bands='all'使它们也参与本征态收敛检查。

3 指定目标占据并启动MOM#

本例把上自旋通道中索引4的电子移到索引5，下自旋通道保持原来的占据。Python轨道索引从0开始：

import numpy as np
from gpaw.mom import prepare_mom_calculation

target = np.array([
    [1, 1, 1, 1, 0, 1, 0],   # 上自旋：轨道4空出，轨道5占据
    [1, 1, 1, 1, 1, 0, 0],   # 下自旋
], dtype=float)
calc = GPAW('ground.gpw', legacy_gpaw=True, txt='MOM.txt')
atoms = calc.get_atoms()
atoms.calc = calc
prepare_mom_calculation(calc, atoms, target,
                        use_projections=True)
excited_energy = atoms.get_potential_energy()
forces = atoms.get_forces()

use_projections=True用目标占据子空间的投影；改为False时按逐轨道重叠匹配。附件让两种方式分别从同一个基态检查点开始，没有把第二种方式直接当作第一种方式的重复打印。

这里轨道索引适用于本例CO及其参数。换成自己的体系时，应先检查轨道能量、对称性和空间形状，再选择要移走和加入电子的位置，不能把4→5当作任意分子的固定设置。

4 实际能量与两种算法比较#

两种方式在本例落到了相同的电子态。相对能量是在固定结构上的垂直能量差。

**本例是混合自旋行列式，不是已经自旋纯化的单重激发态。**上下自旋电子数相等，只说明自旋投影对应的电子数差为0；它不能保证总自旋的平方对应纯单重态。因此不能直接把6.7398 eV标成实验单重态吸收峰。

图中两种算法的曲线重合，是实际计算结果；没有人为错开数值以制造差别。左图以中心结构的基态能量为统一零点，只有MOM分支做了这三个键长的计算。

5 两个真实位移检查受力#

附件从各自实际收敛的MOM检查点出发，沿键方向对称移动两原子，把键长改为1.127、1.133 Å；每个结构重新收敛电子态。中心差分为：

F(向外径向) ≈ −[E(1.133)−E(1.127)] / 0.006 Å。

两种MOM跟踪方式都完成了上述检查，差值小于本例0.02 eV/Å比较标准。较大的正向力说明当前1.13 Å几何在该电子态上有伸长趋势；它不是结构优化完成的标志。

6 看自旋密度与电子数#

总密度展示两个原子的分布；自旋密度为n↑−n↓。全局上下自旋电子数相同，并不要求每个空间位置的自旋密度都为0。

附件用LCAO的完整PAW重叠矩阵检查C*SCᵀ，独立得到每个自旋通道5个价电子。网格积分估计磁矩在本例约为2.1×10⁻⁵，属于有限网格积分结果；不要把这种微小差别误解成多出了一个电子或出现了磁矩翻转。两种数量都保留在结果表中。

7 导出能量、受力和磁矩分量#

source env.sh
"$YS_PY" - <<'PYCODE'
from pathlib import Path
import csv, json
v = json.loads(Path('results/VALIDATION.json').read_text())
with open('co-mom.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['condition', 'bond_A', 'energy_eV', 'radial_force_eV_A',
                'overlap_spin_difference', 'grid_estimated_moment'])
    for r in v['results']:
        w.writerow([r['name'], r['bond_A'], r['energy_eV'],
                    r['radial_force_eV_A'], r['overlap_net_spin'],
                    r['grid_estimated_moment']])
        print(r['name'], r['occupations'])
PYCODE

只重画附件中的实际数据时：

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

reference/保留完整两通道三维密度、占据、受力以及AO系数和重叠矩阵。轨道系数本身可以有相位或简并子空间的规范自由度，附件用完整密度和C*SCᵀ约束核对结果，不把系数逐分量的正负号当作物理差异。

8 环境与自己的模型#

使用GPAW 26.7.0的legacy LCAO接口、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。下载包不包含PAW数据集、GPW文件或程序二进制。

对于自己的分子，先确定电荷、自旋投影和目标轨道，检查空轨道数与表示范围，再设置MOM。做位移或优化时，应从目标态的收敛检查点接续，检查是否保持了同一电子态；不能只凭一次自洽收敛就认定目标态正确。弥散激发、强电荷转移或接近态交叉时，还需检查基组、盒子及适合的直接优化方法。

接口与方法：GPAW最大重叠法说明。