分子动力学模拟结果怎么分析：RMSD、RMSF、回转半径、氢键与 SASA

# 单链可溶蛋白（如溶菌酶）：两次调用即可
# 第 1 次：以 Protein 居中，整个体系按分子放回盒子；输出 System
printf "Protein\nSystem\n" | gmx trjconv -s md.tpr -f md.xtc -o md_center.xtc -pbc mol -center -ur compact
# 第 2 次（只用于可视化和 MDAnalysis）：以 Backbone 叠合，输出 System；此后不再做任何 PBC 处理
printf "Backbone\nSystem\n" | gmx trjconv -s md.tpr -f md_center.xtc -o md_fit.xtc -fit rot+trans

# 多链蛋白或蛋白-配体复合物：按官方顺序 whole -> nojump -> center -> fit
# 0) 建一个“蛋白+配体”组（配体残基名以 LIG 为例）
printf '"Protein" | "LIG"\nq\n' | gmx make_ndx -f md.tpr -o index.ndx
# 1) 修补被盒子切断的分子
printf "System\n" | gmx trjconv -s md.tpr -f md.xtc -o md_whole.xtc -pbc whole
# 2) 做一个蛋白与配体在同一侧的参考帧（导出后在 VMD/PyMOL 中确认），再以它为参考做 nojump
printf "Protein_LIG\nSystem\n" | gmx trjconv -s md.tpr -f md_whole.xtc -n index.ndx -o frame0.gro -dump 0 -pbc mol -center -ur compact
printf "System\n" | gmx trjconv -s frame0.gro -f md_whole.xtc -o md_nojump.xtc -pbc nojump
# 3) 以蛋白+配体居中；居中之后不能再用 -pbc nojump
printf "Protein_LIG\nSystem\n" | gmx trjconv -s md.tpr -f md_nojump.xtc -n index.ndx -o md_center.xtc -center
# 4) 只用于可视化和 MDAnalysis：以蛋白骨架叠合；之后不再使用任何 PBC 选项
printf "Backbone\nSystem\n" | gmx trjconv -s md.tpr -f md_center.xtc -n index.ndx -o md_fit.xtc -fit rot+trans

# 配体是真的解离还是显示问题：最小距离按周期性边界计算，不受 trjconv 处理影响
printf "Protein\nLIG\n" | gmx mindist -s md.tpr -f md.xtc -n index.ndx -od mindist_pl.xvg -tu ns
# 蛋白与自身周期性镜像的最小距离，应始终大于非键截断（rcoulomb、rvdw）
printf "Protein\n" | gmx mindist -s md.tpr -f md.xtc -pi -od mindist_pi.xvg -tu ns

# 蛋白骨架 RMSD：第一个提示是叠合组，第二个是计算组
printf "Backbone\nBackbone\n" | gmx rms -s md.tpr -f md_center.xtc -o rmsd_bb.xvg -tu ns
# 相对晶体结构（能量最小化前后的 em.tpr）的 RMSD
printf "Backbone\nBackbone\n" | gmx rms -s em.tpr -f md_center.xtc -o rmsd_xtal.xvg -tu ns
# 配体 RMSD：先按蛋白骨架叠合，再算配体，不对配体自身叠合
printf "Backbone\nLIG\n" | gmx rms -s md.tpr -f md_center.xtc -n index.ndx -o rmsd_lig.xvg -tu ns
# 只看配体自身构象变化（与上一条含义不同，论文中要分开写）
printf "LIG\nLIG\n" | gmx rms -s md.tpr -f md_center.xtc -n index.ndx -o rmsd_lig_self.xvg -tu ns
# all-to-all RMSD 矩阵，用于判断是否回到已采样构象
printf "Backbone\nBackbone\n" | gmx rms -s md.tpr -f md_center.xtc -m rmsd_matrix.xpm -dt 100

