Fourthorder计算与收集真实受力

下载Si真实结构、受力、力常数、8张图和终端脚本，保存为 si-fourphonon.zip。所有本系列文章使用同一个包：

python3 -m zipfile -e si-fourphonon.zip .
cd fourthorder-fourphonon-silicon-wave288
bash preview.sh
ls preview/figures

包中 reference/ 保存本次实算参考；work/ 是您重新生成的文件，results/data/ 是提取后的CSV。受力来自真实VASP计算；下载包不含POTCAR、波函数、完整电子输出或第三方程序。

1 两条使用路径#

想学习力常数重构，可直接使用包内真实受力，执行 bash prepare.sh、bash reconstruct.sh。想为自己的体系计算新的受力，则按本篇准备DFT任务，不能把Si的受力套到别的元素或结构。

本例三个H各48个位移，另有三阶8个，共152组逻辑受力。逐个核对晶胞、原子顺序、几何和电子参数后，33组完全相同的构型复用了既有结果，实际完成119个独立静态DFT；不是把152个文件都算成152次独立计算。全部使用鸢算科技集群P100与V100 16G的单GPU任务。

2 静态INCAR#

SYSTEM=Si16-Fourthorder-H020-01
ENCUT=450
PREC=Accurate
EDIFF=1E-6
NBANDS=48
KPAR=1
ISMEAR=0
SIGMA=0.05
ALGO=Normal
NELM=160
LREAL=.FALSE.
LASPH=.TRUE.
LWAVE=.FALSE.
LCHARG=.FALSE.
IBRION=-1
NSW=0
ISYM=0

这里 IBRION=-1、NSW=0 表示仅求固定构型的电子态与受力，不允许把位移结构重新优化回平衡位置。ISYM=0保留受力计算的完整自由度，LREAL=.FALSE.有利于力的一致性；电子收敛为1E-6。EDIFFG不是这个静态任务的停止条件。所有位移应使用相同PAW、ENCUT、k点、NBANDS及展宽。

KPOINTS：

Si16 electronic mesh
0
Gamma
3 3 3
0 0 0

这是16原子超胞的电子k点网格，与后面FourPhonon的声子q网格是两回事。

3 准备48个受力目录并提交#

以下以H040为例；下载包已有真实受力时，不必再提交这些DFT任务。

bash prepare-dft.sh
ls dft-forces/01
cat dft-forces/01/INCAR
sbatch forces.sh
squeue -u "$USER"

forces.sh是48个数组子任务，每个申请P100一张卡、一个MPI进程、一个CPU核。POTCAR仅在授权VASP作业启动时从集群势库复制。更换元素时必须同时更换结构、合适的PAW及电子参数。

#!/bin/bash
#SBATCH -J si-fourthorder-force
#SBATCH -p gpu-8p100-hs
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -c 1
#SBATCH --gres=gpu:p100-16g:1
#SBATCH --mem=10G
#SBATCH -t 00:15:00
#SBATCH --array=1-48
set -eo pipefail
source /etc/profile
module purge
module load vasp/6.6.1-p100-hdf5
export OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1
printf -v id '%02d' "$SLURM_ARRAY_TASK_ID"
cd "$SLURM_SUBMIT_DIR/dft-forces/$id"
test ! -e vasprun.xml
cp /public/software/vasp/potcar/version64/PBE/Si/POTCAR POTCAR
ulimit -s unlimited
mpirun --bind-to none -np 1 vasp_std > vasp.log

本例另用模块 vasp/6.6.1-v100-hdf5 在V100 16G执行了相应单卡任务；切换平台时同时修改分区、GPU型号和模块，不能只改一个名称。本例没有用V100 32G。

4 检查完成并收集#

tail -n 20 dft-forces/01/OUTCAR
grep 'EDIFF is reached' dft-forces/01/OUTCAR
module load pymatgen/2026.9.24
python collect_forces.py
wc -l fourthorder-xml.list
head -n 3 fourthorder-xml.list

应得到48行。配套收集脚本检查每个构型只有一次静态计算、电子循环未到NELM、受力形状为16×3且数值有限，并检查总力漂移。本例最大电子迭代14步，119个独立任务均保留了原生受力。

收集到真实 vasprun.xml 后，回到与sow相同H的工具环境，将列表传给reap。不要用按修改时间排列的文件列表，也不要缺几个结构就填零。具体见重构与位移量检查。

5 下载包里的XML是什么#

为了不分发完整电子输出，配套重建脚本从实际数值受力生成只含forces数组的最小XML载体，再交给原生reap。这些载体不是新的DFT，也不冒充完整VASP输出。您自己提交的新DFT则使用实际生成的完整 vasprun.xml 收集受力。