py4vasp 激子权重与空间密度

从Si真实BSE本征矢和条件密度读取跃迁权重、实际空穴位置与54原子超胞。下面按准备、读取、核对与画图的顺序操作。

1. 下载数据，先提取再画图#

下载响应与激子终端实操包，解压后先看现成数据：

unzip py4vasp-response-practical-wave290-*.zip
cd py4vasp-response-practical-wave290
bash preview.sh

再完整运行读取、原生工具与绘图：

bash run.sh
ls results/data
ls plots/figures

run.sh依次执行prepare.sh、extract.sh、tools.sh、plot.sh。结果写入results/data/，原生数据报告和视图矢量在tools-output/，七张PNG在plots/figures/；输出目录受保护，重复运行请换一个新解压目录。随包的figures/可直接查看。

后处理用py4vasp/0.11.3的Python3.11环境，静态图用pymatgen/2026.9.24环境的NumPy2.3.5、Matplotlib3.11.2，只需已分配资源内1个CPU，不重新运行VASP。carriers/保存实际计算导出的自有数值与字段配方；运行时在work/重建最小HDF读取载体，不是原始完整VASP输出。这里全部选定数组保留完整分辨率，无插值、无抽样，不包含PAW、波函数或程序。

inputs/保留本例完整INCAR、POSCAR、KPOINTS。真实GPU计算使用VASP6.6.1/HDF5、NVHPC25.3、CUDA12.8、P100或V10016G；BSE使用VASP6.6.1 CPU/HDF5、Intel oneAPI2024.2、MPI、HDF51.14.6。提取工具不替代VASP求解。

2. 先有真实GW前置文件，再做BSE#

本例沿用Si两原子原胞真实G₀W₀计算保存的64带WAVECAR、WAVEDER与屏蔽文件，原GW的16条准粒子带覆盖本次所需4条占据带和8条空带。这里没有重新计算基态或GW，新增BSE本征矢与固定空穴条件下的电子空间密度。

mkdir Si-BSE
cp inputs/Si-BSE/{INCAR,POSCAR,KPOINTS} Si-BSE/
cp gw/{WAVECAR,WAVEDER,POTCAR} Si-BSE/
cp gw/W0*.tmp gw/WFULL0*.tmp Si-BSE/

gw/应换成自己的、与本例设置匹配的前置计算目录；不能拼接不同晶胞、k点、带数或PAW的波函数与屏蔽文件。关键设置：

ALGO = BSE
IBSE = 2
ANTIRES = 0
NBANDS = 64
NBANDSO = 4
NBANDSV = 8
NBSEEIG = 2
BSEHOLE = 0.11 0.13 0.17
LH5 = .FALSE.
LCHARGH5 = .TRUE.

本安装版本的条件密度使用IBSE=2的新精确对角化驱动和Tamm–Dancoff近似；NBSEEIG选择输出的最低两个态。这里输入是普通WAVECAR与W*.tmp，所以全局LH5保持false，单独用LCHARGH5输出密度HDF；全局开启LH5会改为寻找对应HDF前置文件。详见官方IBSE和BSEHOLE坐标说明。

3. 读取本征矢，先归一化再分析贡献#

module load py4vasp/0.11.3
python - <<'PYCODE'
import numpy as np
import py4vasp
c = py4vasp.Calculation.from_path('work/Si-BSE')
d = c.exciton.eigenvector.read()
A = d['eigenvectors']
print('本征矢:', A.shape)
print('原生模平方:', np.sum(abs(A)**2, axis=1))
A = A / np.sqrt(np.sum(abs(A)**2, axis=1))[:,None]
index = d['bse_index'][0]
weight = abs(A[:, index])**2
print('跃迁映射:', index.shape)
print('每态权重和:', weight.sum(axis=(1,2,3)))
np.save('normalized-transition-weights.npy', weight)
PYCODE

本例为27个k点×8条空带×4条占据带=864跃迁，本征矢形状(2,864)。py4vasp0.11.3读出的原生向量模平方是 729=27²，不能直接把模平方当已归一的概率。上面先对每态归一化，权重和才是1。BSEFATBAND最后两列的复振幅已经采用归一化约定，不能和HDF原生振幅直接比数值；其带宽列又使用原生振幅的绝对值。

API的bse_index、first_valence_band、first_conduction_band已经转为Python零起点。本文CSV改回VASP习惯，价带1–4、导带5–12，且保留实际跃迁索引。

