cMD(普通MD）-amber

蛋白质处理

需要先对蛋白质进行质子化预测：

#激活pdb2pqr:  （gp4n)
conda activate pdb2pqr预测
#质子化
pdb2pqr30 ${protein}.pdb ${protein}.pqr --ff AMBER --ffout AMBER --with-ph 7.4 --pdb-output ${protein}_out.pdb


#Tleap处理蛋白质：准备好空白蛋白质.pdb
#进入tleap：
tleap
#加载力场：
source oldff/leaprc.ff14SB
source leaprc.gaff
#locate一下目录
loadAmberParams frcmod.ionsjc_tip3p
#加载定义蛋白质
pro=loadpdb ${protein}_out.pdb     
#load蛋白和配体的结构之后，有二硫键的需要添加二硫键，无则不需要
bond pro.22.SG pro.96.SG
（这条命令是指pro中第残基22的SG原子与pro中残基96的SG原子形成共价键，注意每一个二硫键都要单独用一条命令进行bon

#溶剂化，加盒子，默认10

solvatebox pro TIP3PBOX 10
#检查电荷
charge pro
#charge com 是检查复合物体系带电情况，用于后面添加电荷，使体系为电中性，计算电荷结果：氯离子是64个，为了保证中和体系，添加70个Na+,64个Cl-
#（添加正负电荷的数目需要根据盒子的体积进行计算，使得最后的离子浓度为0.15M，同时也要保证体系为电中性）
计算所需要的Cl电荷：1.5 Calculating Salt Molarity in an Explicit Water System (ambermd.org)
#加正负电荷
addionsrand pro Na+ 70 Cl- 64 
#再次检查电荷为0：
charge pro
#保存蛋白质受体的拓扑文件和坐标文件
saveamberparm pro anti.prmtop anti.inpcrd 
#保存蛋白质加完盒子的pdb文件
savepdb pro pro-${protein}.pdb
ctrl+C退出tleap
一般使用tleap处理之后的结构，会自动对残基进行重新编号，这里导出pdb是为了检查重新编号的情况，后面的步骤会用到。
模拟小组公共空间上传了后面用于能量最小化等步骤的四个.in参数文件，链接如下：AMBER-MD参数文件

#能量最小化

pmemd.cuda -O -i min.in -p anti.prmtop -c anti.inpcrd -o min.out -r min.rst -ref anti.inpcrd
min.in中的参数，重点关注：
cut=8.0  （cut指cutoff，一般用10.0埃米，所有.in参数文件中尽量保持一致）
restraintmask=':1-439'    (这里的1-439修改成你的体系复合物实际的残基数目编号)
格式后查看：
cpptraj -p anti.prmtop -y min.rst -x min.rst7
ambpdb -p anti.prmtop < min.rst7 > min.pd

# 升温 

pmemd.cuda -O -i heat.in -p anti.prmtop -c min.rst -o heat.out -r heat.rst -ref min.rst
是内存被占满了，换了208，但还是报错，怀疑从能量最小化开始，可能就有问题，将min.rst文件转成PDB，发现还是体系有问题，蛋白质有断裂，体系崩溃了，与老师交流后，选择不加Zn离子，因为它里口袋很远，也不是发挥功能所必须的，而且之前彭师兄跑MD额时候就没有加Zn，那4个CYS都跑开了


注意检查和修改的部分同min.in

# 预平衡 

pmemd.cuda -O -i density.in -p anti.prmtop -c heat.rst -o density.out -r density.rst -ref heat.rst
注意检查和修改的部分同min.in

# 提交MD，注意调整md.in文件里的nstlim

nohup pmemd.cuda -O -i md.in -p anti.prmtop -c density.rst -o md.out -x md.crd -r md.rst -ref density.rst &
nohup pmemd.cuda -O -i md.in -p anti.prmtop -c md.rst -o md_2.out -x md_2.crd -r md_2.rst -ref md.rst &
restraintmask=':ZN，167，170，218，220' (这部分用于位置限制的残基修改为ZN2+以及与其形成配位键的残基)
nstlim=100000000， （步长，该参数决定要模拟的时长，参考手册进行调整）
cut=10.0， (与min.in, heat.in, density.in保持一致，建议统一用10.0)
后边采用了不加Zn离子的蛋白质，restraintmask=':ZN，167，170，218，220'和ntr那两行删掉才正唱运行