QE：线性响应计算Hubbard U

下载完整输入、实际日志、响应矩阵、能级CSV和绘图脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd qe-hubbard-U-response
module load ase/3.29.0 qe/7.6-8173m
bash read.sh
bash plot.sh

本例使用鸢算科技集群的 QE 7.6 CPU 版本，构建工具为 Intel oneAPI 2024.2、Intel MPI 2021.13 和 MKL。材料为四原子菱方 LiCoO₂；固定原子位置和晶胞，采用非磁 PBEsol、USPP、50/400 Ry 截断能、20 条能带、2×2×2 电子 k 网格和 1×1×1 响应 q 网格。小网格适合先熟悉流程；正式使用需要加密并检查 U 及所需性质。

1. 先做含Hubbard投影子的基态#

顺序是：准备结构和赝势 → pw.x 基态 → hp.x 响应 → 读取 U → 新的含 U 计算。hp.x 需要前一步保留的 tmp/LiCoO2.save 等文件，仅复制 scf.out 不够。

完整 inputs/scf.in：

&control
 calculation='scf', restart_mode='from_scratch',
 prefix='LiCoO2', pseudo_dir='../pseudo', outdir='./tmp',
 verbosity='high', tprnfor=.true., tstress=.true.
/
&system
 ibrav=5, celldm(1)=9.3705, celldm(4)=0.83874,
 nat=4, ntyp=3, ecutwfc=50, ecutrho=400, nbnd=20
/
&electrons
 conv_thr=1d-10, mixing_beta=0.4, electron_maxstep=150
/
ATOMIC_SPECIES
Co 59.0 Co.pbesol-spn-rrkjus_psl.0.3.1.UPF
O  16.0 O.pbesol-n-rrkjus_psl.0.1.UPF
Li 7.0 Li.pbesol-s-rrkjus_psl.0.2.1.UPF
ATOMIC_POSITIONS crystal
Co 0.0 0.0 0.0
O 0.2604885 0.2604885 0.2604885
O 0.7395115 0.7395115 0.7395115
Li 0.5 0.5 0.5
K_POINTS automatic
2 2 2 0 0 0
HUBBARD ortho-atomic
U Co-3d 0.0001

HUBBARD ortho-atomic 定义正交化原子投影子。极小的 U Co-3d 0.0001 用于启用对应投影子路线；它是 0.0001 eV，不能称为严格 U=0。这里的 Co-3d、赝势轨道和后续 hp.x 的相关子空间必须一致。

2. 准备三种赝势#

资料包未打包赝势，提供官方文件网址和 SHA256。运行：

python3 setup-pseudo.py
ls pseudo
cat pseudopotentials.json

本例使用 Co.pbesol-spn-rrkjus_psl.0.3.1.UPF、O.pbesol-n-rrkjus_psl.0.1.UPF、Li.pbesol-s-rrkjus_psl.0.2.1.UPF。更换赝势、交换关联泛函、投影子或结构后，应重新计算 U，不能直接把本例数值复制到其他材料。

3. 提交自动流程#

bash submit.sh
squeue -u "$USER"

模板申请 cpu-8173m 的 56 核整节点，短例使用 1 MPI 进程和 1 线程。脚本依次执行基态、响应和代入实际 U 后的新基态。已获分配的 CPU 资源内，也可执行 bash run.sh；不要同时运行两份。

希望理解每一步时，可在分配好的 CPU 资源内手动执行：

mkdir -p manual/pseudo manual/baseline
cp pseudo/*.UPF manual/pseudo/
cd manual/baseline
pw.x < ../../inputs/scf.in > scf.out
grep 'convergence has been achieved' scf.out
ls tmp/LiCoO2.save
hp.x < ../../inputs/hp.in > hp.out
cd ../..

两个程序使用同一工作目录、prefix 和 outdir；不要在运行 hp.x 之前清理 tmp。

4. 响应输入与完成信号#

完整 inputs/hp.in：

&inputhp
 prefix='LiCoO2', outdir='./tmp',
 nq1=1, nq2=1, nq3=1,
 conv_thr_chi=1d-7, iverbosity=2,
 niter_max=100, max_seconds=550
/

conv_thr_chi 控制响应自洽，不等于 pw.x 的电子收敛阈值。nq1/nq2/nq3 是响应 q 网格，不能拿电子 k 网格替代。检查：

grep -E 'iter #|chi:|CONVERGENCE HAS BEEN REACHED' work-new/baseline/hp.out
cat work-new/baseline/LiCoO2.Hubbard_parameters.dat

本次实际响应经过 13 次迭代，末次残差为 5.40×10⁻⁸，得到 U=7.775865 eV。

5. 从响应矩阵理解U#

本例只有一个相关 Co 位点，可直接看标量形式：U=χ₀⁻¹−χ⁻¹。实际文件包含：

chi0 matrix:       -0.261000
chi matrix:        -0.086153
chi0 inverse:      -3.831416
chi inverse:      -11.607281
Hubbard matrix:     7.775865

用打印出来的 χ₀、χ 重新求倒数，约为 7.775839 eV；差异来自响应矩阵仅保留六位小数。取实际 Hubbard matrix 的值用于下一步，而不是用四舍五入后的响应重复求 U。

多相关位点是矩阵求逆，不能逐个元素取倒数；含 V、磁性或多个 Hubbard manifold 时，应按对应设置处理。本例不等价于这些更复杂情况。

下一篇：将实际 U 代回计算并比较能级。