MDAnalysis距离矩阵、配对距离与接触图
用真实蛋白快照区分完整距离、自距离上三角、逐对距离和稀疏接触矩阵，并展示编号偏移。
阅读本篇目录
“算距离”可能是所有原子两两距离，也可能是一一配对的距离，或只判断接触。本篇把这些常见需求分清楚，并给出可直接导出的真实例子。
下载三个真实快照、完整拓扑、脚本及结果图包，保存为 example.zip：
python3 -m zipfile -e example.zip .
cd mdanalysis-distance-tools
sbatch job.slurm
使用现有MDAnalysis2.9.0及qvasp2.26 Python环境，仅做轨迹分析，不调用qvasp/VASP。CPU分区按56核整节点申请，实际1核；分区维护时等待恢复，在已分配CPU作业内可 bash run.sh。两篇共用一次运行；已有 VALIDATION.json时拒绝覆盖，重新计算请新目录。
输入是真实泛素轨迹0、50、100ps的三个完整分子快照，每帧1231原子、76个Cα；只是从已有轨迹提取坐标，没有新增模拟。原始子集在 inputs/protein-three-snapshots.npz，PDB_CONECT与原生GROMACS拓扑元数据同时保留。结果为 VALIDATION.json、完整 analysis-fields.npz 和3图，原件在 reference-results/。
1. 读取三个完整快照#
import numpy as np
import MDAnalysis as mda
from MDAnalysis.coordinates.memory import MemoryReader
from MDAnalysis.analysis import distances
x=np.load('inputs/protein-three-snapshots.npz')
u=mda.Universe('inputs/protein-conect.pdb')
u.load_new(x['positions_A'],format=MemoryReader,dt=50)
ca=u.select_atoms('name CA')
本例快照已经按真实连通关系使蛋白完整；坐标单位Å、时间ps。自己的XTC或其他轨迹应先核对拓扑顺序、分子完整性与盒信息，不能只凭PDB原子数相同就认为顺序正确。
2. 完整矩阵与自距离压缩数组#
full=distances.distance_array(ca.positions,ca.positions)
condensed=distances.self_distance_array(ca.positions)
76个Cα的完整矩阵是76×76，包含对角0和对称的重复项；自距离数组只存i<j，共76×75/2=2850项。它不是76个距离，也不能不经索引还原就reshape成76×76。
3. 一一配对的dist及编号偏移#
paired=distances.dist(ca[:-1],ca[1:],offset=(100,200))
两组各75个原子，第一个与第一个配对，以此类推。返回形状为3×75：第一行A的resid+100，第二行B的resid+200，第三行为距离。offset=10则给两边同时加10；偏移只改变返回的编号，不改变坐标和距离。本包实际比较了标量与双偏移两种调用。
4. 接触布尔矩阵与稀疏格式#
dense=distances.contact_matrix(ca.positions,cutoff=4.5,returntype='numpy')
sparse=distances.contact_matrix(ca.positions,cutoff=4.5,returntype='sparse')
两者都判断4.5Å内接触，输出全矩阵或稀疏表示。本例对三帧的全部元素核对一致；稀疏数组的CSR数据、列索引、行指针也保存在科学数组中。
原生接触矩阵包含对角自接触。 要统计无向、非自身的接触数，用 np.triu(dense,1).sum()，不要直接求整个矩阵和；后者会包含对角，还把每对接触数两次。
| 真实时间/ps | 无向非自接触数 | 几何交集原子数（下篇说明） |
|---|---|---|
| 0 | 82 | 4 |
| 50 | 81 | 2 |
| 100 | 83 | 2 |
5. 在终端导出矩阵和配对表#
export PATH="/public/toolkit/miniconda3/envs/qvasp-2.26/bin:$PATH"
python - <<'PYDATA'
import numpy as np
x=np.load('analysis-fields.npz')
np.savetxt('CA-distance-frame0.dat',x['frame0_distance_A'])
np.savetxt('CA-contact-frame0.dat',x['frame0_dense_contacts'],fmt='%d')
np.savetxt('CA-paired-frame0.csv',x['frame0_paired_dist'].T,delimiter=',',header='resid_A_plus100,resid_B_plus200,distance_A',comments='')
PYDATA
head -n 5 CA-paired-frame0.csv
距离、接触布尔量和编号均逐项与原始坐标独立核对，不把“脚本退出成功”当成数据正确。几何桥连与PBC边界见桥连选区篇。参考MDAnalysis2.9距离辅助工具。