GROMACS 碳链取向序参量操作

下载输入、真实结果与绘图资料包，解压进入资料包目录。inputs 可修改，reference 是集群已计算的参考结果；脚本每次在 work 下建立独立目录。计算在分配好的资源内执行，可用 sbatch run.sbatch 跑整组，结果路径写入 LAST-RUN.txt。

先用已知几何理解输出#

本例有 8 条、每条 12 个碳位点的锯齿骨架，按 0° 到 30° 构造 21 帧刚性倾转轨迹。这是人为构造的几何测试，不是 MD、真实脂膜或已经平衡的烷烃体系。 这样可以独立知道方向，核对 order 的索引和定义。

bash run.sh tools

1. 按等价碳位置分组#

inputs/order/carbons.ndx 有 C1–C12 共 12 组，每组包含 8 条链同一位置的碳。不要混入 System、Protein 等通用组，也不要把一整条链放成一组后期待得到逐碳结果。

cat inputs/order/carbons.ndx
cat inputs/order/prescribed-tilt.csv

本例分析结构的完整拓扑在 topol.top，analysis.mdp 只用于生成零步分析 TPR，没有启动 MD。实际客户使用脂膜时应换成自己真实拓扑和轨迹。

2. 运行完整输出与 Z 分量#

在复制了 inputs/order 的新目录执行：

gmx grompp -f analysis.mdp -c chain.gro -p topol.top -o chain.tpr
gmx order -s chain.tpr -f constructed.xtc -n carbons.ndx -d z -o order-tensor.xvg -od deuterium.xvg
gmx order -s chain.tpr -f constructed.xtc -n carbons.ndx -d z -szonly -o order-z.xvg -os order-slices.xvg -sl 1

对内部碳 i，工具通过 i−1 与 i+1 定义方向，因此两端不作为完整内部位点输出。本例原生文件的 1–10 对应实际 C2–C11。

3. 独立检查 Z 分量#

u = normalize(r(i+1) − r(i−1))
Szz = 0.5 × (3 × uz² − 1)

资料包从压缩后的 XTC 读坐标，逐帧、逐链、逐内部碳重新计算。原生结果平均约 0.867368，逐碳与独立坐标计算最大差约 1.4e-6。压缩坐标有有限精度，所以应按实际 XTC 核对，而不强求所有碳与未经压缩的理想解析值完全相同。

order-tensor.xvg 的分量和 deuterium.xvg 的 SCD 按工具的局部几何约定定义；不要将任意分量直接当成实验测得的脂膜氘序参量。实际膜还要确认法向、上下叶层、尾链映射、PBC、温压控制及充分平衡。

见 GROMACS order。