AMBER 轨迹分析与作图（上）：从轨迹文件到第一张图

本文收录于专栏分子动力学模拟-Amber—— 专栏系统覆盖Amber 分子动力学模拟全流程，点击订阅可跟踪后续更新。

摘要：本文系统梳理 Amber 分子动力学模拟轨迹的常用分析工具与完整流程，涵盖 cpptraj/ptraj/VMD 的选型对比，以及 RMSD、RMSF、回旋半径、氢键、距离/二面角、PCA 六类常见分析的 cpptraj 输入文件模板。文章还给出从数据导出到论文级作图的四种方式（Matplotlib、xmgrace、VMD、Jupyter+NGLview+MDAnalysis），并提供一个从对齐骨架到 RMSD/RMSF/Rg/关键残基距离出图的完整案例教程，帮助读者跑完模拟后直接套用命令与脚本生成分析图。

你跑完一条 MD 轨迹，面对md.nc几万个帧不知道从哪下手：RMSD、RMSF、Rg、氢键、PCA 各代表什么、cpptraj 命令怎么写、算完的 .dat 又怎么画成论文级图？这篇上篇把常用分析工具选型（cpptraj/ptraj/VMD）、六类常见分析的 cpptraj 输入文件模板（rmsd.in/rmsf.in/rg.in/hbond.in/distance.in/pca.in）、以及一个从对齐骨架到 RMSD/RMSF/Rg/关键残基距离出图的完整案例教程一次给你。你跑完模拟后照着抄命令、换 mask 就能出自己的分析图；参数级详解与 DSSP/聚类等进阶内容见下篇。

相关教程与核心文献

一、Amber 轨迹分析常用工具

1.cpptraj（推荐）

• 自 AmberTools 13（Amber14 时代）起 cpptraj 就是默认轨迹分析工具（C++ 重写版，比 ptraj 更快、功能更强）
• 支持脚本化、批处理、多线程
• 可直接输出数据供绘图

2.ptraj（旧版，已逐渐淘汰）

• 功能类似 cpptraj，但速度慢、语法略不同

3.VMD + MMTK / NGLview（Jupyter）

• 用于可视化轨迹（非定量分析）
• 可配合 Python 脚本提取数据

二、常见分析类型与 cpptraj 命令示例

假设你已有：

• 拓扑文件：system.prmtop
• 轨迹文件：md.nc（NetCDF 格式）或md.mdcrd

1. RMSD（均方根偏差）

# rmsd.in parm system.prmtop trajin md.nc rms first :1-100&!@H= # 对残基1-100的重原子计算相对于第一帧的RMSD rmsd out rmsd.dat first run

运行：

cpptraj -i rmsd.in

输出rmsd.dat，格式为：帧号 RMSD(Å)

2. RMSF（均方根涨落）

# rmsf.in parm system.prmtop trajin md.nc rms first :1-100&!@H= # 先对结构做对齐 atomicfluct out rmsf.dat :1-100&!@H= byres run

输出每残基的 RMSF（Å）

3. 回旋半径（Radius of Gyration, Rg）

# rg.in parm system.prmtop trajin md.nc radgyr out rg.dat :1-100 run

4. 氢键分析

# hbond.in parm system.prmtop trajin md.nc hbond out hbonds.dat series run

可输出氢键数量随时间变化，或具体氢键列表

5. 距离/二面角/角度

# distance.in parm system.prmtop trajin md.nc distance d1 :10@CA :20@CA out dist_CA.dat run

6. PCA（主成分分析）

# pca.in parm system.prmtop trajin md.nc rms first :1-100&!@H= matrix covar name mycovar :1-100&!@H= run runanalysis diagmatrix mycovar out evecs.dat vecs 3 name pc projection modes pc beg 1 end 3 :1-100&!@H= out pca.dat run

后续可用projection投影到主成分空间

三、作图方式

cpptraj 本身不绘图，需导出数据后用其他工具绘图。

