constant velocity SMD

轨迹可视化

# Trypsin 体系需要生成Protein_CA_MOL的index
gmx_mpi make_ndx -f smd0.gro -o index.ndx
> 1|14|15
21 Protein_CA_MOL:  3338 atoms
# Hiv 体系
gmx_mpi make_ndx -f smd0.gro -o index.ndx
> 1|13
18 Protein_MOL  :  3242 atoms


# 轨迹处理（去周期性，去水）
echo 18|gmx_mpi trjconv -f ../smd1.xtc -s ../smd1.tpr -n index.ndx -o dry1.xtc -pbc mol -ur compact


# 生成参考构象
gmx_mpi trjconv -f dry1.xtc -s ../smd1.tpr -n index.ndx -o dry.pdb -dump 0


# cpptraj 将xtc转为pdb，即可可视化观察轨迹

SMD 的问题

# 20240913 乐云师姐发现：
之前 plumed 的所谓的smd, 因为设置了asp-配体 距离为cv,在smd进行的过程中，ASP发生了剧烈的形变，这使得构象不合理，从而导致了配体很快的解离

plumed

# plumed SLOPE/KAPPA 参数的设置
/home/databank_70t/zzy/project/koff/hiv/hiv-B409/smd-test/
# 测试下来发现 plumed 没办法做恒速的smd,只能恒力，因此转向gmx自带的方案

gmx-smd

# 22224 上测试 gmx 自带的 smd
# smd 只需要在mdp里加一些参数即可，先测试下师姐之前写的参数
# gmx-smd 有一个问题是水盒子需要扩大，10A会报错

# 故重新准备体系，水盒子先设置到15A，还要预平衡
# /home/databank/zzy/project/MD/koff/hiv/hiv-B409/smd/0913-constant-velocity/tleap
# tleap
source leaprc.protein.ff19SB
source leaprc.gaff2
source leaprc.water.opc
loadamberparams /home/databank/zzy/project/MD/koff/hiv/hiv-B409/sys-pre/ligand/tleap-input/antechamber/lig.frcmod
loadamberprep /home/databank/zzy/project/MD/koff/hiv/hiv-B409/sys-pre/ligand/tleap-input/antechamber/lig.prep
mol=loadpdb lig.pdb 
pro=loadpdb 1ec1-pdb2pqr-tleap.pdb
com=combine {pro mol}
savepdb com com-dry.pdb
charge com
solvatebox com OPCBOX 15.0
addionsrand com Cl- 6
charge com
saveamberparm com complex.prmtop complex.inpcrd
savepdb com com.pdb

# amber参数转为gmx参数，并略作修改为gmx认的得的标准化版本，添加位置限制信息
mkdir gmx
cd gmx 


python 
import parmed as pmd
amber = pmd.load_file('../complex.prmtop','../complex.inpcrd')
amber.save('gmx.top')
amber.save('gmx.gro')


sed -i '17597 a ; Include Position restraint file\n#ifdef  POSRES\n#include "posre1.itp\n#endif\n' gmx.top
sed -i '18603 a ; Include Position restraint file\n#ifdef  POSRES\n#include "posre2.itp\n#endif\n' gmx.top


# 位置限制文件只是个索引，规定了哪些原子被限制
gmx_mpi genrestr -f gmx.gro -o posre1.itp 
sed -n '1,108p' posre1.itp > posre2.itp

体系平衡

# 85.24 /home/databank/zzy/project/MD/koff/hiv/hiv-B409/smd/0913-constant-velocity/pe


mkdir em nvt npt mdp
cp -r ../../../pre-equilibrium/mdp ./


### 能量最小化
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 &


gmx_mpi grompp -f ../mdp/npt.mdp -c ../em/em.gro -r ../em/em.gro -p ../gmx.top -o npt.tpr
nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm npt &


# nvt,npt均从很短的，不报错的步数开始，循环以下步骤，直至nvt 1fs 能跑100ps
# 85.24 /home/databank/zzy/project/MD/koff/hiv/hiv-B409/smd/0913-constant-velocity/pe/loop-till-nvt1fs-done
##################################################################################
gmx_mpi grompp -f ../mdp/em.mdp -c npt.gro -r npt.gro -p ../gmx.top -o em.tpr
nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm em &


gmx_mpi grompp -f ../mdp/nvt.mdp -c em.gro -r em.gro -p ../gmx.top -o nvt.tpr
nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm nvt &


gmx_mpi grompp -f ../mdp/em.mdp -c npt.gro -r npt.gro -p ../gmx.top -o em.tpr
nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm em &


gmx_mpi grompp -f ../mdp/npt.mdp -c em.gro -r em.gro -p ../gmx.top -o npt.tpr
nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm npt &
##################################################################################