# 去掉前 10 ns（-b 单位为 ps，此命令不加 -tu），按残基输出 Cα RMSF
printf "C-alpha\n" | gmx rmsf -s md.tpr -f md_center.xtc -b 10000 -res -o rmsf.xvg -oq bfac.pdb

# 回转半径（GROMACS 2024 起的新实现，用 -sel 选择）
gmx gyrate -s md.tpr -f md_center.xtc -sel Protein -o gyrate.xvg -tu ns
# 溶剂可及表面积：-surface 包含全部非溶剂原子，-output 取其子集
gmx sasa -s md.tpr -f md_center.xtc -surface Protein -o sasa.xvg -or resarea.xvg -tu ns
# 蛋白-配体体系
gmx sasa -s md.tpr -f md_center.xtc -surface 'group "Protein" or resname LIG' \
         -output 'group "Protein"' 'resname LIG' -o sasa_pl.xvg -tu ns

# GROMACS 2024+ 新实现：-r 与 -t 都必须给；-num 必须显式写出才会输出
gmx hbond -s md.tpr -f md_center.xtc -r Protein -t Protein -num hbnum.xvg -tu ns
# 蛋白-配体氢键，-o 写出每对氢键的原子索引
gmx hbond -s md.tpr -f md_center.xtc -r Protein -t 'resname LIG' -num hbnum_pl.xvg -o hbond_pl.ndx -tu ns
# 需要氢键寿命或存在矩阵时用旧实现：-ac 给出 Luzar-Chandler 速率常数，-hbm 给出每对氢键逐帧存在矩阵
printf "Protein\nLIG\n" | gmx hbond-legacy -s md.tpr -f md_center.xtc -n index.ndx \
         -num hbnum_legacy.xvg -ac hbac.xvg -hbn hbond.ndx -hbm hbmap.xpm

# GROMACS 2023+ 内置 DSSP v4，不再需要外部 dssp 程序和 DSSP 环境变量
gmx dssp -s md.tpr -f md_center.xtc -sel Protein -o dssp.dat -num dssp_num.xvg -tu ns
# 结构不含氢（如粗处理后的 PDB）：由 C、O 生成氢的伪原子
gmx dssp -s protein_noH.pdb -f protein_noH.pdb -sel Protein -hmode dssp -clear -o dssp_noH.dat

# 块平均误差估计（Hess 2002 方法），只用生产段
gmx analyze -f gyrate.xvg -b 10 -ee gyrate_errest.xvg
# 余弦含量：每次只投影一个主成分，再交给 gmx analyze -cc
printf "C-alpha\nC-alpha\n" | gmx covar -s md.tpr -f md_center.xtc -b 10000 -o eigenval.xvg -v eigenvec.trr -av average.pdb
printf "C-alpha\nC-alpha\n" | gmx anaeig -s md.tpr -f md_center.xtc -b 10000 -v eigenvec.trr -eig eigenval.xvg -first 1 -last 1 -proj pc1.xvg
gmx analyze -f pc1.xvg -cc pc1_cosine.xvg
# 前后两半各做一次 covar（此例生产段 10-60 ns），再比较前 10 个本征向量张成的子空间
printf "C-alpha\nC-alpha\n" | gmx covar -s md.tpr -f md_center.xtc -b 10000 -e 35000 -o eigval_h1.xvg -v eigvec_h1.trr
printf "C-alpha\nC-alpha\n" | gmx covar -s md.tpr -f md_center.xtc -b 35000 -e 60000 -o eigval_h2.xvg -v eigvec_h2.trr
gmx anaeig -v eigvec_h1.trr -v2 eigvec_h2.trr -first 1 -last 10 -over overlap.xvg

# 延长已结束的模拟（单位 ps，此处加 100 ns）
gmx convert-tpr -s md.tpr -extend 100000 -o md_ext.tpr
gmx mdrun -deffnm md -s md_ext.tpr -cpi md.cpt   # 默认追加写入原文件
# 若用了 -noappend，会得到 md.part0002.xtc 等文件（编号是模拟段号，之后每段递增），分析前先拼接
gmx trjcat -f md.xtc md.part0002.xtc -o md_all.xtc