方式1：Python + Matplotlib（推荐）

示例：绘制 RMSD

import matplotlib.pyplot as plt import numpy as np data = np.loadtxt('rmsd.dat', comments=['#', '@']) time = data[:, 0] * 0.002 # 假设每帧间隔2 ps，单位 ns rmsd = data[:, 1] plt.figure(figsize=(8, 5)) plt.plot(time, rmsd, color='blue') plt.xlabel('Time (ns)') plt.ylabel('RMSD (Å)') plt.title('Backbone RMSD over Time') plt.grid(True) plt.tight_layout() plt.savefig('rmsd.png', dpi=300) plt.show()

示例：RMSF（按残基）

data = np.loadtxt('rmsf.dat', comments=['#', '@']) residues = data[:, 0] rmsf = data[:, 1] plt.bar(residues, rmsf, width=1.0, color='lightcoral') plt.xlabel('Residue Number') plt.ylabel('RMSF (Å)') plt.title('Per-residue RMSF') plt.savefig('rmsf.png', dpi=300)

方式2：xmgrace（Grace）

• 适用于快速查看.dat文件
• 在 Linux 下安装：sudo apt install grace
• 命令：xmgrace rmsd.dat

可手动调整样式、导出 EPS/PNG

方式3：VMD（可视化+简单分析）

• 打开 VMD → Loadsystem.prmtop和md.nc
• 使用 Extensions → Analysis 中的工具（如 RMSD Trajectory Tool）
• 可导出数据，但不如 cpptraj 灵活

方式4：Jupyter Notebook + NGLview + MDAnalysis（高级）

适合交互式分析，但需额外安装库：

pip install nglview mdanalysis matplotlib

import MDAnalysis as mda import matplotlib.pyplot as plt from MDAnalysis.analysis import rms, rmsf u = mda.Universe('system.prmtop', 'md.nc') R = rms.RMSD(u, select="backbone") R.run() plt.plot(R.results.rmsd[:, 1], R.results.rmsd[:, 2])

四、完整案例教程：蛋白质 MD 轨迹分析与绘图

步骤1：准备文件

protein.prmtop

• md.nc（100 ns，每 10 ps 保存一帧，共 10,000 帧）

步骤2：编写 cpptraj 脚本analysis.in

parm protein.prmtop trajin md.nc 对齐骨架 rms first :1-150&!@H= RMSD rmsd out rmsd.dat first :1-150&!@H= RMSF atomicfluct out rmsf.dat :1-150&!@H= byres Rg radgyr out rg.dat :1-150 关键残基距离（如活性位点） distance d1 :45@CA :120@CA out dist_45_120.dat run

步骤3：运行分析

cpptraj -i analysis.in > analysis.log

步骤4：Python 绘图（plot_all.py）

import numpy as np import matplotlib.pyplot as plt 设置全局字体 plt.rcParams.update({'font.size': 12}) RMSD t, rmsd = np.loadtxt('rmsd.dat', unpack=True, comments=['#','@']) t_ns = t * 0.01 # 假设每帧 10 ps plt.figure() plt.plot(t_ns, rmsd) plt.xlabel('Time (ns)'); plt.ylabel('RMSD (Å)') plt.title('RMSD') plt.savefig('rmsd.png') RMSF res, rmsf = np.loadtxt('rmsf.dat', unpack=True, comments=['#','@']) plt.figure() plt.plot(res, rmsf, 'r-') plt.xlabel('Residue'); plt.ylabel('RMSF (Å)') plt.title('RMSF') plt.savefig('rmsf.png') Rg t, rg = np.loadtxt('rg.dat', unpack=True, comments=['#','@']) plt.figure() plt.plot(t*0.01, rg) plt.xlabel('Time (ns)'); plt.ylabel('Rg (Å)') plt.title('Radius of Gyration') plt.savefig('rg.png') Distance t, d = np.loadtxt('dist_45_120.dat', unpack=True, comments=['#','@']) plt.figure() plt.plot(t*0.01, d) plt.xlabel('Time (ns)'); plt.ylabel('Distance (Å)') plt.title('Distance between Res 45 and 120 CA') plt.savefig('dist_45_120.png')

