MDAnalysis扩散映射与结果检查

扩散映射用构象之间的相似程度构造低维坐标。本例输入是同一真实泛素轨迹的51帧、每帧76个Cα，统一对齐同一个参考。图中的低维点是构象几何的投影，不直接代表自由能，也不自动给出物理时间尺度。

1 下载并运行#

下载构象集合、扩散映射与并行密度完整示例，保存为ensemble-map.zip：

python3 -m zipfile -e ensemble-map.zip .
cd mdanalysis-ensemble-map
bash run.sh
bash verify.sh
bash plot.sh

使用集群现有MDAnalysis 2.9.0，不安装或修改共享环境。密度分析会启动2个工作进程；自己的大型轨迹应在已分配的CPU资源内运行。输出位于results/和plots/，数值保存在results/FULL-FIELDS.npz。再次运行请指定新目录，例如bash run.sh --output results-new。

2 先计算距离，再构造映射#

配套脚本实际运行：

from MDAnalysis.analysis.diffusionmap import DistanceMatrix, DiffusionMap
matrix = DistanceMatrix(universe, select='all', cutoff=0.0).run()
dmap = DiffusionMap(matrix, epsilon=1.0).run()
native_embedding = dmap.transform(3, time=1)

epsilon=1.0的单位是本例距离单位的平方，即Ų。构象距离为共同参考下的RMSD，没有为每一对构象再做一次最优叠合。高斯核为Kij=exp(-RMSDij**2/epsilon)，再用每行的核之和归一化，得到每行和为1的Markov矩阵P。time=1是对特征值取一次幂，不等于轨迹的1 ps。

3 2.9.0的原生特征向量需要特别核对#

本例的原生距离、完整51×51 Markov矩阵及全部特征值都能与独立计算对应，但原生特征向量未满足P @ V = V @ diag(lambda)。相对残差约0.0228982，大于本例固定的1e-10数值检查门限。因此包内保留原生输出作诊断，不把原生transform的坐标作为合格的科学结果。

这项说明针对本例所用的MDAnalysis 2.9.0及实际数据，不能据此假定所有版本都有同样行为。检查文件为results/DIFFUSION-NATIVE-VALIDATION.json；集群安装的库未作任何修改。

source env.sh
"$YS_PY" - <<'PYCODE'
import json
from pathlib import Path
r = json.loads(Path('results/DIFFUSION-NATIVE-VALIDATION.json').read_text())
print('native vectors qualified:', r['passed'])
print('relative eigen residual:', r['relative_eigen_residual'])
PYCODE

4 配套的独立求解如何使用#

为了仍能展示可靠的计算流程，analyze.py另行求解对称相似矩阵S=D^(-1/2) K D^(-1/2)，再转换为P的右特征向量。特征向量按列存放，特征值按降序排列；跳过第一条平凡模态后得到三维映射。这是配套脚本自身的求解结果，没有覆盖原生方法。

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
f = np.load('results/FULL-FIELDS.npz')
P = f['diffusion_Markov']
V = f['diffusion_own_eigenvectors_columns']
value = f['diffusion_own_eigenvalues']
residual = np.linalg.norm(P @ V - V * value[None, :]) / np.linalg.norm(V)
print('own solver residual:', residual)
print('row-sum error:', abs(P.sum(axis=1) - 1).max())
np.savetxt('diffusion-own-coordinates.dat', f['diffusion_own_embedding'])
PYCODE

本例独立求解的相对残差约1.01e-16。低维坐标的正负号本身可以变化，退化模态还允许子空间旋转，因此比较不同运行时应先比较方程残差、特征值及子空间，而不是仅凭散点图方向判断是否一致。

用于自己的数据时，先核对原子选择、参考对齐和距离单位，再测试不同epsilon、采样间隔与轨迹长度。本例没有把51帧的几何投影当成已收敛的构象自由能面。

参考：MDAnalysis 2.9.0 DiffusionMap文档与源码。