GROMACS 蛋白螺旋几何与局部取向

下载本组输入、真实参考结果和绘图脚本，解压后进入资料包目录。inputs 是可以修改的输入，reference 是鸢算科技集群已经计算得到的参考结果。每次运行在 work 下建立新目录，路径写入 LAST-RUN.txt。计算请在已经分配的资源中执行；整组可用 sbatch run.sbatch 提交。

数据与一键操作#

资料包包含泛素实际轨迹的纯蛋白子集：1231 个原子，0–100 ps，共 51 帧；原始蛋白建模流程见蛋白质教程。本例研究第 23–34 残基的 12 个 Cα。

bash run.sh protein

这一步运行 helixorient、helix 与配套的 FRET 几何演示。原始模拟的水与离子已从下载子集中去掉，不必为几何分析下载整个溶剂盒。

1. 准备自己轨迹的蛋白子集#

若使用完整体系，先截取相同的蛋白原子，并保持 TPR、GRO、XTC 的原子顺序一致。下面是来源计算中蛋白组为 1 的情况，自己的组号请先核对：

printf '1\n' | gmx convert-tpr -s md.tpr -o protein.tpr
printf '1\n1\n' | gmx trjconv -s md.tpr -f md.xtc -o protein.xtc -pbc mol -center
printf '1\n' | gmx trjconv -s md.tpr -f md.xtc -o protein.gro -dump 0

资料包已备好 inputs/protein，手动分析可直接复制：

module purge
module load gromacs/2026.4-p100-threadmpi
mkdir helix-manual
cp inputs/protein/* helix-manual/
cd helix-manual
gmx select -s protein.tpr -select 'resnr 23 to 34 and name CA' -on helix-ca.ndx
gmx select -s protein.tpr -select 'protein' -on helix-protein.ndx

2. 计算局部与整体螺旋几何#

printf '0\n' | gmx helixorient -s protein.tpr -f protein.xtc -n helix-ca.ndx -oaxis helix-axis.dat -ocenter helix-center.dat -orise helix-rise.xvg -oradius helix-radius.xvg -otwist helix-twist.xvg -obending helix-bend.xvg -otilt helix-tilt.xvg -orot helix-rotation.xvg
printf '0\n' | gmx helix -s protein.tpr -f protein.xtc -n helix-protein.ndx -ahxstart 23 -ahxend 34 -cz helix-z.gro

helixorient 用连续 Cα 定义局部几何；helix 需要完整蛋白主链来构造相应属性，不能用只含 Cα 的组代替。本例人工指定 23–34 残基，换成其他蛋白前应先确认二级结构与编号。

3. 读取半径、升高与扭转#

head -30 radius.xvg
head -30 helix-radius.xvg

helix-radius.xvg 为每帧多个位置的局部量，本例是时间列加 12 列 Cα 数据。两端由于局部定义缺少邻居而出现零值，统计平均时排除这些未定义端点；不要把它们当作真实半径为零。

4. 两个容易读错的旧输出#

helicity.xvg 在本例是逐残基的 74 行统计，不是 51 帧时间序列。当前原生横坐标相对这份输入的实际残基编号偏移一位；配套 helicity-residue-corrected.csv 已按实际主链残基 2–75 映射。其他结构须重新核对编号，不能普遍执行“加一”。

cd222.xvg 是该工具旧实现中的经验二面角权重指标；原生单位标签不能当作经过标定的实验圆二色谱单位。可以辅助比较同一模型的构象变化，不应直接写成实验 CD 强度或完整 CD 光谱。

100 ps 是操作演示的短轨迹。研究螺旋稳定性需要更长采样、独立重复和统计误差。

另见 helix 官方说明、helixorient 官方说明。