运行：

python plot_all.py

五、进阶建议

1. 轨迹预处理：使用strip去除水/离子，减小文件体积

parm system.prmtop trajin md.nc strip :WAT,Na+,Cl- trajout md_noWat.nc netcdf run

1. 多重复实验：对多个轨迹分别分析，取平均±标准差
2. 自由能景观：结合 PCA 与二维直方图（np.histogram2d）
3. 自动化脚本：用 Bash/Python 批量处理多个体系

六、cpptraj 批处理与交互模式

1.CPPTRAJ- Amber官方轨迹分析工具

CPPTRAJ是AmberTools中处理坐标轨迹和数据文件的主要程序，支持批量处理和交互式两种模式。

启动方式：

bash

# 批量模式（推荐用于生产分析） cpptraj -p topology.prmtop -i analysis.in 交互模式（适合调试和测试） cpptraj -p topology.prmtop cpptraj> [命令]

七、RMSD/RMSF 的批处理写法与美化图

1. RMSD（均方根偏差）分析

衡量结构相对于参考结构的偏离程度，评估模拟稳定性。

CPPTRAJ输入文件 (rmsd.in)：

bash

# 加载拓扑和轨迹 parm protein.prmtop trajin production.nc 自动镜像处理（周期性边界条件） autoimage 计算骨架原子RMSD（以第一帧为参考） rmsd rmsd_backbone :1-300@CA,C,N out rmsd_backbone.dat time 0.1 计算全原子RMSD（以晶体结构为参考） reference crystal.pdb [crystal] rmsd rmsd_all :* reference [crystal] out rmsd_all.dat time 0.1 拟合并输出叠加后的轨迹（可选） trajout fitted.nc netcdf

Python作图 (Matplotlib)：

Python

import pandas as pd import matplotlib.pyplot as plt import numpy as np 读取CPPTRAJ输出 data = pd.read_csv('rmsd_backbone.dat', delim_whitespace=True, comment='#', names=['Frame', 'Time', 'RMSD']) fig, ax = plt.subplots(figsize=(10, 6)) 绘制RMSD曲线 ax.plot(data['Time'], data['RMSD'], color='#2E86AB', linewidth=1.5, alpha=0.8) 添加移动平均（平滑曲线） window = 50 ma = data['RMSD'].rolling(window=window).mean() ax.plot(data['Time'], ma, color='#A23B72', linewidth=2.5, label=f'{window}帧移动平均') 标注关键区域 ax.axhline(y=data['RMSD'].mean(), color='gray', linestyle='--', alpha=0.5, label=f'平均值: {data["RMSD"].mean():.2f} Å') ax.fill_between(data['Time'], data['RMSD'].mean() - data['RMSD'].std(), data['RMSD'].mean() + data['RMSD'].std(), alpha=0.2, color='gray') ax.set_xlabel('Time (ns)', fontsize=12) ax.set_ylabel('RMSD (Å)', fontsize=12) ax.set_title('Backbone RMSD vs Time', fontsize=14, fontweight='bold') ax.legend() ax.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('rmsd_analysis.png', dpi=300, bbox_inches='tight') plt.show()

2. RMSF（均方根涨落）分析

识别蛋白质柔性区域，常用于B因子计算。

CPPTRAJ输入文件 (rmsf.in)：

bash

parm protein.prmtop trajin production.nc autoimage rms reference [average] :1-300@CA,C,N 按残基计算RMSF（需先计算平均结构） average crdset MyAverage run rms ref MyAverage :1-300@CA,C,N atomicfluct out rmsf_byres.dat :1-300@CA byres 或计算B因子（与晶体学B因子比较） atomicfluct out bfactor.dat :1-300@CA byres bfactor

Python作图（热力图+柱状图）：

