QE电荷密度与数据导出

下载完整输入、赝势、真实结果和绘图脚本。解压后，inputs/是可修改输入，reference/是本次真实参考输出，复跑结果写入 work/。下面的计算命令在已分配的相应资源内运行。

1 电荷密度从哪来#

先完成电子结构流程，用pp.x读取同一prefix/outdir中的电荷密度，再导出cube文件。本篇示例是Si的价电子密度，不是包含全部冻结芯电子的全电子密度。

bash run.sh si-electronic
cat inputs/si-electronic/density.in
ls work/si-electronic/si-density.cube

2 pp.x完整输入#

&INPUTPP prefix='si', outdir='./tmp', filplot='si-rho.dat', plot_num=0 /
&PLOT nfile=1, filepp(1)='si-rho.dat', weight(1)=1, iflag=3, output_format=6, fileout='si-density.cube' /

plot_num=0选择电荷密度；iflag=3导出三维；output_format=6写cube。filepp对应前一段filplot，fileout是最终文件名。若只手工改其中一处名字，后处理可能读错前置文件。

单独在计算目录执行：

module load qe/7.6-8173m
cd work/si-electronic
pp.x -in density.in > density.in.out
head -n 8 si-density.cube

3 用工具核对网格与电子数#

from ase.io.cube import read_cube
from ase.units import Bohr
with open('si-density.cube') as f:
    cube = read_cube(f)
rho = cube['data']
volume_bohr3 = cube['atoms'].get_volume() / Bohr**3
print('grid:', rho.shape)
print('valence electrons:', rho.mean() * volume_bohr3)

本次网格60×60×60，密度单位electron/Bohr³，积分为7.999968，接近Si2的8个价电子。ASE把几何长度转换成Å，并不自动改变密度数组的物理单位；积分时漏掉Bohr³转换会得到错误电子数。

图中密度切片为第三个原胞分数坐标0.5平面，横纵是前两个分数坐标；Si原胞并非正交立方盒，不能把这张图直接标成笛卡尔x–y距离。cube可用VESTA等三维工具查看，等值面、切片和配色应注明实际数值。

4 重绘与替换自己的数据#

cd ../..
module load calorine/4.0
python extract.py charge
python plot.py
ls qe-pdos-density.png

附件plot.py读取reference中的原始输出。客户使用自己的结果时把相应数据路径改为work；保持单位、原胞定义和目标物理量一致。差分电荷需要相同晶胞、网格、坐标与电子数口径，不能将不同网格数组直接相减。