在云服务器AutoDL实现分子动力学模拟全流程

如果是在自己的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这种LBlackWell架构的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文件及数据处理脚本分享）

# 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 &