Python

import matplotlib.pyplot as plt import seaborn as sns 读取RMSF数据 rmsf_data = pd.read_csv('rmsf_byres.dat', delim_whitespace=True, comment='#', names=['Residue', 'RMSF']) fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8), height_ratios=[3, 1]) 上图：RMSF曲线 colors = plt.cm.RdYlBu_r(rmsf_data['RMSF'] / rmsf_data['RMSF'].max()) ax1.bar(rmsf_data['Residue'], rmsf_data['RMSF'], color=colors, width=0.8, edgecolor='none') ax1.plot(rmsf_data['Residue'], rmsf_data['RMSF'], color='black', linewidth=1, alpha=0.5) 标注高柔性区域（如loop区） threshold = rmsf_data['RMSF'].quantile(0.9) high_flex = rmsf_data[rmsf_data['RMSF'] > threshold] ax1.scatter(high_flex['Residue'], high_flex['RMSF'], color='red', s=50, zorder=5, label=f'高柔性区域 (>90%分位数)') ax1.set_ylabel('RMSF (Å)', fontsize=12) ax1.set_title('Per-Residue Root Mean Square Fluctuation', fontsize=14, fontweight='bold') ax1.legend() ax1.grid(True, alpha=0.3, axis='y') 下图：二级结构示意图（示例） 实际应用中可从DSSP分析导入 ss_data = np.random.choice([0, 1, 2], size=len(rmsf_data)) # 0=coil, 1=helix, 2=sheet ss_colors = {0: '#FFFFFF', 1: '#FF6B6B', 2: '#4ECDC4'} for i, ss in enumerate(ss_data): ax2.barh(0, 1, left=i, color=ss_colors[ss], height=0.5, edgecolor='black', linewidth=0.5) ax2.set_xlim(0, len(rmsf_data)) ax2.set_ylim(-0.5, 0.5) ax2.set_xlabel('Residue Number', fontsize=12) ax2.set_yticks([]) ax2.set_title('Secondary Structure (Red=Helix, Cyan=Sheet, White=Coil)', fontsize=10) plt.tight_layout() plt.savefig('rmsf_analysis.png', dpi=300, bbox_inches='tight')

八、氢键网络与 PCA 进阶

3. 氢键分析

分析蛋白质内部或蛋白-配体间的氢键网络。

CPPTRAJ输入文件 (hbond.in)：

bash

parm complex.prmtop trajin production.nc autoimage 蛋白质内部氢键 hbond hb_protein :1-300 out hbond_protein.dat avgout hbond_avg_protein.dat series uuseries hbonds_uu.gnu 蛋白-配体氢键（假设配体为残基301） hbond hb_ligand :1-300,301 out hbond_ligand.dat avgout hbond_avg_ligand.dat series uuseries hbonds_lig.gnu 溶剂桥接氢键 hbond hb_solvent :1-300 solventdonor :WAT solventacceptor :WAT@O out hbond_solvent.dat avgout hbond_avg_solvent.dat

Python作图（氢键时间序列）：

Python

# 读取氢键序列数据 hb_data = pd.read_csv('hbond_protein.dat', delim_whitespace=True, comment='#') fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8)) 上图：氢键数量随时间变化 total_hb = hb_data.iloc[:, 1:].sum(axis=1) # 假设第一列是时间 ax1.plot(hb_data.iloc[:, 0], total_hb, color='#2E86AB', linewidth=1.5) ax1.set_ylabel('Number of H-bonds', fontsize=12) ax1.set_title('Total Hydrogen Bonds Over Time', fontsize=14) ax1.grid(True, alpha=0.3) 下图：特定氢键存在矩阵（热力图） 选择出现频率最高的10个氢键 hb_freq = hb_data.iloc[:, 1:].mean().sort_values(ascending=False).head(10) top_hb = hb_data[hb_freq.index] im = ax2.imshow(top_hb.T, aspect='auto', cmap='Blues', interpolation='nearest') ax2.set_xlabel('Frame', fontsize=12) ax2.set_ylabel('H-bond ID', fontsize=12) ax2.set_title('Top 10 Most Frequent H-bonds Presence Matrix', fontsize=14) plt.colorbar(im, ax=ax2, label='Present (1) / Absent (0)') plt.tight_layout() plt.savefig('hbond_analysis.png', dpi=300)

