py4vasp 选带电荷密度

从Si第4、5带的实际LPARD输出读取分波电荷，说明归一化、空带及原生负值。先用实际数据完成提取和绘图，再按相同顺序准备自己的计算。

1. 先用现成数据走通流程#

下载真实空间数据与终端实操包，解压进入资料目录：

unzip py4vasp-realspace-practical-wave289-*.zip
cd py4vasp-realspace-practical-wave289
bash preview.sh

先在 preview/figures/ 看5张静态图，再执行：

bash run.sh
ls results/data
ls tools-output

run.sh依次准备读取载体、提取数值、调用py4vasp原生绘图、生成静态PNG。四步也可分开运行 bash prepare.sh、bash extract.sh、bash tools.sh、bash plot.sh。脚本保护已存在的输出，重复操作请换一个新解压目录。后处理使用已分配资源内1个CPU；不会重新运行VASP，也不需要GPU。

后处理模块为 py4vasp/0.11.3（Python3.11独立环境）；静态图用 pymatgen/2026.9.24 环境的NumPy2.3.5、Matplotlib3.11.2。HTML用浏览器打开，通过CDN加载Plotly；离线查看可直接用PNG。

包内 carriers/*.npz 是本次实际计算导出的自有数值。prepare.sh在work/重建只包含指定数组的最小HDF读取载体，并写明来源；它不是原始完整VASP输出，也不能用来读取未提供的物性。包中不含PAW、波函数或VASP程序。reference/保留独立提取的数据，inputs/提供各例完整INCAR、POSCAR、KPOINTS。

2. 先选择带，再生成分波电荷#

本例沿用Si两原子原胞真实G₀W₀计算保存的WAVECAR，64条带、4个不可约k点、550 eV。这里只从已存在波函数生成第4、5带的选带电荷，不重新做SCF或GW。该后处理实际使用VASP6.6.1 GPU/HDF5（NVHPC25.3、CUDA12.8），单张P100；py4vasp读取本身只用CPU。

LPARD = .TRUE.
LPARDH5 = .TRUE.
IBAND = 4 5
NBMOD = 0
LSEPB = .TRUE.
LSEPK = .FALSE.

保留原波函数的晶胞、PAW、k点与对称性设置。LSEPB使两条带分别保存，LSEPK=.FALSE.汇总k点；所以读取时k点标号为0。LPARDH5的文件与版本前提见VASP说明，选带语义见LPARD。

准备自己的工作目录时可按以下顺序操作，PAW按原计算自行提供：

mkdir partial-charge
cp scf/POSCAR scf/KPOINTS partial-charge/
cp scf/WAVECAR scf/POTCAR partial-charge/
cp inputs/Si-partial/INCAR partial-charge/
# 按自己的波函数带数和目标能带检查INCAR后，再通过Slurm提交

3. 明确指定第4、5带#

第1节运行后，在资料根目录执行：

module load py4vasp/0.11.3
python - <<'PYCODE'
import numpy as np
import py4vasp
c = py4vasp.Calculation.from_path('work/Si-partial')
d = c.partial_density.read()
V = abs(np.linalg.det(d['structure']['lattice_vectors']))
print('bands:', d['bands'], 'kpoints:', d['kpoints'])
for band in [4, 5]:
    raw = c.partial_density.to_numpy(band=band)
    rho = raw / V
    print(band, rho.shape, '积分/e:', rho.mean()*V)
    np.save('Si-band%d-density.npy' % band, rho)
PYCODE

分别保存带4、5时，默认band=0代表的汇总数组并不存在，必须显式写band=4、band=5。本例网格60³、体积40.89273224 ų，两带积分分别 2.0000000087、2.0000000210 e。

4. 归一化和负值怎样理解#

本次VASP6.6.1的HDF网格与同一波函数生成的原生PARCHG逐点比较，二者使用相同数值约定，均有mean(raw)≈2。因此本例密度为 raw/V，单位e/ų；再除60³会错误地把积分缩成约9.26×10⁻⁶ e。处理其他版本/其他格式时，先用输出守恒关系确认，不要机械叠加网格因子。

选中自旋简并的空带后也得到2的归一化，这不是实际导带占据或激发载流子数。不能拿这张图直接计算激发态电子总数。本例原生带分解表示有负值；第4带范围约−0.05190～0.16925 e/ų，第5带约−0.02892～0.12082 e/ų，图和数据均保留原符号。它不是纯平面波概率密度或完整全电子电荷，不能先裁成零再归一化。

图取c=0.25，a、b为晶格分数坐标，两幅共用色标；比较的是相同截面与表示。积分表在results/data/partial-integrals.csv，三维数值在Si-band-partial-density.npz。换自己的体系时，先依据能带与占据确定目标IBAND，再比较多个截面或等值面，不能只凭一张切片认定整个轨道形状。

完整计算输入#

下列为对应原生计算的输入，PAW和前置电荷/波函数请按正文说明自行准备。后处理包默认不运行这些计算。

SYSTEM=Si-true-wave-band-density
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=.FALSE.
NELM=1
ALGO=Normal
ISYM=2
NBANDS=64
ISTART=1
LPARD=.TRUE.
LPARDH5=.TRUE.
IBAND=4 5
NBMOD=0
LSEPB=.TRUE.
LSEPK=.FALSE.

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 表面功函数
• py4vasp 空间屏蔽与诱导电流
• 全部quantity与文件前提