Amber分子动力学模拟16: Amber轨迹分析与作图-1

本文收录于专栏 分子动力学模拟-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 → Load system.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

五、进阶建议

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

六、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

更多专栏：