FourPhonon读取TDEP格式力常数

已有TDEP结果时，不能只把文件改名为FORCE_CONSTANTS_2ND。FourPhonon有专门的TDEP二阶常数读取入口。本例用已经计算的真实Si力常数演示这条输入路径；这里只做无损格式表达，不是重新拟合温度依赖力常数。

下载实际Si力常数、完整数值、CONTROL、三个结构/数据图及终端脚本。在已分配的CPU资源中计算，登录节点用于准备和提交。

python3 -m zipfile -e fourphonon-interfaces-wave321-*.zip .
cd fourphonon-interfaces
module load pymatgen/2026.9.24
bash run.sh

这一步读取已有输出并绘图，不重新计算DFT或输运。新的原生计算使用：

bash prepare.sh
sbatch cpu.sh
squeue -u "$USER"

CPU脚本采用当前集群的整节点资源申请约定；队列应选择当时可用、支持AVX512的兼容CPU资源。FourPhonon1.3-fix20260924固定源码5b814b3，oneAPI2024.2、Intel MPI2021.13与顺序MKL。这里不提供CPU节点性能结论。

1 两种输入各自用什么文件#

本例infile.forceconstant来自同一组16原子超胞的完整二阶常数。等长最短像分别分配权重，原胞原子索引和晶格平移写入文件；按超胞周期重新合并后，16×16×3×3全部分量与原数据一致。不是只复制一个3×3弹簧矩阵或只对齐一条Gamma频率。

head -n 12 inputs/infile.forceconstant
cat cases/Si-q5-TDEP-harmonic/CONTROL

2 完整CONTROL#

&allocations
 nelements=1, natoms=2, ngrid(:)=5 5 5
&end
&crystal
 lfactor=0.1
 lattvec(:,1)=0 2.715 2.715
 lattvec(:,2)=2.715 0 2.715
 lattvec(:,3)=2.715 2.715 0
 elements="Si", types=1 1
 positions(:,1)=0 0 0
 positions(:,2)=0.25 0.25 0.25
 masses(1)=28.085, gfactors(1)=0.00000000000000000e+00
 scell(:)=2 2 2
&end
&parameters
 T=300, scalebroad=1, maxiter=200, eps=1e-6
 num_sample_process_3ph=-1
 num_sample_process_3ph_phase_space=-1
 num_sample_process_4ph=-1
 num_sample_process_4ph_phase_space=-1
&end
&flags
 nonanalytic=.false., espresso=.false., tdep=.true.
 convergence=.false., four_phonon=.false.
 four_phonon_iteration=.false.
 isotopes=.false., autoisotopes=.false.
 nanowires=.false., onlyharmonic=.true.
&end

这里是Si非极性模型，nonanalytic=.false.。onlyharmonic=.true.在谐波/相空间阶段结束，不读取三四阶常数，也不输出完整的热导计算结果；因此本目录不需要放三、四阶文件。

3 在自己的计算资源中运行#

准备脚本已创建work/Si-q5-TDEP-harmonic/。单目录的原生命令如下，与批处理脚本选一种执行：

module load fourphonon/1.3-fix20260924-8173m
export OMP_NUM_THREADS=8 MKL_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1
cd work/Si-q5-TDEP-harmonic
mpirun -np 1 fourphonon
head BTE.omega

BTE.omega是角频率rad/ps，除以2π才是THz；不能直接当成THz画图。回到附件根目录，读取自己新算的结果：

module purge
module load pymatgen/2026.9.24
python verify.py --cases work --output new-results
python plot.py --data new-results --output new-plots

4 同一力常数的实际频率#

两条原生读取路径的单位常数精度略有区别，原始频率约相差0.00113%；包内按固定源码常量及其实际Fortran类型换算后，完整频率数组最大差约8.5E-9 rad/ps。这是输入读法的对应关系，不是通过修改力常数去拟合频率。

5 换成真实TDEP温度模型时#

真实温度相关结果需要先从相应温度的轨迹和受力拟合TDEP常数，再检查结构、原胞、质量、单位和实空间截断。极性材料还涉及Born电荷、介电张量与对应Ewald参数，并应使用匹配的裸短程常数，避免重复加入非解析修正。本Si输入没有演示这些极性步骤，不能将开关直接改成true就视为等价。

参考：固定版本TDEP输入说明。