clean version
配体准备
- # 22224
- # 172.21.85.24
- # 使用薛定谔优化结构（图形化界面操作）
- # pH = 8.0，手性从3D结构定（晶体结构），力场为2017版的默认力场OPLS3
- # 高斯优化
- # 原始文件 /home/databank/lywu/trypsin/inhibitiors/gaussian_job
- # 注意检查体系的电荷和自选多重度
- antechamber -i lig.sdf -fi sdf -o lig.gjf -fo gcrt -at gaff2 -gn "%nproc=4" -gm "%mem=4GB" -gk "#B3LYP/6-31G* em=gd3bj pop=MK iop(6/33=2,6/42=6) opt" -rn MOL -nc 0
- # 上述工作乐云师姐已完成，高斯优化已完成，拿到了log文件
- ###################################################################################
- # 以下工作由我接手
- # 85.24 /home/databank/zzy/project/MD/koff/trypsin/ligprep-opt
- # 基于高斯优化日志文件，拟合RESP电荷，生成AMBER需要的参数文件（prep)
- conda activate ambertools
- for i in
`ls`;do cd $i;antechamber -i lig.log -fi gout -o lig.prep -fo prepi -c resp;parmchk2 -i lig.prep -f prepi -o lig.frcmod -a y;cd ../;done - # 这一步本质上是获得坐标，尝试直接从晶体结构中抠出小分子，提取坐标
- # gv里分别打开 NEWPDB.PDB 与 5eg4_ligand_native_277.mol2（晶体结构配体文件）
- # 打开GV tools-atom list-鼠标点击symbol列-rows-sort selected -拖动鼠标选择所有H-edit-delete-selected atoms），接着仍在atom list中，edit-z matrix-standardize，分别保存<步骤一>。
- # 得到 g16-gv.pdb crystal-gv.pdb
- # GV里 检查 g16-gv.pdb crystal-gv.pdb，更改crystal-gv.pdb （atom list - Tag)，使其Tag 与g16-gv.pdb 完全一致。
- # 20240826 经乐云师姐指正，参考子涧师兄20230703的提出的方法，<步骤一>完成后，两个文件的原子index是一致的，无需再对（atom list - Tag)做修改。
- # 注意：输入的文件晶体结构文件需要是sdf文件！mol2无法成功
- dos2unix *
- grep -v "CONECT" g16-gv.pdb |cut -b 1-26 > symbol
- cat crystal-gv.pdb | cut -b 29-78|sed '1,2d'|sed '1i Combineing the coordinates from the crystal structure with the atomic numbers from the Gaussian optimized structure'> coordinate
- paste -d " " symbol coordinate|sed '1i TITLE lig.pdb'> lig.pdb
- ### 可视化检查
受体准备
- # PDBBIND
- # /home/databank/zzy/project/MD/koff/trypsin/pro-pre/pdbbind
- conda activate pdb2pqr
- pdb2pqr30 5eg4_ligand_native_277_protein.pdb 5eg4_pred.pqr --ff AMBER --ffout AMBER --with-ph 8.0 --pdb-output 5eg4_pred.pdb
- # 直接从PDBBIND蛋白结构出发，用tleap读取，存成新的蛋白文件，再进行二硫键相关的处理
- # /home/databank/zzy/project/MD/koff/trypsin/pro-pre
- tleap -f tleap-pdbclean.in
- # 得到 5eg4-pdb2pqr-tleap.pdb
- # 二硫键的特殊处理
- 1）找到形成二硫键的氨基酸，把残基名字改为CYX；
- 2）删掉CYX的所有H；
- 3）用tleap中的bond连接二硫键，可以在cpptraj中的trajout output.pdb conect命令中输出CONET信息检查
- grep -v "H CYX" 5eg4-pdb2pqr-tleap.pdb|grep -v "HA CYX"|grep -v "HB2 CYX"|grep -v "HB3 CYX" >pro-tleap.pdb
- (ambertools) [lywu@localhost pro-pre]$ wc -l 5eg4_pred.pdb 5eg4-pdb2pqr-tleap.pdb pro-tleap.pdb
- 3221 5eg4_pred.pdb
- 3222 5eg4-pdb2pqr-tleap.pdb
- 3174 pro-tleap.pdb
- # 0827 加入钙离子
- # /home/databank/zzy/project/MD/koff/trypsin/pro-pre
- # 参考 /home/databank/lywu/trypsin/com/tleap/pro-tleap.pdb 最后几行，在最后添加钙离子信息
复合物准备
- # /home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-31
- # tleap
- source leaprc.protein.ff19SB
- source leaprc.gaff2
- source leaprc.water.opc
- loadamberparams /home/databank/zzy/project/MD/koff/trypsin/ligprep-opt/g16/31/lig.frcmod
- loadamberprep /home/databank/zzy/project/MD/koff/trypsin/ligprep-opt/g16/31/lig.prep
- mol=loadpdb lig.pdb
- pro=loadpdb pro-tleap.pdb
- bond pro.7.SG pro.137.SG
- bond pro.25.SG pro.41.SG
- bond pro.173.SG pro.197.SG
- bond pro.162.SG pro.148.SG
- bond pro.116.SG pro.183.SG
- bond pro.109.SG pro.210.SG
- com=combine {pro mol}
- saveamberparm pro pro.prmtop pro.inpcrd
- saveamberparm mol lig.prmtop lig.inpcrd
- saveamberparm com native.prmtop native.inpcrd
- savepdb com com-dry.pdb
- charge com
- solvatebox com OPCBOX 10
- addionsrand com Na+ 38 Cl- 48
- 和师姐讨论后，改为 addionsrand com Cl- 10
- charge com
- saveamberparm com complex.prmtop complex.inpcrd
- savepdb com com.pdb
- # 二硫键的特殊处理
- 1）找到形成二硫键的氨基酸，把残基名字改为CYX；
- 2）删掉CYX的所有H；
- 3）用tleap中的bond连接二硫键，可以在cpptraj中的trajout output.pdb conect命令中输出CONET信息检查
- # /home/databank/zzy/project/MD/koff/trypsin/complex-pre/31/tleap/check-SS-bond
- cpptraj
- parm ../native.prmtop
- trajin ../com-dry.pdb
- trajout test.pdb conect
- run
- # 是检查prmtop文件的键连信息
- # 除此以外，也可以检查topol文件的键连信息（bond部分）
- #进入python
- import parmed as pmd
- amber = pmd.load_file('complex.prmtop','complex.inpcrd')
- amber.save('gmx.top')
- amber.save('gmx.gro')
- #修改gmx.top 与 gmx.gro
- sed -i "s/WAT/SOL/g" gmx.top
- sed -i "s/system1/Protein/g" gmx.top
- # 0827 冲！
- # 在蛋白质结束和小分子开始前添加位置限制（moleculetype）
- sed -i '36213 a ; Include Position restraint file\n#ifdef POSRES\n#include "posre1.itp\n#endif\n' gmx.top
- # 在小分子结束和离子开始前添加位置限制（moleculetype）
- sed -i '37443 a ; Include Position restraint file\n#ifdef POSRES\n#include "posre2.itp\n#endif\n' gmx.top
- sed -i "s/WAT/SOL/g" gmx.gro
- # /home/databank_70t/zzy/project/koff/trypsin-31/tleap/gmx
- # Protein_MOL_CA Water_Cl-
- gmx_mpi make_ndx -f gmx.gro -o index.ndx
- > 1|15|14
- > 17|16
- # 位置限制文件只是个索引，规定了哪些原子被限制
- gmx_mpi genrestr -f gmx.gro -o posre1.itp
- sed -n '1,62p' posre1.itp > posre2.itp
体系平衡
- # 85.24 /home/databank/zzy/project/MD/koff/trypsin/complex-pre/31/pre-equilibrium
- # 85.23 /home/databank_70t/zzy/project/koff/trypsin-31/pre-equilibrium
- mkdir em nvt npt mdp
- cp /home/databank/lywu/vsremd_plus/protein-ligand/trypsin-benzamidine/pre-equilibrium/mdp-file/* ./mdp/
- cp ../tleap/gmx/* ./
- ### 能量最小化
- gmx_mpi grompp -f ../mdp/em.mdp -c ../gmx.gro -r ../gmx.gro -p ../gmx.top -o em.tpr
- nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm em &
- ### 等温等容预平衡
- # 修改为1fs
- gmx_mpi grompp -f ../mdp/nvt.mdp -c ../em/em.gro -r ../em/em.gro -p ../gmx.top -o nvt.tpr
- nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm nvt &
- # 修改为2fs
- gmx_mpi grompp -f ../mdp/nvt.mdp -c nvt.gro -r nvt.gro -p ../gmx.top -o nvt2.tpr
- nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm nvt2 &
- ### 等温等压预平衡
- # 修改为1fs
- gmx_mpi grompp -f ../mdp/npt.mdp -c ../nvt/nvt2.gro -r ../nvt/nvt2.gro -p ../gmx.top -o npt.tpr
- nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm npt &
- # 修改为2fs
- gmx_mpi grompp -f ../mdp/npt.mdp -c npt.gro -r npt.gro -p ../gmx.top -o npt2.tpr
- nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm npt2 &
SMD
- # 24上的SMD很慢，去23上跑
- # 85.23 /home/databank_70t/zzy/project/koff/trypsin-31/smd
- scp -r lywu@172.21.85.24:/home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-31 ./
- cp ../pre-equilibrium/npt/npt.gro ./
- cp ../pre-equilibrium/topol.top ./gmx.top
- 在npt.gro里找到对应的原子index,修改 dat文件里的index
- # 如上图，师姐的dat里，171ASP是一样的，只需要改配体的index
- sed -i 's/3225/3328/g' smd.*
- # 检查 smd.sh 的调用核数
- (base) [dddc@localhost smd]$ nohup sh smd.sh &
- [1] 89571 (done)
- # 苯甲脒 在SMD中翻了个身出去了，但是分子本身变化不是很明显，增加副本试一下
- scp lywu@172.21.85.24:/home/databank/lywu/vsremd_plus/protein-ligand/trypsin-benzamidine/test/splumed.in ./
- for i in {0..39};do m=
`echo $i|awk '{print 0.20+$i*0.10}'`; sed "s/temp1/${m}/g" splumed.in > smd.${i}.dat;sed -i "s/COLVAR/COLVAR.${i}/g" smd.${i}.dat; done - (base) [dddc@localhost 40loop]$ nohup sh smd.sh &
- [1] 212774 (done)
- # /home/databank_70t/zzy/project/koff/with_NA_0827_trypsin-31/smd/40loop-mol-centerofmass
- # /home/databank_70t/zzy/project/koff/trypsin-31/smd
- # 40个SMD的结果看了，抑制剂31被拽着头（苯甲脒）拉出去了，怀疑这个不合理，经讨论，尝试将ligand feather 改为 质心，即抑制剂的全原子index
- # 初始文件：gmx.top md.mdp npt.gro smd.sh splumed.in（dat文件模板）
- for i in {0..39};do m=
`echo $i|awk '{print 0.20+$i*0.10}'`; sed "s/temp1/${m}/g" splumed.in > smd.${i}.dat;sed -i "s/COLVAR/COLVAR.${i}/g" smd.${i}.dat; done - (ambertools) [dddc@localhost smd]$ nohup sh smd.sh &
- [1] 214342
vsREMD
- # pymol 里看SMD的结果，挑选20个构象
- # 85.23 /home/databank_70t/zzy/project/koff/trypsin-31/smd/40loop/gro-check/vsremd-select
- for i in {0,10,13,15,17,21,24,25,27,28,29,31,32,33,34,35,36,37,38,39};do cp ../smd${i}.gro ./;done
- ls > index
- for i in {1..20};do b=
`sed -n "${i}p" index`;mv $b smd-select${i}.gro;done - # 到24 上跑vsREMD 100个核
- # 85.24 /home/databank/zzy/project/MD/koff/trypsin/vsremd/trypsin-31/smd-select
- for i in {0..19};do a=
`expr ${i} + 1`;mv smd-select${a}.gro smd-select${i}.gro;done - # /home/databank/zzy/project/MD/koff/trypsin/vsremd/trypsin-31
- nohup sh -c 'for i in {0..19};do a=
`expr ${i} + 1`; b=`sed -n "${a}p" temp`; mkdir md$i; sed "s/300.00/${b}/g" md.mdp > md$i/md.mdp; cd md$i; gmx_mpi grompp -f md.mdp -p ../gmx.top -n ../index.ndx -r ../smd-select/smd-select${i}.gro -c ../smd-select/smd-select${i}.gro -o md.tpr; cd ../; done' - source ~/.gmx2022.5-vsremd.sh
- nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm md -multidir md0 md1 md2 md3 md4 md5 md6 md7 md8 md9 md10 md11 md12 md13 md14 md15 md16 md17 md18 md19 -replex 1000 &
- # 85.24 资源不足，放到4090（75.1）上跑
- # cd
- source ~/.gmx-2022.5-vsremd.sh
- # 资源有限，想在周五拿到初步结果，只跑了10个副本
- nohup mpirun -np 50 gmx_mpi mdrun -v -deffnm md -multidir md0 md1 md2 md3 md4 md5 md6 md7 md8 md9 -replex 1000 &