Phonopy 中子与X射线散射

下载本节数值输入、脚本与绘图程序，在自己的CPU工作目录运行：

unzip phonopy-dsf-dos-nac-files-wave249.zip -d phonopy-practice
cd phonopy-practice
bash setup.sh
bash run.sh

setup.sh联网创建本目录个人.env，安装固定的Linux x86_64预编译环境，并下载、校验ASE3.29.0的纯Python包到extras/。不改公共模块或用户全局Python。环境包含Phonopy4.8.3、Python3.11.18、NumPy2.4.6、symfc1.7.3、pypolymlp0.22.1；C/C++运行库为libgcc/libstdcxx16.2.0，绘图使用集群phonopy/4.8.3模块。流程单线程运行，不需要GPU，也不启动新的DFT。

source/保存原生命令结果，analysis/、dsf/、xray/、debye-windows/保存提取数据，figures/保存四张高清图。为保护已有结果，重复练习请换新目录解压。下面手动操作都另建工作目录。

1 本教程计算什么#

这里计算谐近似下、相干散射的单声子产生项：中子探针比较0 K与300 K，X射线探针演示300 K。它不是含分辨率、多声子、无序和仪器几何的完整实验谱，也不包含声子湮灭的全部通道。

Si力常数来自集群实际16原子VASP DFPT，设置见文件与输入说明。脚本读取原始数值，用12×12×12声子网格计算Debye–Waller因子；电子k网格与声子q网格是两回事。

2 设置Q、温度和中子散射长度#

cat dsf.py
head -n 8 dsf/onephonon-data.csv

核心API：

import phonopy
p = phonopy.load("inputs/Si-geometry.yaml",
                 force_constants_filename="inputs/Si-FORCE_CONSTANTS",
                 symmetrize_fc=False)
p.run_mesh([12,12,12], is_gamma_center=True,
           with_eigenvectors=True, is_mesh_symmetry=False)
Q = [[1.1,0,0], [1.2,0,0], [1.3,0,0]]
ds = p.run_dynamic_structure_factor(
    Q, 300, scattering_lengths={"Si": 4.1491e-5}, freq_min=0.1)
print(ds.frequencies)
print(ds.dynamic_structure_factors)

这里Q使用声子原胞倒易基矢下的约化坐标，不能直接当成立方常规胞的(hkl)。程序再分解为第一布里渊区内q与倒格矢G。CSV同时保存Q、q、分支、THz、meV和强度，便于自己画能量轴。

NIST表中的天然Si相干束缚散射长度为4.1491 fm，传给当前Phonopy接口前换算成4.1491×10⁻⁵ Å。直接填4.1491会让强度大10¹⁰倍；能量单位换算与散射长度换算不能混用。

freq_min=0.1避开Γ附近零频率。真实材料如果存在显著虚频，应先检查结构与谐模型，不能靠设置阈值就当作稳定。

3 运行与结果图#

完整bash run.sh已经运行两种探针。单独练习中子部分时，在新解压目录先完成环境准备，再运行：

export PYTHONPATH="$PWD/extras"
.env/bin/python dsf.py
head -n 8 dsf/onephonon-data.csv

谱线位置由声子频率决定；强弱还受本征矢、原子质量、相位干涉、散射长度或形状因子、Bose占据和Debye–Waller因子影响。因此能带中的一条分支可能在某个Q处非常弱，不等于程序漏算。

本例同时保存0 K结果。声子产生项含n+1，所以0 K仍可有强度。简并子空间中的本征矢可以旋转，比较程序结果时应比较简并模式总强度，而不是要求每个任意本征矢的强度分别相同。

4 X射线：形状因子和物理单位#

cat inputs/Si-WK1995.json
cat xray.py
head -n 8 xray/Si-form-factors.csv
head -n 8 xray/Xray-onephonon.csv

本例采用ESRF DABAX提供的Waasmaier–Kirfel 1995中性Si系数，用五个高斯项加常数计算f₀(s)，s=|Q|/(4π)。f₀(0)约13.998917，接近Si的14个电子。

原始f₀是以电子数表示的无量纲形状因子。若要与中子分支一样以m²/J表示，本例传入的是经典电子半径rₑ[Å]×f₀(s)；没有乘rₑ的数组另行保留为无量纲输入下的API结果，不贴成绝对散射截面。

这采用中性独立原子的非共振近似，未加入反常色散、偏振和具体仪器几何因子；归一化取N·k′/k=1。图中两种探针的相对强弱反映此模型及单位约定，不能直接当作实验探测器计数。

需要改材料时，先准备该结构的可靠力常数，再替换元素散射参数、Q、温度和网格；不要保留Si参数只改标题。接口与归一化定义见Phonopy动态结构因子文档，本包实际运行版本为4.8.3。