MDAnalysis构象集合相似性

本例比较已有泛素轨迹的前26帧与后26帧。两段共享中间的一帧，全部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 HES比较什么#

谐和集合相似性HES把选中的笛卡尔坐标描述为多元高斯分布，比较均值和协方差。本例选前三个Cα，共9维坐标；它是一个清楚可计算的小型描述符，不是整个蛋白的谐振动或自由能计算。协方差采用最大似然定义，分母是样本数N，而不是N−1。

配套脚本使用原生数值接口，显式采用双精度：

from MDAnalysis.analysis.encore.covariance import ml_covariance_estimator
from MDAnalysis.analysis.encore.similarity import harmonic_ensemble_similarity
x = paths[0][:, :3].astype(float).reshape(26, 9)
y = paths[1][:, :3].astype(float).reshape(26, 9)
mean_x, mean_y = x.mean(axis=0), y.mean(axis=0)
cov_x = ml_covariance_estimator(x)
cov_y = ml_covariance_estimator(y)
hes_value = harmonic_ensemble_similarity(cov_x, cov_y, mean_x, mean_y)

本次HES为约4.401540，无量纲。完整均值、9×9协方差及独立公式的结果都在下载包中。独立公式为两组逆协方差的均值差二次型，加上交叉协方差的迹项，再乘1/4；它与原生数值函数一致。

原生便利接口hes(..., cov_estimator='ml', weights=None, align=False)也实际运行了，得到约4.401528。该版本便利接口读取的默认单精度轨迹均值存在约1.47e-6 Å的舍入差别；因此配套精确数值采用上述双精度数值接口，同时保留便利接口结果供比较。双精度不会提高原始轨迹本身的坐标精度。

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
f = np.load('results/FULL-FIELDS.npz')
print('double HES:', f['HES_matrix'][0, 1])
print('convenience HES:', f['HES_convenience_float32_matrix'][0, 1])
print('first ensemble mean:', f['HES_mean_0'])
np.savetxt('ensemble-first-covariance.dat', f['HES_cov_0'])
PYCODE

3 CES比较哪些构象出现过#

聚类集合相似性CES先把两段轨迹放到同一组几何聚类中，再比较各组占比。这里使用全部76个Cα的共同参考RMSD矩阵、原生AffinityPropagationNative：

method = AffinityPropagationNative(
    damping=0.9, preference=-0.5, max_iter=500,
    convergence_iter=50, add_noise=False
)
values, details = ces(
    ensembles, select='all', clustering_method=method,
    distance_matrix=triangular_rmsd_matrix, ncores=1
)

本次产生34个几何聚类，CES为0.666488。它采用自然对数的Jensen–Shannon散度，范围为0至ln(2)≈0.693147；不要把它和开平方后的Jensen–Shannon距离混为同一数值。本例关闭聚类加噪以便重算，聚类偏好参数是−0.5，换参数或增删帧可能改变分组。

source env.sh
"$YS_PY" - <<'PYCODE'
import numpy as np
f = np.load('results/FULL-FIELDS.npz')
p = f['CES_cluster_probabilities']
print('clusters:', len(p))
print('each ensemble probability sum:', p.sum(axis=0))
print('CES:', f['CES_matrix'][0, 1])
np.savetxt('cluster-occupancy.dat', p, header='first26 last26')
PYCODE

HES与CES采用的描述符和假设不同，数值尺度不能横向比较。34个几何聚类不等于34个亚稳态；本例也没有用较大的CES证明两个时间片属于不同热力学相。自己的研究应增加独立轨迹、比较分段结果、调整描述符和聚类参数，并单独检查采样是否充分。

参考：MDAnalysis 2.9.0 ENCORE相似性接口。该分支在2.9仍可使用，后续新版本的独立MDAKit接口需按所用版本查看。