FourPhonon 自动同位素参数与热导

计算同位素散射时，需要分别决定平均原子质量、质量方差参数 g，以及是否启用同位素散射。本篇用四组 Si 实算，演示自动读取和手动设置之间的关系。

1. 下载并先看结果#

下载力常数、输入、实际输出与绘图脚本，保存为 example.zip：

python3 -m zipfile -e example.zip .
cd fourphonon-auto-isotopes

source ./env.sh
(module load ase/3.29.0; ase info --files inputs/POSCAR-primitive)
cat isotope-table.json
bash show-reference.sh
bash plot-results.sh
ls -lh plots/fourphonon-auto-isotope-mass-scattering.png

力常数沿用 Si真实受力与力常数流程：2 原子原胞、16 原子超胞，三阶采用 H0.01 Å，四阶采用 H0.04 Å 的结果。下载包已包含二、三、四阶力常数和参考结构，可以直接进行下面的四组 CPU 输运计算。

使用 FourPhonon 1.3-fix20260924，固定提交 5b814b337e03e72907a15a2c953723171fba7102，Intel oneAPI 2024.2 Fortran / Intel MPI 2021.13 / MKL 顺序库，CPU 构建。分析使用已有 ASE 3.29.0、NumPy 和 SciPy；ASE 用于读取 POSCAR。示例保持单 MPI 进程、单线程；生产计算请在分配的兼容 CPU 资源内运行。

2. 自动读取会覆盖哪些参数#

自动组的关键输入为：

&allocations
 nelements=1, natoms=2, ngrid(:)=3 3 3
&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
 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=.false.
 convergence=.false., four_phonon=.true., four_phonon_iteration=.false.
 isotopes=.true., autoisotopes=.true.
 nanowires=.false., onlyharmonic=.false.
&end

autoisotopes=.true. 会按 elements 从内置数据库重新计算 masses 和 gfactors 两项。 因此这里写入的 28.085 和 0 会被内存中的自动结果覆盖，原 CONTROL 文件仍保持原文字。

这个固定版本的 Si 内置表为：

按这张表得到平均质量 28.0855099896 u，质量方差 g=0.00020091202549504095。g 是无量纲量：先按丰度求平均质量，再求 Σ cᵢ(1−mᵢ/平均质量)²，其中 cᵢ 为归一化比例。

本系列早期手动自然 Si 例子采用了另一套明确列出的质量/丰度值。比较两组时，应同时核对质量表、平均质量、g 和开关，不能只看到“自然 Si”就假定输入完全相同。

3. 实际运行自动组#

source ./env.sh
mkdir manual
cp inputs/FORCE_CONSTANTS_* manual/
cp reference-results/native/automatic/CONTROL manual/
(cd manual && mpirun -np 1 fourphonon)
cat manual/T300K/BTE.kappa_tensor

fourphonon 无需额外参数，在当前目录读取 CONTROL 和三份力常数。BTE.kappa_tensor 保留完整 3×3 热导张量；本例为立方模型，三个对角分量相同。

4. 手动控制质量与方差#

手动匹配组在 crystal 中设：

masses(1)=28.0855099896
gfactors(1)=0.00020091202549504095

并在 flags 中使用：

isotopes=.true.
autoisotopes=.false.

手动指定质量或 g 时，把 autoisotopes 明确关闭。 当前版本默认允许自动读取，若仍保留自动开关，写入的手动值会被覆盖。

完整匹配输入已在 reference-results/native/manual-matched/CONTROL。这个组与自动组的频率、三/四声子散射、同位素散射及热导均一致。

另外一组 double-variance 保持相同平均质量，仅令 g 加倍。实际结果显示，所有同位素散射率均加倍，而频率、群速度和三/四声子散射保持一致。它是人为控制方差的对照，不对应本文另外指定的一种实际同位素丰度组合。

5. 关闭同位素散射后，自动质量仍可能生效#

第四组采用：

isotopes=.false.
autoisotopes=.true.

这会关闭同位素散射，但自动平均质量仍会被读入，所以频率与前面三组相同。若希望保留自己在 CONTROL 中写入的质量，同时关闭同位素散射，应将两个开关都明确设为 false。

四组均为 300 K、3³ 网格、三加四声子 RTA：

这里全部保持相同有效平均质量和相同力常数，可以清楚地区分额外散射造成的变化。

6. 从输出提取散射与寿命#

根目录 BTE.w_isotopic 为同位素散射；温度目录中的 BTE.w_3ph、BTE.w_4ph 和 BTE.w 分别给出三声子、四声子及合并散射率。当前组合权重已包含在原生通道中，满足 Γ总=Γ3+Γ4+Γ同位素。

第一列为角频率 rad/ps，第二列为散射率 ps⁻¹。转换 THz 时除以 2π；按本程序约定，RTA 寿命为 1/Γ总，单位 ps。零散射模没有在表中伪造一个很大的有限寿命。

配套工具按实际输出的分支/不可约 q 点顺序整理 mode-scattering.csv，并用 Bose 热容、完整 q 映射、速度及寿命在 SI 单位下独立重建完整 3×3 RTA 张量。更详细的通道说明见 散射率与寿命提取。

7. 重算四组并生成图表#

sbatch job.slurm

或在已经获得的兼容 CPU 资源内：

bash run.sh
bash plot-results.sh
cat auto-isotope-summary.csv
head -n 8 mode-scattering.csv

整套会建立四个工作目录，依次进行原生计算、提取和绘图。已有 work/ 时停止；重新解压到另一个目录可以从头重算。

本包用于固定 Si 力常数与 3³ 网格的开关对照，配套分析也针对这套模型。研究新的材料和同位素组成时，应重新准备结构、力常数、质量/丰度与收敛检查。3³ 的这些数值用于方法参考；较细网格的实际操作见 温度网格与求解选项。

内置质量与丰度可对照 该固定提交的官方数据表。