8.29 other inhibitors-18

受体准备

# 含一个钙离子，223个氨基酸，3221个原子
# /home/databank/zzy/project/MD/koff/trypsin/pro-pre/pro-tleap.pdb


# 0826开的这个体系5个抑制剂所用的受体一致，其中第一个做的抑制剂31是复合物晶体结构里原有的配体，其余4个是对接进去的
# 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
# ！PDBBind里的5eg4没有钙离子！实际应该是要有的，将原始结构align到我们处理好的结构上，保存为新文件。查看其文本，将钙离子信息补充到我们的pro-tleap.pdb中
# /home/databank/zzy/project/MD/koff/trypsin/pro-pre/pro-tleap.pdb

配体准备

# 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  


# 上述工作乐云师姐已完成，5个抑制剂的高斯优化已完成，拿到了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


# 我们需要晶体结构里的坐标信息。但是5eg4_ligand_native_277.mol2里的原子名字tleap不认识，因为它读取的是NEWPDB.PDB的原子名称。故需要将5eg4_ligand_native_277.mol2里的原子名称map过去。
# 85.24  /home/databank/zzy/project/MD/koff/trypsin/ligprep-opt/18
cp ../g16/18/NEWPDB.PDB ./
cp ../../complex_from_database/5eg4_ligand_18_278_dock1/5eg4_ligand_18_278_dock1.mol2 ./


# gv里分别打开 NEWPDB.PDB 与 5eg4_ligand_18_278_dock1.mol2（晶体结构配体文件）
# 两个文件下载到本地
# 本地打开GV tools-atom list-鼠标点击symbol列-rows-sort selected -拖动鼠标选择所有H-edit-delete-selected atoms），接着仍在atom list中，edit-z matrix-standardize。
# GV里 检查 两个文件的atom list，更改晶体结构的 atom list - Tag，使其Tag 与高斯优化结果完全一致,另存为crystal-gv.pdb
# 高斯优化结果保存为 g16-gv.pdb


# 85.24  /home/databank/zzy/project/MD/koff/trypsin/ligprep-opt/18
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
rm symbol coordinate
### 可视化检查（lig.pdb是否完整，位置是否与晶体结构重合，与NEWPDB.PDB序号是否一一对应）

复合物准备

# /home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-18/tleap
cp ../../../ligprep-opt/18/lig.pdb ./
cp /home/databank/zzy/project/MD/koff/trypsin/pro-pre/pro-tleap.pdb ./


conda activate ambertools
# tleap检查&更改 配体参数路径，添加离子个数
###########################################################################################
source leaprc.protein.ff19SB
source leaprc.gaff2
source leaprc.water.opc
loadamberparams /home/databank/zzy/project/MD/koff/trypsin/ligprep-opt/g16/18/lig.frcmod
loadamberprep /home/databank/zzy/project/MD/koff/trypsin/ligprep-opt/g16/18/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 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部分）

# /home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-18/tleap/gmx
#进入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 '37252 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


# 位置限制文件只是个索引，规定了哪些原子被限制
gmx_mpi genrestr -f gmx.gro -o posre1.itp
# 配体102个原子，posre文件前四行固定，从第5行开始记录哪些原子需要被限制
sed -n '1,106p' posre1.itp > posre2.itp

体系平衡

# 85.24 /home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-18/pre-equilibrium
mkdir em nvt npt


cp -r /home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-31/pre-equilibrium/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-18/smd
# 初始文件：gmx.top  md.mdp  npt.gro  smd.sh  splumed.in（dat文件模板）
scp lywu@172.21.85.24:/home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-18/pre-equilibrium/gmx.top ./
scp lywu@172.21.85.24:/home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-18/pre-equilibrium/npt.gro ./

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


source ~/.gmx2022.5-vsremd-plumed.sh
nohup sh smd.sh &


mkdir gro-check
cp ../*gro ./
# pymol 里看SMD的结果，挑选20个构象
mkdir select10 select20
# /home/databank_70t/zzy/project/koff/trypsin-18/smd/gro-check/select10
for i in {5,9,13,16,19,24,25,28,30};do cp ../smd${i}.gro ./;done
cp ../npt.gro ./
# /home/databank_70t/zzy/project/koff/trypsin-18/smd/gro-check/select20
for i in {0,2,5,7,8,9,11,13,14,15,16,19,21,22,24,25,28,30,35};do cp ../smd${i}.gro ./;done
cp ../npt.gro ./

vsREMD

85.23 /home/databank_70t/zzy/project/koff/trypsin-18/vsremd
mkdir 10replica 20replica

10replica

# /home/databank_70t/zzy/project/koff/trypsin-18/vsremd/10replica/smd-select
ls > index
for i in {0..9};do a=`expr ${i} + 1`; b=`sed -n "${a}p" index`;mv $b smd-select${i}.gro;done


cp /home/databank_70t/zzy/project/koff/trypsin-31/vsremd/temp ./
cp /home/databank_70t/zzy/project/koff/trypsin-31/vsremd/md.mdp ./


# 制作index文件：Protein_MOL_CA 和 Water_Cl-
gmx_mpi make_ndx -f smd-select/smd-select0.gro -o index.ndx
> 1|15|14
> 17|16


# 准备好文件，暂时不生成tpr-mdp参数中 运行时间可能会根据资源有所修改
#  nohup sh -c 'for i in {0..9};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'


# 4090 172.21.75.1 /home/dddc/zzy/project/koff/trypsin-18/vsremd
scp -r dddc@172.21.85.23:/home/databank_70t/zzy/project/koff/trypsin-18/vsremd/10replica ./
source ~/.gmx-2022.5-vsremd.sh
nohup mpirun -np 60 gmx_mpi mdrun -v -deffnm md -multidir md0 md1 md2 md3 md4 md5 md6 md7 md8 md9 -replex 1000 &


(openmpi) dddc@gpu-4090:~/zzy/project/koff/trypsin-18/vsremd/10replica$ nohup mpirun -np 60 gmx_mpi mdrun -v -deffnm md -multidir md0 md1 md2 md3 md4 md5 md6 md7 md8 md9 -replex 1000 &
[1] 210317
step 4100, will finish Thu Sep  5 09:57:37 2024
(openmpi) dddc@gpu-4090:~/zzy/project/koff/trypsin-18/vsremd/10replica$ date
Fri Aug 30 04:31:49 AM UTC 2024
交了200ns的任务，跑到25000000步杀掉

20replica

# /home/databank_70t/zzy/project/koff/trypsin-18/vsremd/20replica/smd-select
ls > index
for i in {0..19};do a=`expr ${i} + 1`; b=`sed -n "${a}p" index`;mv $b smd-select${i}.gro;done


cp /home/databank_70t/zzy/project/koff/trypsin-31/vsremd/temp ./
cp /home/databank_70t/zzy/project/koff/trypsin-31/vsremd/md.mdp ./


# 制作index文件：Protein_MOL_CA 和 Water_Cl-
gmx_mpi make_ndx -f smd-select/smd-select0.gro -o index.ndx
> 1|15|14
> 17|16


# 准备好文件，暂时不生成tpr-mdp参数中 运行时间可能会根据资源有所修改
#  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 ../select20/smd-select${i}.gro -c ../select20/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 &


(base) [lywu@localhost 20replica]$ 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 &
[1] 125618

sleep（可选）

# 可选
############################"sleep.sh" 3L, 195C #############################################
source ~/.gmx2022.5-vsremd.sh
sleep 6h
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
#############################################################################################