4. 主成分分析 (PCA)

提取构象变化的主要模式。

CPPTRAJ输入文件 (pca.in)：

bash

parm protein.prmtop trajin production.nc autoimage rms reference :1-300@CA,C,N 协方差矩阵计算 matrix covar name covmat :1-300@CA run 对角化（runanalysis 让它在解析阶段立即执行） runanalysis diagmatrix covmat out evecs.dat vecs 10 name pc 投影到PC空间（modes 接数据集名，不是文件名） projection proj out projections.dat modes pc beg 1 end 3 :1-300@CA run

Python作图（PCA自由能景观）：

Python

from scipy.stats import gaussian_kde proj = pd.read_csv('projections.dat', delim_whitespace=True, comment='#', names=['Frame', 'PC1', 'PC2', 'PC3']) fig, axes = plt.subplots(1, 3, figsize=(18, 5)) PC1 vs PC2 自由能图 xy = np.vstack([proj['PC1'], proj['PC2']]) kde = gaussian_kde(xy) xmin, xmax = proj['PC1'].min(), proj['PC1'].max() ymin, ymax = proj['PC2'].min(), proj['PC2'].max() X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] positions = np.vstack([X.ravel(), Y.ravel()]) Z = np.reshape(kde(positions).T, X.shape) Z = -np.log(Z + 1e-10) # 转换为自由能 (kT) im1 = axes[0].contourf(X, Y, Z, levels=20, cmap='viridis') axes[0].scatter(proj['PC1'], proj['PC2'], c='white', s=1, alpha=0.3) axes[0].set_xlabel('PC1', fontsize=12) axes[0].set_ylabel('PC2', fontsize=12) axes[0].set_title('Free Energy Landscape (PC1 vs PC2)', fontsize=14) plt.colorbar(im1, ax=axes[0], label='Free Energy (kT)') PC1时间序列 axes[1].plot(proj['Frame'], proj['PC1'], color='#E63946', linewidth=1) axes[1].set_xlabel('Frame', fontsize=12) axes[1].set_ylabel('PC1', fontsize=12) axes[1].set_title('PC1 Projection Over Time', fontsize=14) axes[1].grid(True, alpha=0.3) 特征值分布（方差解释率） eigvals = pd.read_csv('evecs.dat', skiprows=1, nrows=10, delim_whitespace=True, header=None) axes[2].bar(range(1, 11), eigvals[0], color='#457B9D') axes[2].set_xlabel('Principal Component', fontsize=12) axes[2].set_ylabel('Eigenvalue', fontsize=12) axes[2].set_title('Variance Explained by PCs', fontsize=14) plt.tight_layout() plt.savefig('pca_analysis.png', dpi=300)

九、能量组分提取与多面板图

5. 能量组分提取与作图

提取能量（使用process_mdout.perl或cpptraj）：

bash

# 使用Amber工具提取 process_mdout.perl production.out 或使用cpptraj cpptraj -p topology.prmtop << EOF readdata production.out writedata energy.dat production.out[Etot] production.out[TEMP] production.out[PRESS] production.out[VOLUME] time 0.1 EOF

Python多面板能量图：

Python

