GPAW局域轨道与有效模型
从水分子实际LCAO计算得到局域轨道，比较3、20、23轨道模型，给出终端流程、矩阵、谱线和局域轨道投影图。
阅读本篇目录
LCAO得到的哈密顿量可以按原子分块，分别对角化，得到更适合分析局域能级和构建有效模型的轨道。本篇用一个实际水分子计算，解释如何得到局域轨道，以及最小模型、扩展模型与完整模型之间的关系。
这里使用GPAW原生LocalOrbitals，它与跨原胞构建的Wannier函数是不同的分析方法。
1 运行与结果目录#
下载水分子LCAO操作包，保存为gpaw-lcao-workflows.zip。普通作业在解压目录执行sbatch job.slurm。已分配好CPU资源时依次执行：
python3 -m zipfile -e gpaw-lcao-workflows.zip .
cd gpaw-lcao-workflows
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh
程序完成六个小分子条件后，读取results/A-own/state.gpw，在results/local/进行轨道后处理。读取和变换过程中不重新做SCF。其他条件分别用于私有基组和重叠误差修正教程。
2 准备LCAO结果#
本例水分子采用固定坐标、PBE、DZP、0.22 Å网格，1个O和2个H共有23个原子轨道基函数。使用nbands='nao'保留完整基组空间，并用calc.write('state.gpw', mode='all')保存系数。
在已经完成计算的calc上，主要接口为：
from gpaw.lcao.local_orbitals import LocalOrbitals
los = LocalOrbitals(calc)
los.subdiagonalize('OH', groupby='energy')
print(los.groups)
los.take_model(minimal=True)
H = los.get_hamiltonian() # eV
S = los.get_overlap()
对于从GPW恢复的旧后端LCAO，重叠矩阵可能延后初始化。附件local.py使用原生位置初始化恢复重叠矩阵，并保留原波函数系数；直接运行附件即可，无需手工改检查点。
3 分块变换做了什么#
每个原子对应一个基函数块：O为13维，两个H各5维。程序在各个原子的H、S子块中解广义本征值问题，得到变换矩阵U。完整变换后的矩阵为U† H U和U† S U。
每个原子内部的重叠块成为单位阵，哈密顿块成为对角阵。不同原子之间仍可能有重叠和耦合，所以后续一般仍应求解广义本征值问题，不能直接把整个S丢掉。
4 三类实际模型#
| 模型 | 本例维度 | 实际含义 |
|---|---|---|
| 自动最小模型 | 3 | 每个原子选择局域能级靠近费米能级的轨道，并加入原生静态杂化修正 |
| 扩展模型 | 20 | 将与最小空间有足够耦合的轨道组加入模型 |
| 完整模型 | 23 | 保留全部局域轨道，恢复完整LCAO谱 |
from scipy.linalg import eigh
los.take_model(minimal=False, cutoff=1e-3)
e_extended = eigh(los.get_hamiltonian(), los.get_overlap(),
eigvals_only=True)
los.take_model(indices=list(range(23)), minimal=True)
e_full = eigh(los.get_hamiltonian(), los.get_overlap(),
eigvals_only=True)
cutoff在这个扩展判断中对应耦合矩阵元的阈值，单位eV。3、20、23是这组结构和设置的实际结果；换结构、基组或费米能级后，自动选出的维度可以改变。最小模型没有保证恢复全部能级，不能把三条谱线当作完整电子结构。
查看自己的模型：
python3 -c "import json; v=json.load(open('results/local/VALIDATION.json')); print([(r['name'],r['dimension'],r['indices']) for r in v['models']]); print(v['blocks'])"
5 局域轨道形状与正交化#
附件还把实际局域轨道映射到实空间网格。下图按每个原子的轨道块对赝轨道模平方求和，再沿一个方向积分，显示轨道空间的大致分布。
这是赝轨道权重，不是全电子电荷密度或电子数；轨道可有不同相位，同一简并空间内也可以旋转，不能只看单个系数的正负判断两次结果不同。
当前安装版本的原生take_model(ortho=True)接口存在初始化限制，本例不用该标志。附件另列出一个显式的对称Löwdin辅助变换：对所选空间的S分解得到S^(-1/2)，同时变换H和S。它是示范脚本中的数学辅助步骤，不作为原生标志已成功运行的结果。
source env.sh
python3 -c "import numpy as np; d=np.load('results/local/fields.npz'); print('original H:',d['original_H_eV'].shape); print('orthogonal S:',d['explicit_Lowdin_S'])"
ls -lh plots/lcao-local-models.png plots/lcao-local-orbital-projectors.png
6 用于自己的体系#
先完成LCAO自洽，保留足够带和完整检查点；选择分析元素或明确的基函数块，再检查模型窗口、原子间重叠、扩展阈值和完整谱。将得到的模型用于传输或其他计算之前，还需验证目标能区的误差。
环境为GPAW26.7.0旧后端，GCC13.3、OpenMPI5.0.10rc2、MKL2023.2和Python3.11.16。方法与原生接口见GPAW Local Orbitals教程。