# 1) 协方差矩阵：只用 Cα，去掉平衡段（-b 单位 ps）
printf "C-alpha\nC-alpha\n" | gmx covar -s md.tpr -f md_center.xtc -b 10000 -o eigenval.xvg -v eigenvec.trr -av average.pdb
# 2) 投影到 PC1、PC2；-2d 输出两列（PC1 PC2），不含时间列
printf "C-alpha\nC-alpha\n" | gmx anaeig -s md.tpr -f md_center.xtc -b 10000 -v eigenvec.trr -eig eigenval.xvg \
       -first 1 -last 2 -2d 2dproj.xvg
# 3) 沿 PC1 的两端极值结构，用于看运动方向
printf "C-alpha\nC-alpha\n" | gmx anaeig -s md.tpr -f md_center.xtc -v eigenvec.trr -first 1 -last 1 -extr pc1_extreme.pdb -nframes 10
# 4) 自由能形貌：必须 -notime；-tsham 改成模拟温度
gmx sham -f 2dproj.xvg -notime -tsham 300 -ngrid 40 40 40 -nlevels 50 -ls fel.xpm -lp prob.xpm

# 水的自扩散系数（默认在 MSD 曲线 10%-90% 区间拟合）
gmx msd -s md.tpr -f md_nojump.xtc -sel 'resname SOL and name OW' -o msd.xvg
# 配体周围水氧的径向分布函数
gmx rdf -s md.tpr -f md_center.xtc -ref 'resname LIG' -sel 'resname SOL and name OW' -selrpos res_com -o rdf.xvg

# 用法：python plot_replicas.py "rep*/rmsd.xvg" 10 RMSD_nm
# 参数：xvg 通配符、丢弃的平衡段（与 xvg 时间列同单位）、y 轴标签
import glob
import sys

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np


def read_xvg(path):
    """读取 GROMACS .xvg：跳过 # 和 @ 注释行，遇到 & 停止（只读第一个数据集）"""
    rows = []
    with open(path) as fh:
        for line in fh:
            s = line.strip()
            if not s or s[0] in "#@":
                continue
            if s.startswith("&"):
                break
            rows.append([float(v) for v in s.split()])
    return np.array(rows)


def block_se(x, n_blocks):
    """块平均标准误：把序列等分成 n_blocks 块，返回块均值的标准误"""
    m = len(x) // n_blocks
    blocks = x[: m * n_blocks].reshape(n_blocks, m).mean(axis=1)
    return blocks.std(ddof=1) / np.sqrt(n_blocks)


pattern, t_equil, ylabel = sys.argv[1], float(sys.argv[2]), sys.argv[3]
files = sorted(glob.glob(pattern))
if len(files) < 2:
    sys.exit("need at least 2 replicas, found %d" % len(files))
data = [read_xvg(f) for f in files]
n = min(len(d) for d in data)
t = data[0][:n, 0]
for f, d in zip(files, data):
    if not np.allclose(d[:n, 0], t):
        sys.exit("time axis differs: %s" % f)
Y = np.stack([d[:n, 1] for d in data])  # 第 2 列；Rg 的总值也在第 2 列
mean, sd = Y.mean(axis=0), Y.std(axis=0, ddof=1)

fig, ax = plt.subplots(figsize=(6, 3.5), dpi=200)
for f, y in zip(files, Y):
    ax.plot(t, y, lw=0.6, alpha=0.5, label=f.split("/")[0])
ax.plot(t, mean, lw=1.5, color="black", label="mean")
ax.fill_between(t, mean - sd, mean + sd, color="gray", alpha=0.3, label="±1 SD")
ax.axvline(t_equil, ls="--", lw=0.8, color="gray")
ax.set_xlabel("Time (same unit as xvg)")
ax.set_ylabel(ylabel)
ax.legend(fontsize=6, frameon=False)
fig.tight_layout()
fig.savefig("replicas.png")

