QE：孤立原子的全电子径向计算

下载三组完整输入、实际径向结果、CSV与绘图脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd qe-atomic-radial
module load ase/3.29.0 qe/7.6-8173m
bash read.sh
bash plot.sh

本例使用 QE 7.6 的 ld1.x，CPU 构建为 Intel oneAPI 2024.2 / Intel MPI 2021.13 与 MKL；后处理使用 NumPy、SciPy 和 Matplotlib。ld1 处理孤立原子的径向问题，与周期体系的 pw.x 平面波计算采用不同的数据和输入方式。

1. 选择全电子计算#

iswitch=1 选择全电子原子计算。本例 H 使用自旋极化 LSDA，He 使用闭壳层 LDA，泛函均为 PZ；rel=0 为非相对论径向方程，不依赖外部赝势文件。

H 的完整输入：

&input
 title='H-LSDA',
 zed=1,
 config='1s1',
 iswitch=1,
 rel=0,
 lsd=1,
 isic=0,
 latt=0,
 dft='PZ',
 xmin=-8.0,
 dx=0.008,
 rmax=40,
 tr2=1d-12,
 prefix='atom',
 file_charge='density.dat',
 noscf=.false.
/

zed 是核电荷，config 是电子组态，lsd=1 开启自旋极化；isic=0、latt=0 表示没有附加自相互作用或 Latter 修正。径向长度使用 Bohr。

tr2 控制 ld1 的内部自洽误差，不应直接当作 VASP 的 EDIFF 或 pw.x 的 conv_thr。对数径向网格由 xmin、dx、rmax 决定，最后一个实际网格点可能略大于请求的 rmax。

2. 提交或在已分配CPU资源内运行#

资料包提供 CPU 作业模板：

bash submit.sh
squeue -u "$USER"

模板遵循本集群 cpu-8173m 分区的 56 核整节点申请规则，实际每个原子计算使用 1 个 MPI 进程、1 线程，三组顺序执行；本例没有做并行性能测试。

已有相应 CPU 资源时，也可直接调用单个输入：

mkdir H-new
cd H-new
ld1.x < ../inputs/H-LSDA.in > native.out
grep -E 'final scf error|Etot|Total magnetization' native.out
cd ..

只查看现成结果时，第一节的 read/plot 即可，不需要重新提交计算。

3. He输入及实际结果#

He 使用 zed=2、config='1s2'、lsd=0，其完整输入在包内：

cat inputs/He-LDA.in
grep -E 'Etot|Ekin|Encl|Eh|Exc' \
  reference-results/He-LDA/native.out

本例实际总能：

H 的总磁矩为 1 μB。这些属于对应近似泛函的径向原子结果，不是精确原子能，也不能直接与使用不同冻结芯约定的赝势总能相减。

4. 径向电荷怎样积分#

density.dat 第一列为 r，后面的列为径向电子分布；H 有两个自旋通道，He 为一个合计通道。它已经包含球面角向测度：

nᵣ(r)=4πr²ρ(r)，电子数 N=∫nᵣ(r)dr。

因此不能再把该列乘一次 4πr²后积分。需要三维电子密度时，反过来除以 4πr²，单位为 electron/Bohr³。

head -n 5 reference-results/H-LSDA/density.dat
head -n 6 analysis-read/atomic-summary.csv

本例积分分别恢复 H 的 1 个和 He 的 2 个电子。

5. 能量分解与新结果读取#

Ekin 为动能，Encl 为核吸引，Eh 为 Hartree，Exc 为交换关联能；它们之和对应 Etot。不要把单个 Kohn–Sham 轨道本征值当作这个多项相加后的总能。

bash read.sh work-new analysis-new
bash plot.sh analysis-new
cat analysis-new/atomic-energies.csv

先等待提交的作业完成再读取 work-new。对照氢样解析轨道和波函数单位，继续阅读径向波函数与氢样参考。参数定义见 QE 官方 ld1 文档。