QE：电荷密度的平面平均与宏观平均

下载完整输入、九组工具输出、CSV 和绘图脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd qe-density-pdos-tools
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；后处理使用 NumPy、SciPy 和 Matplotlib。原始数据来自Si 电荷密度导出和Si 能带与投影的实际计算。本篇只做后处理，不重新运行 SCF、NSCF 或 projwfc。

1. 准备正确的输入文件#

average.x 读取 pp.x 的 filplot 中间文件。本例为 si-rho.dat，不是最终 Cube 文件。原 pp.x 操作为：

&INPUTPP
 prefix='si',
 outdir='./tmp',
 filplot='si-rho.dat',
 plot_num=0
/

plot_num=0 表示价电子密度；把文件名改成 LOCPOT 或 Cube 不会自动转换格式。分析静电势时也可采用对应的 pp.x 中间文件，但纵轴的物理量与单位会不同。

2. 六行输入分别代表什么#

mkdir average-new
cd average-new
cp ../inputs/si-rho.dat .
cat > average.in <<'EOF'
1
si-rho.dat
1.0
60
3
2.0
EOF
average.x < average.in > native.out
head -n 5 avg.dat
cd ..

依次为：文件数量、文件名、该文件权重、输出点数、晶胞方向编号、宏观平均窗口长度。窗口长度的单位是原子长度单位 Bohr。本例沿第三个原胞方向处理，原 FFT 网格为 60×60×60。

avg.dat 有三列：

程序在当前目录写固定文件名 avg.dat，不同条件应分别建目录。

3. 提高输出点数意味着插值#

mkdir average-120-new
cd average-120-new
cp ../inputs/si-rho.dat .
average.x < ../inputs/average-a3-120.in > native.out
wc -l avg.dat
cd ..

输入仅把点数改为 120，使用 Fourier 插值；它没有产生更多 SCF 密度数据，也不会改善原始电子结构的收敛。

本例还实际测试了方向 1、2、3，以及请求 20 点的情况。20 小于该方向原网格 60 点时，程序会提升到 60，不能把它当作直接降采样命令。

4. 斜晶胞中坐标要这样理解#

Si 的 FCC 原胞并非正交立方盒。程序按分数网格对另外两个方向求平均，并用所选原胞矢量长度标记横坐标；不能直接把它写成笛卡尔 z，也不能当成相邻平均平面之间的垂直距离。

因此资料包图片用“第三个原胞方向的分数坐标”表示横轴。对通常沿真空方向正交的表面模型，平面平均静电势更容易直接对应表面法向；斜晶胞则应先核对具体几何和单位。

5. 检查守恒并重画#

head -n 8 analysis-read/planar-density.csv
bash run.sh
bash read.sh work-new analysis-new
bash plot.sh analysis-new

本例独立核对了三维密度网格平均、平面平均、Fourier 插值和周期窗口，价电子密度积分约为 8 个电子。积分需使用 Bohr³ 晶胞体积，不能混用 ASE 的 ų 几何单位。

窗口取 2 Bohr 是工具演示，不是所有材料的推荐宏观平均长度。用于表面势或异质结时，应结合层周期、真空平台和目标物理量选择窗口。当前版本的输入与坐标约定可对照 QE 7.6 average 源码。