# 生产段统计：每个副本先求均值，再用副本间标准差 / sqrt(n) 得标准误
prod = Y[:, t >= t_equil]
rep_means = prod.mean(axis=1)
se = rep_means.std(ddof=1) / np.sqrt(len(rep_means))
print("replica means:", np.round(rep_means, 4))
print("mean = %.4f, SE across replicas = %.4f (n=%d)" % (rep_means.mean(), se, len(rep_means)))
# 块平均：标准误随块变大应趋于平台；不出现平台说明采样不足
for nb in (40, 20, 10, 5):
    if prod.shape[1] < 2 * nb:  # 每块至少 2 帧，帧数太少时跳过
        continue
    print("rep1 blocks=%2d  SE=%.4f" % (nb, block_se(prod[0], nb)))

# 用法：python fel_from_2dproj.py 2dproj.xvg 300
# 与 gmx sham 相同的定义：G = -kT ln P，再平移使最低点为 0
import sys

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

path, T = sys.argv[1], float(sys.argv[2])
xy = np.array([[float(v) for v in l.split()[:2]] for l in open(path)
               if l.strip() and l.strip()[0] not in "#@&"])
H, xe, ye = np.histogram2d(xy[:, 0], xy[:, 1], bins=40)
P = H / H.sum()
kT = 0.0083144626 * T  # kJ/mol
with np.errstate(divide="ignore"):
    G = -kT * np.log(P)
G -= G[np.isfinite(G)].min()
fig, ax = plt.subplots(figsize=(4.5, 3.8), dpi=200)
im = ax.pcolormesh(xe, ye, np.ma.masked_invalid(G).T, cmap="viridis", shading="flat")
fig.colorbar(im, label="G (kJ/mol)")
ax.set_xlabel("PC1 (nm)")
ax.set_ylabel("PC2 (nm)")
fig.tight_layout()
fig.savefig("fel.png")
print("frames: %d, empty bins: %d / %d" % (len(xy), int((H == 0).sum()), H.size))
# 最低自由能格子中的帧序号（从 anaeig -b 指定的起点开始计数），用 gmx trjconv -dump 导出代表构象
i, j = np.unravel_index(np.nanargmin(np.where(np.isfinite(G), G, np.nan)), G.shape)
ix = np.clip(np.digitize(xy[:, 0], xe) - 1, 0, len(xe) - 2)
iy = np.clip(np.digitize(xy[:, 1], ye) - 1, 0, len(ye) - 2)
print("frames in the minimum bin:", np.where((ix == i) & (iy == j))[0][:20])

# 用法：python mda_analysis.py frame0.pdb md_fit.xtc
# MDAnalysis 2.10 读不了 GROMACS 2026 的 tpr（报 Your tpx version is 138, which this parser does not support），
# 先用 gmx trjconv -dump 0 导出 frame0.pdb 作拓扑；能读 tpr 的组合也可以直接传 md.tpr
import sys

import MDAnalysis as mda
import mdtraj as md
import numpy as np
from MDAnalysis.analysis import align, rms
from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis

top, traj = sys.argv[1], sys.argv[2]
u = mda.Universe(top, traj)
protein = u.select_atoms("protein")
has_lig = len(u.select_atoms("resname LIG")) > 0

# 1) RMSD：按骨架叠合；配体 RMSD 用 groupselections，在骨架叠合后计算且不再单独叠合
R = rms.RMSD(u, u, select="backbone",
             groupselections=["resname LIG"] if has_lig else None, ref_frame=0).run()
np.savetxt("rmsd_mda.dat", R.results.rmsd, fmt="%.4f",
           header="frame time_ps backbone_A" + (" ligand_A" if has_lig else ""))

# 2) 回转半径（Å，质量加权）
rg = np.array([(ts.time, protein.radius_of_gyration()) for ts in u.trajectory])
np.savetxt("rg_mda.dat", rg, fmt="%.4f", header="time_ps rg_A")