两个输出态能量均约 2.96001912 eV，实际光学强度非零；后面也有接近零强度的暗态，不能将“第几个本征态”直接当“第几个亮峰”。前两态处于简并子空间，单态图形依赖所选本征基底；比较简并态应考虑整个子空间的组合，不以单态形状得出独立物理结论。

4. 固定空穴实际落在哪#

请求的BSEHOLE=(0.11,0.13,0.17)是原胞分数坐标。实际程序使用原胞网格上的约(0.083333,0.125000,0.166667)，即(2,3,4)/24，原OUTCAR写出该位置。不能拿请求坐标直接标注图中的空穴。

module load py4vasp/0.11.3
python - <<'PYCODE'
import py4vasp
c = py4vasp.Calculation.from_path('work/Si-BSE')
x = c.exciton.density.read()
print(x['charge'].shape)
print(x['structure']['lattice_vectors'])
print('密度超胞原子数:', len(x['structure']['positions']))
PYCODE

密度形状是 (2,72,72,72)，对应实际3×3×3超胞、54个Si原子。这里结构要从exciton.density.read()['structure']取；不能把普通c.structure的2原子原胞坐标盒直接套到密度上。空穴在超胞分数坐标中还要除3。

5. 条件空间密度如何画#

图取超胞c=4/72的完整平面，横纵为斜晶格的超胞分数a、b坐标，红色十字为实际空穴。原生网格平均值约4.06224、2.37926，并非一电子归一化密度。为比较空间形状，图中对每态显式除以mean(raw)×Vsupercell，使显示空间权重体积分为1，单位Å⁻³；这项展示归一化不能解释成实际激发电子数。原始未缩放网格完整保留在NPZ。

head -5 results/data/BSE-optical-transitions.csv
head -5 results/data/BSE-transition-weights.csv
ls tools-output/Si-BSE-exciton.*.txt

tools.sh还实际调用exciton.density.to_view('1', isolevel=1.0)生成原生等值面视图对象；浏览器三维显示需要自己的对应前端环境，静态截面不依赖该环境。这里密度采用原生HDF输出，没有额外宣称CHG文本跨格式一致性。小带窗、稀k网格的Si例用于说明流程，不作为收敛激子束缚能或光谱的材料数据。

提交自己的CPU BSE任务#

将匹配的前置文件准备好后，在自己的Si-BSE目录保存submit.sh。CPU分区按56核整节点申请；该HDF5版本为MPI构建；这个小模型模板采用1个MPI进程，整节点申请规则保持。更大任务应按算法与输入兼容性调整MPI规模，设置OMP_NUM_THREADS不会启用OpenMP。

#!/bin/bash
#SBATCH --partition=cpu-8173m
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=56
#SBATCH --time=00:30:00
set -eo pipefail
source /etc/profile
module purge
module load vasp/6.6.1-8173m-hdf5
ulimit -s unlimited
export OMP_STACKSIZE=256M
export OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1
mpirun -np 1 vasp_std > vasp.log

sbatch submit.sh
squeue -u "$USER"

分区可用时才提交新的CPU计算。资料包的bash run.sh只做后处理，不提交该任务。

完整计算输入#

PAW和前置文件按正文自行提供；下列输入也在下载包中，不会由后处理脚本自动提交。

Si-BSE#

SYSTEM=Si-new-BSE-exciton-observables
ENCUT=550
PREC=Accurate
EDIFF=1E-6
EDIFFG=-0.02
IBRION=-1
NSW=0
ISIF=2
ISMEAR=0
SIGMA=0.05
LREAL=.FALSE.
LASPH=.TRUE.
LWAVE=.FALSE.
LCHARG=.TRUE.
NELM=1
ALGO=BSE
ISYM=2
NBANDS=64
ENCUTGW=150
ISTART=1
IBSE=2
NBANDSO=4
NBANDSV=8
CSHIFT=.1
NBSEEIG=2
BSEHOLE=.11 .13 .17
ANTIRES=0
LH5=.FALSE.
LCHARGH5=.TRUE.

Si
 1.0000000000000000
     0.0000000000000000    2.7344629562984917    2.7344629562984917
     2.7344629562984917    0.0000000000000000    2.7344629562984917
     2.7344629562984917    2.7344629562984917    0.0000000000000000
 Si 
   2
Direct
  0.0000000000000000  0.0000000000000000 -0.0000000000000000
  0.2500000000000000  0.2500000000000000  0.2500000000000000

Gamma mesh
0
Gamma
3 3 3
0 0 0

继续阅读#

• py4vasp 完整弹性模量
• py4vasp Berry极化与位移响应
• py4vasp 原子局域磁矩
• py4vasp MD速度与温度
• 全部quantity与文件前提