如果是在自己的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-tmpMD全流程脚本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 # 选6tip3p是经典三点式水分子表示适合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文件如A1-100res B101-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