energy = pd.read_csv('energy.dat', delim_whitespace=True, names=['Time', 'Etot', 'Temp', 'Press', 'Volume']) fig, axes = plt.subplots(2, 2, figsize=(14, 10)) 总能量 axes[0,0].plot(energy['Time'], energy['Etot']/1000, color='#264653', linewidth=1) axes[0,0].set_ylabel('Total Energy (10³ kcal/mol)', fontsize=11) axes[0,0].set_title('Total Energy', fontsize=12, fontweight='bold') axes[0,0].grid(True, alpha=0.3) 温度 axes[0,1].plot(energy['Time'], energy['Temp'], color='#E76F51', linewidth=1) axes[0,1].axhline(y=300, color='gray', linestyle='--', label='Target (300K)') axes[0,1].set_ylabel('Temperature (K)', fontsize=11) axes[0,1].set_title('System Temperature', fontsize=12, fontweight='bold') axes[0,1].legend() axes[0,1].grid(True, alpha=0.3) 压力 axes[1,0].plot(energy['Time'], energy['Press'], color='#2A9D8F', linewidth=1, alpha=0.7) 添加移动平均 press_ma = energy['Press'].rolling(window=100).mean() axes[1,0].plot(energy['Time'], press_ma, color='#264653', linewidth=2, label='100-pt MA') axes[1,0].set_ylabel('Pressure (bar)', fontsize=11) axes[1,0].set_xlabel('Time (ps)', fontsize=11) axes[1,0].set_title('System Pressure', fontsize=12, fontweight='bold') axes[1,0].legend() axes[1,0].grid(True, alpha=0.3) 密度（从体积计算，假设已知质量和盒子转换） 示例：简单体积图 axes[1,1].plot(energy['Time'], energy['Volume']/1000, color='#F4A261', linewidth=1) axes[1,1].set_ylabel('Volume (10³ ų)', fontsize=11) axes[1,1].set_xlabel('Time (ps)', fontsize=11) axes[1,1].set_title('System Volume', fontsize=12, fontweight='bold') axes[1,1].grid(True, alpha=0.3) plt.suptitle('Molecular Dynamics Equilibration Monitoring', fontsize=16, fontweight='bold', y=1.02) plt.tight_layout() plt.savefig('energy_components.png', dpi=300, bbox_inches='tight')

十、3D 可视化：VMD 与 PyMOL

6. VMD可视化脚本

生成VMD可视化状态的Tcl脚本：

tcl

# save_visualization.vmd # 加载结构 mol new protein.prmtop type parm7 mol addfile production.nc type netcdf first 0 last -1 step 10 waitfor all 显示设置 mol delrep 0 top mol representation NewCartoon 0.3 10.0 4.1 0 mol color ColorID 0 mol selection {protein} mol addrep top 按RMSF着色（需先导入RMSF数据） mol representation NewCartoon 0.3 10.0 4.1 0 mol color Beta mol selection {protein} mol addrep top 设置颜色范围（低RMSF=蓝，高RMSF=红） mol scaleminmax top 1 0.5 3.0 配体显示（如果有） mol representation Licorice 0.3 12.0 12.0 mol color Name mol selection {resname LIG} mol addrep top 保存图片 render TachyonInternal rmsf_colored.png

命令行渲染：

bash

vmd -e save_visualization.vmd -dispdev text

7. PyMOL高级可视化

PyMOL脚本 (visualize.pml)：

Python

# 加载结构（PyMOL 不认 prmtop，先在 cpptraj 里 trajout 导出 PDB 再加载） load protein.pdb, complex load_traj production.nc, complex 设置样式 as cartoon color marine, complex 按B因子（RMSF）着色 spectrum b, blue_white_red, minimum=0.5, maximum=3.0 显示配体 show sticks, resn LIG color atomic, resn LIG 设置视角 set_view ( 0.5, 0.3, 0.8, -0.2, 0.9, -0.4, -0.8, 0.3, 0.5, 0.0, 0.0, -50.0) 保存高质量图片 set ray_trace_mode, 1 set ray_shadows, off bg_color white ray 2400, 2400 save rmsf_pymol.png, dpi=300 制作动画（可选） mset 1 x100 mplay

十一、径向分布函数（RDF）

CPPTRAJ输入：

bash

