MDAnalysis蛋白内坐标与螺旋几何

BAT可以在笛卡尔坐标和非冗余内坐标间转换，HELANAL可以从连续Cα坐标描述螺旋局部几何。本文使用同一真实蛋白轨迹，展示完整坐标重建和螺旋参数提取。

1 下载并运行真实示例#

下载泛素蛋白八类分析输入、轨迹、结果与绘图包，保存为 mdanalysis-protein-wave326-v3.zip：

module load pymatgen/2026.9.24
python -m zipfile -e mdanalysis-protein-wave326-v3.zip .
cd mdanalysis-protein
bash run.sh
source env.sh

run.sh在新目录中完成八类分析并生成results/和plots/，已有同名输出时先换一个解压目录。环境使用集群已安装的Python 3.11 / MDAnalysis 2.9.0；env.sh选择相应Python并限制数值库为单线程，没有调用qvasp或VASP。其他机器可把ANALYSIS_PYTHON设为含相同依赖的Python路径用于run.sh。

本包使用鸢算科技集群已完成的GROMACS泛素轨迹：1231个原子、76个残基、1237条键。原始轨迹有1001帧，间隔0.1 ps；本例每20帧取一帧，得到0、2、…、100 ps共51帧。这里重用真实轨迹进行分析，没有重新做分子动力学；100 ps的变化不能说明蛋白折叠、长时稳定性或构象自由能已经收敛。

2 先处理拓扑和周期性边界#

本例GROMACS 2026.4生成的TPR格式为138，MDAnalysis 2.9.0不能直接读取。包中已提供由原生GROMACS导出的protein-conect.pdb及质量、电荷、元素元数据，不靠PDB猜测来代替拓扑。需要重做导出时：

bash export-topology.sh
source env.sh
python - <<'PYCODE'
from load import load
u = load()
print(len(u.atoms), len(u.residues), len(u.bonds))
print(len(u.trajectory), u.trajectory.dt)
print(len(u.atoms.fragments), u.atoms.charges.sum())
PYCODE

前两行分别得到1231 76 1237及1001 0.1；只有一个连通蛋白片段，总电荷约为数值零。export-topology.sh实际使用：

module load gromacs/2026.4-p100-threadmpi
# export-topology.sh已建立topology-export目录
gmx editconf -f inputs/protein.tpr -o topology-export/protein-conect.pdb -conect
gmx dump -s inputs/protein.tpr > topology-export/native-topology.dump
python topology_metadata.py

重新导出的文件在topology-export/，不会覆盖已有inputs/。topology_metadata.py从原生TPR转储提取质量、电荷和原子序号；本例脚本限定为这个单蛋白拓扑，不要直接当成任意多分子溶剂体系的通用转换器。load()随后加载PDB和XTC、补齐明确元数据，并用键连接执行unwrap让蛋白完整。自己的轨迹必须先确认原子顺序、残基命名、键和盒信息匹配；配体、断链或多链体系要相应修改选择。

3 BAT：3N坐标如何转为内坐标#

python - <<'PYCODE'
import numpy as np
from load import load
from MDAnalysis.analysis.bat import BAT
u = load()
a = BAT(u.atoms).run(stop=1001, step=20)
x = a.results.bat
print(x.shape)
y = a.Cartesian(x[0])
u.trajectory[0]
print(np.max(np.abs(y-u.atoms.positions)))
np.save('my-bat.npy', x)
PYCODE

1231个原子的BAT向量长度为 3693=3N，本例数组为(51,3693)。它除了键长、键角、扭转角，还包含描述整体平移和取向的自由度，不是3693个键长。转换需要正确的连通键图，因此不能删掉PDB的CONECT再依赖默认猜键。

全51帧重建的最大坐标误差为 7.10253×10⁻⁶ Å。这是同一套有限精度输入的数值往返误差，不能据此声称原分子动力学轨迹具有微埃物理精度。不同蛋白、环状结构、断链或多片段应先检查BAT适用条件和选区。

4 HELANAL：选一段真实的连续螺旋#

本例选择第23–34残基的12个Cα，直接调用原生helix_analysis并使用float64运算：

python - <<'PYCODE'
import numpy as np
from load import load
from MDAnalysis.analysis.helix_analysis import helix_analysis
u = load()
ca = u.select_atoms('name CA and resid 23:34')
twist, rise, turns = [], [], []
for ts in u.trajectory[0:1001:20]:
    r = helix_analysis(np.asarray(ca.positions, dtype=np.float64))
    twist.extend(r['local_twists'])
    rise.extend(r['local_heights'])
    turns.extend(r['local_nres_per_turn'])
print('twist / degree:', np.mean(twist))
print('rise / Angstrom:', np.mean(rise))
print('residues per turn:', np.mean(turns))
PYCODE

实际结果：局部扭转平均 99.65709°、局部升高 1.453786 Å、每圈残基数平均 3.630998。平均每圈残基数不是360/平均扭转角，因为先对每个局部几何取倒数再平均，与先平均再取倒数不同。

local_twists按四个连续Cα窗口定义，local_bends反映局部轴间弯折，global_axis描述整体轴方向。只有连续、合适的螺旋选区才有明确意义；不要把整个蛋白的Cα都交给螺旋算法，随后把β折叠和无规区域的输出也解释成螺旋。HELANAL官方手册列出每项几何量的定义。

5 运算精度和坐标精度分开看#

同一几何整体旋转90°后，局部扭转和弯折理应不变。本例MDAnalysis类接口的float32计算会出现约 0.000384927° 的最大差异；上面的原生函数采用float64运算后，差异约为 4×10⁻¹³°。因此配套脚本选择float64用于精细角度比较。

这只减少后处理舍入误差，不会增加XTC原有的坐标信息。常规结构展示也不需要把如此多小数位作为物理结论。precision-limits.json保留两种算法路径的对照；完整字段、逐帧结果和图仍可通过以下命令读取：

head -5 results/summary.csv
python - <<'PYCODE'
import numpy as np
z = np.load('results/FULL-FIELDS.npz')
print(z['BAT'].shape)
print(z['HELANAL64_local_twists'].shape)
print(z['HELANAL64_global_axis'].shape)
PYCODE

继续阅读二级结构与二面角和接触、GNM与PCA。换成自己的体系时先改加载器、选择和参考数据，不能保留本例固定的残基数和结构参照。