GROMACS蛋白SAXS与SANS实操

从MD轨迹计算散射曲线，除了选择SAXS或SANS，还需要处理周期边界、确认原子元素、统一q与坐标单位，并对每一帧的强度求平均。本文使用鸢算科技集群既有的真实泛素轨迹，给出可以直接执行的完整操作。

模型为76残基、1231原子的蛋白子集，包含C378、H629、N105、O118、S1。取0、20、40、60、80、100ps六帧，分别计算全原子与76个Cα的SAXS/SANS；Cα只是几何对照，不代表完整蛋白的散射权重。

1. 下载并执行#

下载真实蛋白子集轨迹、输入、独立Debye工具和结果：

mkdir protein-scattering
cd protein-scattering
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构建信息见软件环境。

当前2026.4的原生 gmx scattering 在元素识别和多帧累计处理上有版本限制。本篇采用配套Python Debye工具，显式读取元素、加入自项、逐帧独立计算后平均；它是本教程的辅助工具，不是另一个官方GROMACS命令。不要直接把不同输出约定的曲线混合使用。

2. 先把蛋白变完整，再导出坐标#

run.sh 先执行以下原生操作；脚本会建立 work/coordinates，这些命令均已实际运行：

module purge
module load gromacs/2026.4-p100-threadmpi
cd work/coordinates
printf '0\n' | gmx trjconv -s ../../inputs/protein.tpr -f ../../inputs/protein.xtc -o whole.xtc -pbc whole -dt 20 -e 100
printf '0\n' | gmx trjconv -s ../../inputs/protein.tpr -f whole.xtc -o structure.pdb -dump 0
printf '0\n' | gmx traj -s ../../inputs/protein.tpr -f whole.xtc -ox all-positions.xvg -fp -xvg none
cd ../..

本包TPR与XTC都只含同一批1231个蛋白原子，组0就是蛋白整体；自己的含水体系不能照抄这个组号，应先提取一致的蛋白拓扑和轨迹，或者调整脚本的选择与元数据顺序。-pbc whole 把同一分子复原，后面的Debye求和使用完整蛋白的普通空间距离，不再对每个原子对做盒子最小镜像。

-dt 20 的单位是ps；本例保存六帧。gmx traj -fp 导出高精度坐标，单位是nm。PDB文件提供原子名称和元素列，PDB中的坐标通常为Å；脚本不使用PDB坐标代替高精度nm坐标，只读取顺序对应的元素与名称。

3. SAXS/SANS的完整Debye求和#

对每一帧计算：

I(q) = Σ_i f_i(q)² + 2Σ_(i<j) f_i(q)f_j(q) sin(qr_ij)/(qr_ij)。

第一项是自散射项，第二项包括对称原子对；q=0按sinc极限1处理。SAXS使用Cromer–Mann原子形状因子，SANS使用相干中子散射长度；本例不加入非相干背景。

q网格是0到2nm⁻¹、间隔0.1nm⁻¹，共21点。Cromer–Mann表中的长度系数以Å为基准，计算形状因子前先把q乘0.1转换成Å⁻¹；计算 sin(qr) 时仍用nm⁻¹×nm。不要只改横轴单位而不改数值。

辅助工具的核心实际实现为：

import numpy as np

def amplitudes(q,kind,species,factors):
 q=np.asarray(q)
 if kind=='sans':return np.array([[factors[e][1] for e in species]]*len(q))
 assert kind=='saxs'
 s2=(q*.1/(4*np.pi))**2
 return np.array([factors[e][-1]+sum(factors[e][2+j]*np.exp(-factors[e][6+j]*s2) for j in range(4)) for e in species]).T

def coherent_intensity(q,positions,amplitudes):
 q=np.asarray(q);amp=np.asarray(amplitudes);i,j=np.triu_indices(len(positions),1);dist=np.linalg.norm(positions[i]-positions[j],axis=1);out=(amp*amp).sum(axis=1)
 for start in range(0,len(i),32768):
  sl=slice(start,start+32768);out+=2*np.sum(amp[:,i[sl]]*amp[:,j[sl]]*np.sinc(q[:,None]*dist[None,sl]/np.pi),axis=1)
 return out

元素参数随包提供，来源为本集群GROMACS2026.4的 scatteringfactors.dat，保留许可证。真实同位素替换应准备明确的同位素组成，不能把所有H随意改成D后解释成实验交换比例。

4. 运行、导出与看图#

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 重现图表。

module purge
module load ase/3.29.0
python - <<'PYCODE'
import csv
with open('work/analysis/all-saxs.csv') as f:
    for r in csv.DictReader(f):
        if float(r['q_nm_inverse']) in [0.0, 1.0, 2.0]:
            print(r['q_nm_inverse'], r['mean_coherent_intensity'], r['normalized_I_over_I0'])
PYCODE

全原子SAXS的I(0)=21077556.8006 electron²，SANS的I(0)=3347990.9177 fm²；这是孤立蛋白坐标模型的相干求和值，不是仪器记录的实验绝对强度。q=2nm⁻¹时，全原子SAXS与SANS的I/I(0)分别约0.131126和0.120374。

5. 迁移到自己的体系#

先确认原子顺序、元素和分子完整性，再修改取帧时间与q网格。配套脚本对元素集和原子数做了本例专用检查；更换结构时需要同步修改这些明确的检查条件，不能直接给新体系沿用蛋白名称和参数。

与实验比较还需要适当的溶剂与排除体积、水合层、浓度/粒子间结构因子、同位素对比度、背景和仪器处理。本例六帧没有覆盖这些因素，也不能证明蛋白系综已经收敛。平均、归一化和Cα对照的区别见散射结果解读。其他相关操作可继续阅读NMR距离约束、FRET几何和偶极与电荷分析。

原生轨迹参数：trjconv、traj；原生散射选项与版本：scattering。