parm system.prmtop trajin production.nc 水分子氧原子围绕蛋白质的RDF radial rdf_water.dat 0.1 10.0 :1-300 :WAT@O volume 离子围绕配体的RDF radial rdf_ion.dat 0.1 10.0 resname LIG :Na+

Python作图：

Python

rdf = pd.read_csv('rdf_water.dat', delim_whitespace=True, comment='#', names=['r', 'g(r)']) fig, ax = plt.subplots(figsize=(10, 6)) ax.plot(rdf['r'], rdf['g(r)'], color='#1D3557', linewidth=2) ax.axhline(y=1, color='gray', linestyle='--', alpha=0.5, label='Bulk density') 标注峰位 from scipy.signal import find_peaks peaks, _ = find_peaks(rdf['g(r)'], height=1.5, distance=10) ax.scatter(rdf['r'].iloc[peaks], rdf['g(r)'].iloc[peaks], color='red', s=100, zorder=5) ax.set_xlabel('r (Å)', fontsize=12) ax.set_ylabel('g(r)', fontsize=12) ax.set_title('Radial Distribution Function: Protein-Water', fontsize=14) ax.legend() ax.grid(True, alpha=0.3) plt.savefig('rdf_analysis.png', dpi=300)

十二、自动化分析流水线

完整Python自动化脚本：

Python

#!/usr/bin/env python3 """ Amber MD Trajectory Analysis Pipeline 自动生成标准分析图表 """ import subprocess import pandas as pd import matplotlib.pyplot as plt import seaborn as sns from pathlib import Path class AmberAnalyzer: def init(self, topology, trajectory, output_dir='analysis'): self.top = topology self.traj = trajectory self.outdir = Path(output_dir) self.outdir.mkdir(exist_ok=True) def run_cpptraj(self, script_name, commands): """运行cpptraj并保存输入脚本""" script_path = self.outdir / f"{script_name}.in" with open(script_path, 'w') as f: f.write(f"parm {self.top}\n") f.write(f"trajin {self.traj}\n") f.write(commands) result = subprocess.run(['cpptraj', '-i', str(script_path)], capture_output=True, text=True) return result def analyze_rmsd(self): """RMSD分析""" commands = """ autoimage rmsd rmsd_protein :1-300@CA,C,N out {outdir}/rmsd.dat time 0.1 """.format(outdir=self.outdir) self.run_cpptraj('rmsd_analysis', commands) self._plot_rmsd() def _plot_rmsd(self): """绘制RMSD图""" data = pd.read_csv(self.outdir/'rmsd.dat', delim_whitespace=True, comment='#', names=['Frame', 'Time', 'RMSD']) plt.figure(figsize=(10, 6)) plt.plot(data['Time'], data['RMSD'], color='#2E86AB', linewidth=1.5) plt.xlabel('Time (ns)') plt.ylabel('RMSD (Å)') plt.title('Protein Backbone RMSD') plt.grid(True, alpha=0.3) plt.savefig(self.outdir/'rmsd_plot.png', dpi=300, bbox_inches='tight') plt.close() 使用示例 analyzer = AmberAnalyzer('protein.prmtop', 'production.nc') analyzer.analyze_rmsd()

十三、工具选型总表

这套分析流程涵盖了Amber分子动力学模拟后处理的绝大多数需求，建议根据具体研究问题选择合适的分析组合。对于大规模轨迹，建议使用pytraj（Python接口）或并行化的CPPTRAJ以提高效率

参考资料

• Amber 官方手册（Amber24）：https://ambermd.org/doc12/Amber24.pdf
• cpptraj 教程：Tutorial C1
• Matplotlib 官网：https://matplotlib.org
• MDAnalysis 用户指南：Redirecting to https://userguide.mdanalysis.org/stable/index.html
• pytraj 文档：pytraj: data analysis package for MD simulation data

关键字

AMBER cpptraj RMSD RMSF 氢键分析 主成分分析 径向分布函数 轨迹作图

系列导航：专栏全集分子动力学模拟-Amber

更多专栏：