MDAnalysis膜叶分类与距离阈值

膜叶分类可以从头基之间的邻接图入手：附近头基连接成图，再找连通分组。它依赖选择、距离阈值和周期性边界，不是简单按z坐标截一刀；本例用原生工具与完整几何图进行对照。

1 下载并运行#

下载RNA与膜叶分类的完整输入、分析脚本、数值和结构图，保存为rna-membrane.zip：

module load pymatgen/2026.9.24
python -m zipfile -e rna-membrane.zip .
cd mdanalysis-rna-membrane
bash run.sh
source env.sh

输出在results/、plots/，请使用新的解压目录。环境是MDAnalysis 2.9.0和NetworkX 3.6.1；setup.sh在示例目录建立NetworkX的局部加载入口，使用集群已有依赖，不修改公共Python环境。env.sh供当前终端中的Python示例使用；没有调用qvasp或VASP。

数据来自MDAnalysisTests 2.9.0官方公开样例，来源、哈希及上游许可随包保留。RNA是PDB 1K5I派生数据中RNAA段的23个核苷酸、739个原子；膜是单个DPPC/胆固醇结构快照。这里在集群上实际执行几何和图分析，没有把这些上游快照说成新做的分子动力学。程序中的90°刚体旋转只检查几何量是否保持不变，不是第二个动力学时间点。

2 选什么原子#

本结构有5040个粗粒化珠子、450个脂质残基，其中360个PO4磷脂头基。本例只用name PO4构图；胆固醇没有被误当成另一个PO4头基。这是Martini粗粒化结构，不是5040个真实化学原子。

python - <<'PYCODE'
import MDAnalysis as mda
u = mda.Universe('inputs/martini_dppc_chol_bilayer.gro')
heads = u.select_atoms('name PO4')
print(len(u.atoms), len(u.residues), len(heads))
print(u.dimensions)
PYCODE

GRO坐标原生单位为nm，MDAnalysis读取后转换为Å；这里的15 Å阈值相当于1.5 nm，不要按15 nm输入。

3 LeafletFinder：周期图中的两组#

python - <<'PYCODE'
import MDAnalysis as mda
from MDAnalysis.analysis.leaflet import LeafletFinder
u = mda.Universe('inputs/martini_dppc_chol_bilayer.gro')
a = LeafletFinder(u, 'name PO4', cutoff=15.0, pbc=True, sparse=False)
for i,g in enumerate(a.groups()):
    print(i, len(g), g.center_of_geometry())
PYCODE

实际得到两组，分别 180、180 个头基。组号只表示连通分组，不能假设0永远是上膜叶。图中按两组头基的平均z值标注上、下组；换成弯曲膜或囊泡时，应按实际几何选择内外膜叶定义。

pbc=True考虑盒边界两侧的邻居，否则一张横跨盒子的膜可能被拆成几个小块。只看“分成两组”也不够，应确认所选确为合适头基，并检查分组大小及空间位置。

4 optimize_cutoff：检查阈值区间#

python - <<'PYCODE'
import MDAnalysis as mda
from MDAnalysis.analysis.leaflet import optimize_cutoff
u = mda.Universe('inputs/martini_dppc_chol_bilayer.gro')
r = optimize_cutoff(u, 'name PO4', dmin=8, dmax=21, step=1,
                    max_imbalance=0.2, pbc=True, sparse=False)
print(r)
PYCODE

在8–20 Å、每步1 Å的本次扫描中，原生推荐值为 11 Å，对应两个连通分组；示例15 Å也得到180/180。这个推荐依据连通分组数量、大小均衡等规则，不是对所有膜或力场通用的物理常数。阈值太小时可能把同一膜叶拆散，太大时可能把不同膜叶连在一起。

python - <<'PYCODE'
import numpy as np
z = np.load('results/FULL-FIELDS.npz')
print(z['membrane_profile'])
print(np.bincount(z['membrane_native_labels15']))
PYCODE

完整13个阈值的分组数与两大组大小都在包中；实际头基图和独立最小镜像几何连通结果一致。本结构只有一帧，因此没有据此计算膜叶交换、翻转速率或时间稳定性。相关接口见Leaflet官方手册。