如果是在自己的ubuntu虚拟机中进行MD模拟，跳过前面的AutoDL部分即可。以下内容为本人使用经验，希望为做MD的朋友推荐高效的计算资源和简易可复现的workflow，无商业行为。

AutoDL是国内先进的GPU算力云平台，主要面向深度学习和人工智能开发，提供高性能GPU租赁及即开即用环境。在AutoDL中RTX4090的租赁价格约2元/小时（视不同地区服务器群配置而有所差异），在行业内极具性价比，并且4090可为GROMACS生物大分子MD模拟提供充足的算力，是运行MD的首选。

GPU开机操作及镜像选择

注册和充值相关事宜参考官方“帮助”文件

首先进入算力市场，选择和自己所在城市距离较近的机房，选择RTX4090或RTX4090D GPU（计算资源较为紧张，夜间容易抢到），点击“n卡可租”，AutoDL为我们提供了现成的GROMACS镜像，如图选择镜像。另外如果MD时间在100ns以上建议扩容数据盘

这里讲解一下为什么选择4090，首先4090的算力足够完成模拟，并且GROMACS的23年版本比较旧，cuda与5090这种BlackWell架构的GPU不兼容，另外该版本本身也没有对BlackWell、Ada架构的单元做适配，亲测该任务上4090是A100-PCIE-40GB计算速度的两倍。如果要使用5090大家需要在服务器内部自行配置最新版GROMACS，也欢迎大家在评论区分享配置流程和自己的镜像。

开机后点击右上角控制台，点击该机器的JupyterLab按钮，进入界面后点击左侧的“autodl-tmp”，再点击右侧“终端“这是数据盘地址，IO速度比较快，我们尽量在这里完成模拟流程。如果直接点击了终端，输入cd /root/autodl-tmp

MD全流程脚本

CHARMM36、Amber99SB-ILDN力场都是同样的流程，只是mdp文件的参数不同，大家可以理解为mdp是个描述设置的文件。MD流程包括：pdb2gmx转化——定义box——注入溶剂环境——能量最小化——NVT温度平衡——NPT压力平衡——MD正式模拟。运行时请确保pdb文件、ions.mdp、minim.mdp、nvt.mdp、npt.mdp、md.mdp六个文件在工作目录下，python脚本及所有mdp文件可在以下两篇文章获取（Amber99SB-ILDN力场MD模拟mdp文件及数据处理脚本分享）（CHARMM36力场MD模拟相关mdp文件）

# Run this for A_B.pdb
gmx pdb2gmx -f A_B.pdb -o A_B_processed.gro -water tip3p -ignh 
# 选6，tip3p是经典三点式水分子表示，适合Amber99SB-ILDN

# Define box
gmx editconf -f A_B_processed.gro -o A_B_box.gro -c -d 1.0 -bt dodecahedron 
# 十二面体，距边界1纳米

# Add water
gmx solvate -cp A_B_box.gro -cs spc216.gro -o A_B_solv.gro -p topol.top 
# spc216是水的标准格式

# Compile the run input file for adding ions
gmx grompp -f ions.mdp -c A_B_solv.gro -p topol.top -o ions.tpr 
# 可视为一个打包机制，需要先写mdp

# Add K+ and Cl- to 0.14 M, and neutralize
gmx genion -s ions.tpr -o A_B_ions.gro -p topol.top -pname K -nname CL -neutral -conc 0.14
# 卵母细胞环境高K低Na 选13 SOL，将部分水分子替换成离子，切忌选protein

# Assuming you have a minim.mdp file
gmx grompp -f minim.mdp -c A_B_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em -ntmpi 1 -ntomp 16 -gpu_id 0
# 还是先有mdp再打包，目的是移开overlap部分确保作用力Fmax小于1000（尽量一开始pdb就删无序区）
# -ntomp 16是CPU线程数，一般4090是16，如果报错就按显示的填