# 3) 氢键：判据写成与 GROMACS 默认接近的 3.5 Å；论文中写明所用判据
sel = "protein or resname LIG" if has_lig else "protein"
if hasattr(u.atoms, "charges"):  # tpr 拓扑：按电荷猜测氢和受体
    hb = HydrogenBondAnalysis(u, between=["protein", "resname LIG"] if has_lig else None,
                              d_a_cutoff=3.5, d_h_a_angle_cutoff=150)
    hb.hydrogens_sel = hb.guess_hydrogens(sel)
    hb.acceptors_sel = hb.guess_acceptors(sel)
else:  # PDB 拓扑没有电荷：按原子名显式指定 N/O 供体、受体和氢
    hb = HydrogenBondAnalysis(u, between=["protein", "resname LIG"] if has_lig else None,
                              donors_sel="(%s) and (name N* or name O*)" % sel,
                              hydrogens_sel="(%s) and name H*" % sel,
                              acceptors_sel="(%s) and (name O* or (resname HI* and name ND1 NE2))" % sel,
                              d_a_cutoff=3.5, d_h_a_angle_cutoff=150)
hb.run()
np.savetxt("hbnum_mda.dat", np.column_stack([hb.times, hb.count_by_time()]),
           fmt="%.3f", header="time_ps n_hbonds")
id2ix = {i: k for k, i in enumerate(u.atoms.ids)}  # count_by_ids 返回原子 id（从 1 开始），不是索引
for d, h, a, n in hb.count_by_ids()[:15]:  # 占有率最高的 15 对
    D, A = u.atoms[id2ix[d]], u.atoms[id2ix[a]]
    print("%s%d:%s -> %s%d:%s  occupancy=%.2f" % (
        D.resname, D.resid, D.name, A.resname, A.resid, A.name, n / u.trajectory.n_frames))
tau, acf = hb.lifetime(tau_max=50)  # 连续型自相关；tau 以帧为单位
np.savetxt("hb_lifetime_acf.dat", np.column_stack([tau, acf]), fmt="%.4f")

# 4) RMSF：MDAnalysis 的 RMSF 不做叠合，先对 Cα 平均结构对齐（in_memory 会改写 u 中的坐标）
avg = align.AverageStructure(u, u, select="protein and name CA", ref_frame=0).run()
align.AlignTraj(u, avg.results.universe, select="protein and name CA", in_memory=True).run()
ca = u.select_atoms("protein and name CA")
F = rms.RMSF(ca).run()
np.savetxt("rmsf_mda.dat", np.column_stack([ca.resids, F.results.rmsf]), fmt="%.4f",
           header="resid rmsf_A")

# 5) SASA：用 MDTraj 的 Shrake-Rupley（探针 0.14 nm，结果单位 nm^2）
t = md.load(traj, top="frame0.pdb")
t = t.atom_slice(t.topology.select("protein"))
sasa = md.shrake_rupley(t, probe_radius=0.14, mode="residue")
np.savetxt("sasa_total_nm2.dat", np.column_stack([t.time, sasa.sum(axis=1)]), fmt="%.3f")

# DuIvyTools 0.6.0（pip 安装，命令名 dit）
pip install DuIvyTools
dit xvg_show -f rmsd.xvg                      # 查看单条曲线
dit xvg_compare -f rep1/rmsd.xvg rep2/rmsd.xvg rep3/rmsd.xvg -c 1 1 1   # 多副本同图比较
dit xpm_show -f fel.xpm                       # 查看 gmx sham 输出的自由能形貌图
dit xpm2csv -f fel.xpm -o fel.csv             # 转成 x, y, z 三列，交给 Origin 等软件作图
dit dssp -f dssp.dat -o dssp.xpm -x "Frame"     # 把 gmx dssp 的 dssp.dat 转成旧式 xpm 和两个 xvg