准备好对接后的蛋白pro.pdb和配体环肽文件test.mo第一步python 00_make_single_lig.py test.mol2 test_LIG_raw.mol200_make_single_lig.py文件如下#!/usr/bin/env python3 import sys if len(sys.argv) ! 3: print(Usage: python 00_make_single_lig.py input.mol2 output.mol2) sys.exit(1) infile sys.argv[1] outfile sys.argv[2] section None with open(infile, r) as fin, open(outfile, w) as fout: for line in fin: if line.startswith(TRIPOS): section line.strip() fout.write(line) if section TRIPOSSUBSTRUCTURE: fout.write( 1 LIG 1 GROUP 0 **** **** 0 ROOT\n) continue if section TRIPOSATOM: parts line.split() if len(parts) 9: atom_id int(parts[0]) atom_name parts[1] x float(parts[2]) y float(parts[3]) z float(parts[4]) atom_type parts[5] charge float(parts[8]) fout.write( f{atom_id:7d} {atom_name:8s} f{x:10.4f} {y:10.4f} {z:10.4f} f{atom_type:8s} {1:5d} {LIG:8s} {charge:10.6f}\n ) else: fout.write(line) elif section TRIPOSSUBSTRUCTURE: continue else: fout.write(line)第二步antechamber由于环肽原子量很多这一步要持续几个小时移植到能量收敛到0.0005.antechamber -i test_LIG_raw.mol2 -fi mol2 \ -o LIG.mol2 -fo mol2 \ -at gaff2 -c bcc -nc 0 \ -rn LIG -s 2第三步parmchk2 -i LIG.mol2 -f mol2 -o LIG.frcmod -s gaff2第四步处理蛋白pdb4amber -i pro.pdb -o pro_amber.pdb --reduce --dry第五步继续处理单边加氢等awk /^ATOM/ || /^HETATM/ { atomsubstr($0,13,4) gsub(/ /,,atom) if (atom ~ /^H/) next } {print} pro_amber.pdb pro_amber_noH.pdb第六步将配体受体合并tleap -f 01_tleap_complex.in | tee 01_tleap_complex.log01_tleap_complex.in文件如下source leaprc.protein.ff19SB source leaprc.gaff2 source leaprc.water.opc loadamberparams LIG.frcmod pro loadpdb pro_amber_noH.pdb lig loadmol2 LIG.mol2 complex combine {pro lig} check complex charge complex solvateoct complex OPCBOX 10.0 addionsrand complex Na 0 addionsrand complex Cl- 0 saveamberparm complex complex.prmtop complex.inpcrd savepdb complex complex_solvated.pdb quit第八步就是正常MD流程1. 限制性最小化02_min1.inMinimization 1: solvent and ions, solute restrained cntrl imin1, maxcyc10000, ncyc5000, cut10.0, ntb1, ntr1, restraint_wt10.0, restraintmask!:WAT,Na,Cl-, / END2. 全体系最小化 03_min2.inMinimization 2: whole system cntrl imin1, maxcyc10000, ncyc5000, cut10.0, ntb1, / END3. 升温04_heat.inHeating from 0 to 300 K cntrl imin0, irest0, ntx1, ntb1, cut10.0, nstlim250000, dt0.002, tempi0.0, temp0300.0, ntc2, ntf2, ntt3, gamma_ln2.0, ntpr1000, ntwx5000, ntwr5000, ntr1, restraint_wt5.0, restraintmask!:WAT,Na,Cl-, nmropt1, / wt TYPETEMP0, istep10, istep2250000, value10.0, value2300.0 / wt TYPEEND / END5. NPT 平衡有弱限制 05_eq1.inNPT equilibration with weak restraint cntrl imin0, irest1, ntx5, ntb2, ntp1, taup2.0, cut10.0, nstlim500000, dt0.002, temp0300.0, ntc2, ntf2, ntt3, gamma_ln2.0, ntpr1000, ntwx5000, ntwr5000, ntr1, restraint_wt1.0, restraintmask!:WAT,Na,Cl-, / END6. NPT 平衡无限制 06_eq2.inNPT equilibration without restraint cntrl imin0, irest1, ntx5, ntb2, ntp1, taup2.0, cut10.0, nstlim500000, dt0.002, temp0300.0, ntc2, ntf2, ntt3, gamma_ln2.0, ntpr1000, ntwx5000, ntwr5000, / END8. 生产模拟单段 10 ns 07_prod.inProduction MD cntrl imin0, irest1, ntx5, ntb2, ntp1, taup2.0, cut10.0, nstlim5000000, dt0.002, temp0300.0, ntc2, ntf2, ntt3, gamma_ln2.0, ntpr5000, ntwx5000, ntwr50000, ioutfm1, ntxo2, / END9.执行运行命令run_md.sh#!/bin/bash set -e PMEMDpmemd.cuda echo Step 1: restrained minimization $PMEMD -O \ -i 02_min1.in \ -o 02_min1.out \ -p complex.prmtop \ -c complex.inpcrd \ -r 02_min1.rst7 \ -ref complex.inpcrd echo Step 2: full minimization $PMEMD -O \ -i 03_min2.in \ -o 03_min2.out \ -p complex.prmtop \ -c 02_min1.rst7 \ -r 03_min2.rst7 echo Step 3: heating $PMEMD -O \ -i 04_heat.in \ -o 04_heat.out \ -p complex.prmtop \ -c 03_min2.rst7 \ -r 04_heat.rst7 \ -x 04_heat.nc \ -ref 03_min2.rst7 echo Step 4: NPT equilibration with restraint $PMEMD -O \ -i 05_eq1.in \ -o 05_eq1.out \ -p complex.prmtop \ -c 04_heat.rst7 \ -r 05_eq1.rst7 \ -x 05_eq1.nc \ -ref 04_heat.rst7 echo Step 5: NPT equilibration without restraint $PMEMD -O \ -i 06_eq2.in \ -o 06_eq2.out \ -p complex.prmtop \ -c 05_eq1.rst7 \ -r 06_eq2.rst7 \ -x 06_eq2.nc echo Step 6: production MD prev06_eq2.rst7 for i in $(seq -w 1 10) do echo Production segment ${i} $PMEMD -O \ -i 07_prod.in \ -o prod_${i}.out \ -p complex.prmtop \ -c ${prev} \ -r prod_${i}.rst7 \ -x prod_${i}.nc \ -inf prod_${i}.info prevprod_${i}.rst7 done echo MD finished.最后运行chmod x run_md.sh ./run_md.sh