MDAnalysis残基近接触网络
在真实泛素全原子轨迹上运行两种closeContactGNM权重，核对完整残基矩阵、本征值和模式方向。
阅读本篇目录
近接触GNM先统计原子接触，再把接触汇总到残基网络。它与只用Cα距离建立的GNM不同：原子选择、接触阈值及残基大小权重都会改变网络。
下载完整真实轨迹、拓扑、分析脚本和数据图包，先按平均结构篇的共用步骤运行一次，本篇读取同一批结果，不重复提交。
1. 选择完整蛋白并检查连通信息#
head -n 6 inputs/topology-metadata.json
cat VALIDATION.json
沿用完整分子预处理后的1231原子、76残基、51帧真实数据。当前实现的残基索引与size权重按残基全部原子计数，本例明确用 select='protein' 包含氢；不要把示例直接改为局部原子子集，再假定索引与权重仍有相同含义。
网络方法本身不使用PBC最小镜像。因此应先让蛋白成为完整分子；如果研究跨周期胞的聚集体或膜体系，不能自动沿用本例的几何网络定义。
2. 两种实际原生调用#
from MDAnalysis.analysis.gnm import closeContactGNMAnalysis
plain=closeContactGNMAnalysis(u,select='protein',cutoff=4.5,weights=None).run()
weighted=closeContactGNMAnalysis(u,select='protein',cutoff=4.5,weights='size').run()
value=weighted.results.eigenvalues
vector=weighted.results.eigenvectors
每帧都建立76×76矩阵。对不同残基i、j，数出4.5Å内的原子对：不加权用接触数n_ij；size权重用 n_ij/sqrt(Ni*Nj)，Ni、Nj是相应残基原子数。非对角元为负接触权重，对角元为该行接触权重之和，行和为0。
两种权重是两种不同网络，不是同一结果的单位换算。 当前输出是几何模型单位，不能直接把本征值开平方写成cm⁻¹，也没有拟合真实弹簧常数。
3. 原生矩阵与本征模式怎么读#
资料包用一个只记录矩阵的子类捕获原生 generate_kirchoff() 返回值，随后仍由原生算法求模式，没有替换库的物理计算。独立检查从全部原子坐标用另一种空间搜索计算接触，核对51张完整矩阵、对称性和行和。
每帧网络有一条均匀平移零模。原生结果是第一个非零本征值及方向，不是输出全部76条原生模式。完整独立谱在数组中单独标记为independent。
| 原生方向检查 | 51帧最大方程残差 |
|---|---|
| 不加权 | 2.5763e−13 |
| size权重 | 2.8222e−14 |
同时检查 K*v=lambda*v、方向单位长度和与独立模式的投影矩阵。方向整体正负号可以翻转；图中为了可读性按首帧符号对齐，不能把这种符号切换当作蛋白运动。
4. 导出第一帧矩阵与全部本征值#
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('contact-size-frame0.dat',x['contacts-size_native_K'][0])
np.savetxt('contact-size-modes.dat',np.column_stack([x['original_times_ps'],x['contacts-size_native_lambda']]),header='time_ps first_nonzero_native_eigenvalue')
PYDATA
读取自己的数据时，先确认蛋白是否完整、原子选择、残基编号、阈值及权重，再讨论网络变化；本例100ps有限轨迹不代表已收敛的蛋白柔性或真实振动频率。参考MDAnalysis2.9 GNM模型说明。