Fourthorder生成四阶位移
从Si原胞生成16原子超胞，理解最近邻截断、48个位移及位移量的单位。
阅读本篇目录
下载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 先分清两个程序#
Fourthorder根据晶体对称性生成有限位移结构，再从DFT受力重构四阶力常数；FourPhonon读取二、三、四阶力常数，计算散射与热导。它们不替代VASP/QE的电子结构计算。上游Fourthorder和FourPhonon分别提供这两段流程。
本例晶格常数5.43 Å、金刚石Si原胞2原子，2×2×2超胞16原子，四阶相互作用先取最近邻，目的是走通操作；正式研究还需检查超胞、截断、电子k网格和声子q网格。
2 准备原胞POSCAR#
Si
1.0000000000000000
0.0000000000000000 2.7149999999999999 2.7149999999999999
2.7149999999999999 0.0000000000000000 2.7149999999999999
2.7149999999999999 2.7149999999999999 0.0000000000000000
Si
2
Direct
0.0000000000000000 0.0000000000000000 -0.0000000000000000
0.2500000000000000 0.2500000000000000 0.2500000000000000
坐标是分数坐标。这个胞不是含8原子的传统立方胞，不能在超胞大小相同的情况下混用两种原胞及其力常数。
3 生成位移#
先用配套工作流建立独立目录，再调用原生工具：
bash prepare.sh
bash reconstruct.sh
ls work/fourthorder/H040/4TH.POSCAR.*
head -n 10 work/fourthorder/H040/4TH.SPOSCAR
reconstruct.sh会实际调用三组Fourthorder sow/reap和Thirdorder sow/reap；这里读取包内已有真实受力，不运行新DFT。
如果您从自己的POSCAR独立开始，原生入口如下：
module load fourthorder/0.1-fix20260618
mkdir my-fourthorder
cp reference/reference/POSCAR-primitive my-fourthorder/POSCAR
cd my-fourthorder
Fourthorder_vasp.py sow 2 2 2 -1
ls 4TH.POSCAR.*
最后的 -1 表示按第一近邻确定截断；正数截断的单位是 nm，不是Å。4TH.SPOSCAR是未位移超胞；4TH.POSCAR.01至.48是48个逻辑位移结构。不要只凭文件名推断原子序号，原胞、超胞和受力的原子排列必须对应。
图取本例H=0.04 Å的第一个实际构型；左侧箭头放大14倍，右侧保留真实长度。三次位移可能落在同一原子/方向，因此单个原子的总移动不一定恰好等于H。
4 默认位移量与本例的三组比较#
安装版 Fourthorder_common.py 的H单位为nm，默认 1e-3 nm = 0.01 Å。本例另比较0.02、0.04 Å；对应脚本值分别为0.002、0.004 nm。sow和reap必须使用同一H，否则四阶常数会按H³产生明显错误。
配套脚本只在 work/programs/H010、H020、H040 的私人副本修改H，不改共享安装，也不把第三方源码随包分发。每组48个逻辑位移，共144组；其中部分完全相同的物理构型可严格核对后复用，不能随意用对称结构的受力替换。
5 软件环境#
模块 fourthorder/0.1-fix20260618：Python 3.11、NumPy 1.26.4、SciPy 1.13.1、GCC 13.3、spglib 2.7，采用官方Python3资源和距离修正。版本后缀为集群构建标识，不是上游新发行号。其编译扩展应由配套Python运行，不要把它导入另一套NumPy 2环境。下一步见计算与收集受力。