GROMACS散射平均与结果解读
六帧真实曲线、全原子与Cα、SAXS/SANS单位、自项和归一化的操作对照。
阅读本篇目录
得到一条散射曲线后，先检查它代表哪些原子、如何平均、是否包含自项，以及单位和归一化。本文使用同一套真实泛素六帧数据，演示逐帧强度、平均曲线、Cα与全原子的差别。
1. 取得实际输入与结果#
mkdir scattering-average
cd scattering-average
curl -fL https://yuansuan.top/assets/downloads/full-function/gromacs-scattering-wave90-38bdcdb34f85.zip -o examples.zip
unzip examples.zip
bash run.sh
这一步只处理既有轨迹，不启动新MD。小包使用1个CPU核即可；普通登录节点可做本例小规模处理，自己的大轨迹请提交CPU作业。重复处理请用新的解压目录。脚本不覆盖已有 work。
轨迹工具模块为 gromacs/2026.4-p100-threadmpi，本例不调用GPU内核；Python分析使用 ase/3.29.0。GROMACS构建信息见软件环境。
坐标提取、完整求和和版本适用范围见SAXS/SANS实操。debye.py 是配套辅助工具，采用显式PDB元素和完整相干Debye公式；不把当前原生 gmx scattering 的多帧输出直接作为这里的平均值。
2. 先平均强度，再按定义归一化#
在本例相同原子与元素条件下，每帧先独立计算I(q)，再按六帧等权平均。图上的粗线为 mean[I(q)]/mean[I(0)]，阴影是各帧强度的样本标准差除以平均I(0)。不是先平均坐标，也不是把平均结构的散射当成系综平均。
本例每帧原子组成相同，I(0)也相同，所以逐帧归一化后平均与“先平均后归一化”在这里一致；变化选择、组成或权重时，不能自动沿用这一等价关系。
终端重现：
module purge
module load ase/3.29.0
python debye.py work
python check_results.py work
python export_structures.py work
python plot_results.py work
坐标和元素元数据在 work/coordinates/，CSV、逐帧NPZ及XYZ在 work/analysis/，图片在 work/figures/。包内 reference/ 是本文实际结果；未重做轨迹提取时，可把以上 work 改为 reference 重现图表。
python - <<'PYCODE'
import numpy as np
a = np.load('work/analysis/all-saxs-frames.npz')
intensity = a['intensity']
print('时间/ps:', a['times_ps'])
print('每帧I(0):', intensity[:, 0])
print('平均I(q)/平均I(0):', intensity.mean(axis=0)/intensity[:, 0].mean())
print('帧间标准差:', intensity.std(axis=0, ddof=1))
PYCODE
六帧来自同一条100ps轨迹，相邻帧未必统计独立；帧间标准差不是独立样本标准误，也不能从它判断所有慢构象变化已经收敛。
3. 为什么不能把Cα曲线当作完整蛋白#
Cα对照保留76个碳位置，适合检查链的几何、原子选择和距离单位。完整蛋白还包含侧链、N/O/S及大量H，不同元素的X射线形状因子与中子散射长度不同；SANS中H的相干散射长度还可以为负。
| 原子选择 | 类型 | I(0) | I(2nm⁻¹)/I(0) |
|---|---|---|---|
| 1231全原子 | SAXS | 21077556.80063674 | 0.131125898 |
| 1231全原子 | SANS | 3347990.91770193 | 0.120374187 |
| 76个Cα几何对照 | SAXS | 207880.55409665 | 0.146971760 |
| 76个Cα几何对照 | SANS | 255121.96921600 | 0.148267439 |
SAXS的原始强度单位是electron²，SANS为fm²，两列I(0)不能跨种类直接比较大小来讨论“哪一种信号更强”。I/I(0)是无量纲形状对照，也不能据此忽略溶剂背景、实验对比度与归一化方式。
4. q和坐标单位要成对检查#
q[nm⁻¹] = 10 × q[Å⁻¹]。例如2nm⁻¹等于0.2Å⁻¹。同一结构若把坐标从nm乘10变成Å，几何核中的q就应除10，才能保持qr不变。形状因子表仍要按其自身的Å系数解释。
包内 check_results.py 实际检查了两碳解析/方向积分、自项、q=0、平移/旋转不变性和nm/Å等价，也核对每个CSV与六帧原始强度的平均关系。q=0的完整相干极限应为 (Σf_i(0))²，不能把只含 i<j 的交叉对项直接当作同一个I(0)。
5. 结果文件怎么用于后续作图#
head -n 5 work/analysis/all-saxs.csv
head -n 5 work/analysis/all-sans.csv
CSV五列依次为q[nm⁻¹]、平均相干强度、帧间标准差、I/I(0)和自项。NPZ保存每一帧的完整强度，可以重新选窗口或检查异常帧；XYZ是从nm转换成Å的真实坐标副本，保留原子序号。
需要实验可比的散射时，应在合理构象采样基础上加入适当的溶剂/排除体积、实验背景和同位素模型，并选择相应专业工具。这里的孤立蛋白求和用于学习可解释的数据流程，不替代这些后续步骤。