3cl共价化合物MD

先将复合物结构去氢保存成pdb格式，然后文本大开将化合物的坐标信息调整到与之相联的氨基酸后，将氨基酸名称、所在的链均修改一致，原子类型需与bcc点和文件中的一致：

然后对蛋白质进行质子化预测，主要是确定HIS的构型，根据out文件修改上述pdb文件的HIS：

#激活pdb2pqr:
conda activate pdb2pqr预测
#质子化
pdb2pqr30 wt.pdb wt.pqr --ff AMBER --ffout AMBER --with-ph 7.4 --pdb-output wt_out.pdb

配体处理：共价配体当做非标准氨基酸处理

将上述文件中的配体与形成共价键的氨基酸取出，对N端加H，C端加OH，用pymol加H后存成pdb格式，然后用高斯进行优化，取的resp电荷：

#gjf文件：
%chk=molecule
%mem=4GB
%nproc=8
#B3LYP/6-31G* scrf=(solvent=water) SCF=tight Test Pop=MK 
iop(6/33=2) iop(6/42=6) opt


nohup g16 4wi.gjf&

从log文件中提取resp电荷，得到ac文件：

antechamber -i 4wi.log -fi gout -o 4wi.ac -fo ac -c resp -at amber
#注意是否有不识别的原子类型DU

还要用上述的pdb文件生成bcc电荷，方便对应原子序号，蛋白质文件、ac文件的原子命名和原子类型需一一对应：

antechamber -fi pdb -i 4wi.pdb -bk 4WI -fo ac -o 4wi-bcc.ac -c bcc -at amber

创建4wi.mc文件：

#文件内容：
HEAD_NAME N
TAIL_NAME C
MAIN_CHAIN CA
OMIT_NAME H1
OMIT_NAME H2
OMIT_NAME H01
OMIT_NAME O01
PRE_HEAD_TYPE C
POST_TAIL_TYPE N
CHARGE 0.0
#HEAD_NAME和TAIL_NAME定义这个残基的头和尾原子（和其他标准残基形成肽键的原子）
MAIN_CHAIN定义氨基酸残基中连接头原子和尾原子的原子名
OMIT_NAME定义形成肽键后，会被删掉的O、H、N
PRE_HEAD_TYPE和POST_TAIL_TYPE氨基酸残基连接的原子类型（这会影响到后面的力场参数）
————————————————

用prepgen生成prepin文件：

prepgen -i 4wi.ac -o 4wi.prepin -m 4wi.mc -rn 4WI

为非标准残基生成力场文件：

parmchk2 -i 4wi.prepin -f prepi -o frcmod.cro -a Y -p $AMBERHOME/dat/leap/parm/parm10.dat

将带有“ATTN”标签的参数去除，用默认的gaff生成参数文件，补上parm10的缺陷:

grep -v "ATTN" frcmod.cro > frcmod1.cro # Strip out ATTN lines
parmchk2 -i 4wi.prepin -f prepi -o frcmod2.cro

tleap生成top文件

source leaprc.protein.ff14SB
source leaprc.gaff
source leaprc.water.tip3p
set default PBRadii mbondi3
#locate一下目录
loadAmberPrep 4wi.prepin
loadAmberParams frcmod2.cro
loadAmberParams frcmod1.cro
com = loadPDB 7rfw_dry.pdb  #所有connect信息都删掉，去氢
saveamberparm com native.prmtop native.inpcrd
#溶剂化，加盒子
solvatebox com TIP3PBOX 10
#检查电荷
charge com
#加载定义蛋白质和配体
addionsrand com Na+ 50 Cl- 46 
#再次检查电荷为0：
charge com
#保存蛋白质和配体复合物的拓扑文件和坐标文件
saveamberparm com complex.prmtop complex.inpcrd
#保存蛋白质和配体复合物的pdb文件
savepdb com com.pdb

用amber预处理体系

#能量最小化
pmemd.cuda -O -i min.in -p complex.prmtop -c complex.inpcrd -o min.out -r min.rst -ref complex.inpcrd
min.in中的参数，重点关注：
cut=10.0  （cut指cutoff，一般用10.0埃米，所有.in参数文件中尽量保持一致）
restraintmask=':1-306'    (这里的1-439修改成你的体系复合物实际的残基数目编号)


# 升温 
pmemd.cuda -O -i heat.in -p complex.prmtop -c min.rst -o heat.out -r heat.rst -ref min.rst


# 预平衡 
pmemd.cuda -O -i density.in -p complex.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 complex.prmtop -c density.rst -o md.out -x md.crd -r md.rst -ref density.rst &
nohup pmemd.cuda -O -i md.in -p complex.prmtop -c md.rst -o md_2.out -x md_2.crd -r md_2.rst -ref md.rst &
cut=10.0， (与min.in, heat.in, density.in保持一致，建议统一用10.0)
后边采用了不加Zn离子的蛋白质，restraintmask=':ZN，167，170，218，220'和ntr那两行删掉才正唱运行

计算mmgbsa