VASP cRPA计算有效相互作用U与J

DFT+U需要输入U，但U不是只由元素名称决定的常数。cRPA从指定局域轨道出发，保留其他电子的屏蔽作用，计算这个低能模型中的有效相互作用。本例在鸢算科技集群实际计算SrVO₃，展示三个V t₂g轨道之间的U、轨道间排斥和交换J怎样取得、怎样阅读。

1. 先明确计算的对象#

模型是五原子的立方SrVO₃，晶格常数3.84652 Å；采用固定晶胞、固定理想原子位置，未做晶胞优化。高对称结构的原子力为0.000000 eV/Å。非磁PBE、PAW64的Sr_sv/V_sv/O，共41个价电子；ENCUT=550 eV，基态电子收敛EDIFF=1E-6，4×4×4 Gamma网格，Fermi展宽0.1 eV。

目标子空间是V的dxy、dxz、dyz三个Wannier轨道。这里对应第21—23条能带；与相邻带的最小同k点间隔为1.432、0.316 eV。换材料、换赝势或换磁性设置，都要重新核对能带编号，不能直接照抄排除列表。

使用vasp/6.6.1-8173m-wannier90，Intel oneAPI 2024.2构建，链接Wannier90 3.1.0。这是CPU程序，示例采用4个MPI进程；不要把它当成GPU加速数据。

2. 工作流：先算基态，再补空带，最后算屏蔽#

下载附件后，在calculations下每个目录准备自己合法的POTCAR，顺序为Sr_sv、V_sv、O。回到附件根目录执行：

sbatch run.sh

脚本依次执行以下步骤，自动传递本次生成的波函数和导数文件：

第二步使用ALGO=Exact、NELM=1和LOPTICS=.TRUE.。这里的单步是读取已有自洽解后做精确对角化，不是从随机初态只做一次普通SCF。

cRPA步骤的关键输入为：

ALGO = CRPA
NBANDS = 96
ENCUTGW = 150
NCSHMEM = 1
NTARGET_STATES = 1 2 3
LWRITE_WANPROJ = .TRUE.
NUM_WANN = 3
WANNIER90_WIN = "
num_bands = 96
num_iter = 200
exclude_bands = 1-20,24-96
begin projections
V:dxy;dxz;dyz
end projections
"

NTARGET_STATES的1、2、3是三个Wannier轨道的编号，不是原始DFT能带的1—3。屏蔽中排除的是目标子空间内部的贡献。这里关注静态ω=0，UIJKL.1文件头中两个频率分量都应为0。

3. 实际得到什么#

U取直接相互作用矩阵的对角平均；U′取非对角平均；J取交换矩阵的非对角平均。原始四指标张量保留了更多信息，不要把它压成一个数字后就忘记轨道定义。

图中左侧对照未屏蔽与屏蔽后的相互作用；中间展示最终组的直接相互作用矩阵；右侧展示三个设置下的U、J变化。三个Wannier中心均位于V原子等价位置，展宽约1.6672 Ų。

这些对照只检查了两个总能带数和两个响应截断能，尚不能证明4³k网格已收敛。实际研究还要增加k网格、空带、ENCUTGW，核对Wannier子空间和泛函，并保持各组结构、赝势及局域轨道定义一致。

4. 如何提取自己的结果#

计算结束后在附件根目录运行：

module load ase/3.29.0
python extract_crpa.py

生成CRPA-DATA.json，包含三组U、U′、J、直接/交换矩阵、Wannier中心与展宽。参考数据在RESULTS.json，图片在figures；附件不含赝势和原始波函数。

阅读OUTCAR时，先看bare Hubbard U/J，再看screened Hubbard U/J。UIJKL.1是本次VASP版本输出的静态屏蔽四指标张量，VIJKL是未屏蔽张量；其他版本或频率设置的文件编号可能不同，应同时核对文件头。

5. U和J拿到后，能直接填进DFT+U吗#

不能只把本例的U搬到任何含V体系。这里的三个t₂g Wannier轨道，与某些DFT+U实现使用的完整原子d投影空间不同；屏蔽排除规则和双计数约定也会改变参数含义。

对于Dudarev形式，输入的有效参数通常按U_eff=U−J组织，但这个减法本身不能解决目标子空间不一致的问题。应先确定后续模型采用什么局域轨道、哪些相互作用项，再决定如何使用cRPA矩阵。

本例可作为“从电子结构到模型相互作用”的流程参考，不是SrVO₃高精度定值，也不包含DMFT求解、关联谱或温度相图。

参考与下载#

方法入口：VASP cRPA官方示例、cRPA形式与目标子空间。本文表格和图片来自上述集群实算。

下载输入、自动顺序脚本、提取器和参考数据