# NVT (Temperature equilibration - e.g., 100 ps)
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr 
# 升到生理温度310K，先验步 如果有warning需要确认无害后添加参数 -maxwarn 1
gmx mdrun -v -deffnm nvt

# NPT (Pressure equilibration - e.g., 100 ps)
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr 
# 升到一个大气压，等压过程一定会改变box体积，该步对overlap敏感
gmx mdrun -v -deffnm npt
# 同上，注释： -f参数 -c主体坐标 -r约束 -t状态

# Production run (e.g., 100 ns)
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr
gmx mdrun -v -deffnm md
# 正式模拟，耗时较长

# GROMACS为了防止原子跑到体系外，有一个机制，如果原子从一端出了盒子会从另一端进入
# 因此如果直接分析，被分开的原子会影响rmsd等数据，该步将被盒子分开的原子拼接回原位 选0
gmx trjconv -s md.tpr -f md.xtc -o md_noPBC.xtc -pbc mol -ur compact
# 再将蛋白居中 选1 再选0
gmx trjconv -s md.tpr -f md_noPBC.xtc -o md_final.xtc -center -pbc none
# 该步骤产生三个相近大小的轨迹文件.xtc 因此数据盘可以扩容

# 基本分析，.xvg文件可以用python做图
gmx energy -f md.edr -o energy.xvg
gmx rms -f md.xtc -s md.tpr -o rmsd.xvg -tu ns
gmx rmsf -f md.xtc -s md.tpr -o rmsf.xvg -res
gmx hbond -f md.xtc -s md.tpr -num hbond.xvg
gmx gyrate -f md.xtc -s md.tpr -o gyrate.xvg

# 互作分析，系统并不知道哪个res是哪个蛋白的，我们需要给它们一个index
# 由于GROMACS对蛋白的编号顺序与Chimerax中看到的可能不同，所以建议先获取最后一帧观察一下
gmx trjconv -f md.xtc -s md.tpr -o last_frame.pdb -dump 100000 # 假如100ns情况下

# 创建index文件（如A：1-100res B：101-300res）进入后应该会有16个现有标签
gmx make_ndx -f md.tpr -o index.ndx # 如果已有index.ndx可以写 -n index.ndx
ri 1-100
name 17 A
ri101-300
name 18 B


# 我个人很喜欢的一种分析，判断A-B间的氢键数量，先选17再选18
gmx hbond -f md.xtc -s md.tpr -n index.ndx -num hbonds_interface.xvg

# AutoDL镜像自带相关包，ubuntu系统需要自行下载pandas等python包
python figure_hbonds_interface.py

# contacts位点数，先选17再选18
gmx mindist -f md.xtc -s md.tpr -n index.ndx -on num_contacts.xvg -d 0.6 # 6埃截断

# 处理contacts的python
python plot_contact.py

# sasa互作面积，整体求分别求再相减除二
gmx sasa -f md.xtc -s md.tpr -n index.ndx -o total_sasa.xvg
gmx sasa -f md.xtc -s md.tpr -n index.ndx -o sasa_A.xvg
gmx sasa -f md.xtc -s md.tpr -n index.ndx -o sasa_B.xvg

# python
python interface_area.py


# 服务器相关操作，后台运行、拼接轨迹等
#续跑
nohup gmx mdrun -v -deffnm md -cpi md.cpt -noappend &
# 合并
gmx trjcat -f md.xtc md.part0002.xtc -o full_100ns.xtc -cat
# 能量文件合并
gmx eneconv -f md.edr md.part0002.edr -o full_md.edr
# 看运行状态
tail -f md.part0002.log     ps aux | grep gmx

不建议后台运行，很容易断，如果像上面我的脚本一样，请直接关闭Jupyter，不要关闭终端，如果要后台运行请自己写一下nohup &

如果我的内容有所帮助，还望友友们点赞收藏支持一下～关注我一起探讨更多生物信息学内容