# 提交 2fs 的nvt任务
# /home/databank/zzy/project/MD/koff/hiv/hiv-B409/smd/0913-constant-velocity/pe/nvt
gmx_mpi grompp -f ../mdp/nvt.mdp -c ../loop-till-nvt1fs-done/nvt.gro -r ../loop-till-nvt1fs-done/nvt.gro -p ../gmx.top -o nvt.tpr
nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm nvt &


### 等温等压预平衡
# 修改为1fs
gmx_mpi grompp -f ../mdp/npt.mdp -c ../nvt/nvt.gro -r ../nvt/nvt.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

# 修改 smd 参数里的 group1，生成相应的Index
gmx_mpi make_ndx -f npt2.gro -o index.ndx


# 跑
gmx_mpi grompp -f smd.mdp -p gmx.top -c npt2.gro -r npt2.gro -n index.ndx -o smd.tpr
nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm smd &


# 检查轨迹
gmx_mpi trjconv -f ../../test2-slower/smd.xtc -s ../../test2-slower/smd.tpr -n index.ndx -o dry.xtc -pbc mol -ur compact


gmx_mpi trjconv -f ../../test2-slower/smd.xtc -s ../../test2-slower/smd.tpr -n index.ndx -o dry.xtc -pbc mol -ur compact


# 结果不错，能跑出去，且 asp171动的不是很明显
# 先提交5个副本/10个副本的测试试一下
# 0913晚上的测试中。没有放解离构象，是无效的测试
# 0914 先对smd结果做进一步的分析 （85.24）
# /home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-31/gmx-smd/vis/test1
# 优化初始构象的选取，根据前期得到的PMF，高能构象可能分布在质心距离为15-25埃米的范围，所以初始构象中，在保证有完全结合和解离的情况下，最好集中在这个区间里
# 根据smd得到的 smd_pullf.xvg；smd_pullx.xvg；再结合本身轨迹

vsREMD

0913 这一套是有钠离子和氯离子的

# /home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-31/gmx-smd/select
# 先把整个轨迹看一下，定位到后面400ps包含了配体解离的过程，把这部分轨迹间隔两帧输出为一个个独立的gro
# 注意，需要选择system(0)
gmx_mpi trjconv -f ../../test1/smd.xtc -s ../../test1/smd.tpr -sep -n index.ndx -o traj.gro -b 600 -skip 2 -pbc mol -ur compact
gmx_mpi trjconv -f ../../../test1/smd.xtc -s ../../../test1/smd.tpr -n index.ndx -o ref.gro -dump 0 -pbc mol -ur compact


# /home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-31/gmx-smd/select/test1
# pymol 里看这一段，挑选20个构象
for i in {0,4,18,21,28,32,34,36,38,43,54};do cp sep/traj${i}.gro ./;done
for i in {45..52};do cp sep/traj${i}.gro ./;done


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


# 到24 上跑 10个副本的 vsREMD 100个核
# 85.24 /home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-31/vsremd/gmx-smd
# md.mdp 和 temp 和之前的测试体系一致，Index为保险，重新生成。gmx.top检查了一致的
gmx_mpi make_ndx -f smd-select/smd-select0.gro -o index.ndx


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'


# 生成20个副本各自的文件
# 跑 10个副本的 vsREMD 100个核
# /home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-31/vsremd/gmx-smd/10replica-50ns
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 -replex 1000 &


# 85.23 上交一个 5 副本的测试
# /home/databank_70t/zzy/project/koff/trypsin-31/vsremd/gmx-smd/5replica-50ns
scp -r lywu@172.21.85.24:/home/databank/zzy/project/MD/koff/trypsin/complex-pre/trypsin-31/vsremd/gmx-smd/10replica-50ns/* ./
source ~/.gmx2022.5-vsremd.sh
nohup mpirun -np 100 gmx_mpi mdrun -v -deffnm md -multidir md0 md1 md2 md3 md4 -replex 1000 &

PMF

# 85.23 /home/databank_70t/zzy/project/koff/trypsin-31/vsremd/gmx-smd/5replica-50ns/analysis
scp lywu@172.21.85.24:/home/databank/lywu/trypsin/inhibitiors/trypsin-18/20replica/zzy-atom_center/pmf-2024-09-06-atomcenter.ipynb ./


gmx_mpi trjconv -f ../md1/smd.xtc -s ../md1/smd.tpr -n index.ndx -o dry.xtc -pbc mol -ur compact

# 选10个副本-温度前10个，初始构象需要放解离的

参考资料：

Steered Molecular Dynamics (uiuc.edu)