Amber分子动力学简例1NPO

Amber分子动力学简例1NPO

#引用+原创

引用：

https://blog.csdn.net/rogerzhanglijie/article/details/8084187?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-1.control&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-1.control

PDB code： 1NPO.pdb
删除ACD，保留B链（有一对二硫键），命名为oxyt.pdb

1、激活amber环境:

conda antivate amber

2、删除的氢原子和水分子（其他软件添加的氢可能无法识别，特别S-S上可能 含有孤对电子，需要清除）

pdb4amber -i oxyt.pdb -o oxyy.pdb -y -d
pdb4amber -i oxyy.pdb -o oxy.pdb --reduce

(利用amber自带的reduce模块加氢，注：加氢同时使用 --add-missing-atoms会有异常）
pdb4amber处理oxy.pdb后，已将CYS改为CYX，并连接好二硫键

2、准备分子top文件和crd文件

(1)启动tleap: tleap

(2)加载力场：source oldff/leaprc.ff03

(3)读入蛋白文件：oxy = loadpdb oxy.pdb

(4)检查文件 ：check oxy(注意是否有警告或报错，最后应该是OK）

(5)添加水模型参数:source leaprc.water.tip3p

(6)添加溶剂solvateoct oxy TIP3PBOX 9.0 (oct是八面体水，box是立方体）

(7)检查电性：charge oxy（该演示例为0，若非电中性，添加抗衡离子）

(8)保存文件：saveamberparm oxy oxy.top oxy.crd

(9)退出tleap：quit

3、水环境下的分子动力学模拟

(1)溶剂环境能量最优化。这一步保持溶质不变，去除溶剂中能量不正常的范德华相互作用

(2)整系统的能量最优化，去除整个系统中能量不正常相互作用。

(3)NVT模拟，温度逐渐从0K升到300K；

(4)NPT模拟，在1atm，300K环境下MD。

①溶剂环境能量最优化，配置文件min1.in
oxytocin:  initial minimisation solvent + ions
&cntrl
imin = 1,
maxcyc =  1000,
ncyc = 500,
ntb = 1,
ntr = 1,
cut = 10
/
Hold the  protein fixed
500.0   #单位kcal/mol，表示作用在肽键上使其不动的力
RES 1  9 #表示肽链残基数目，oxy.pdb只有9个残基
END
END

模拟命令如下：

sander –O –i min1.in  –o min1.out –p oxy.top –c oxy.crd –r oxy_min1.rst –ref oxy.crd  &

②整系统的能量最优化，配置文件min2.in
oxytocin:  initial minimisation whole system
&cntrl
imin = 1,
maxcyc =  2500,
ncyc = 1000,
ntb = 1,
ntr = 0,
cut =  10
/

模拟命令如下：

sander  –O –i min2.in –o min2.out –p oxy.top –c oxy_min1.rst –r oxy_min2.rst  &

③NVT模拟，配置文件md1.in
oxytocin:  20ps MD with res on protein
&cntrl
imin = 0, #表示模拟过程为分子动力学，不是能量最优化。
irest = 0,
ntx =  1,
ntb = 1,  #表示分子动力学过程保持体积固定
cut = 10,
ntr = 1,
ntc = 2,
ntf = 2,
tempi =  0.0, #系统开始时的温度。
temp0 = 300.0,#表示最后系统到达并保持的温度，单位为K。
ntt = 3, #温度转变控制，3表示使用兰格氏动力学
gamma_ln = 1.0,#表示当ntt＝3时的碰撞频率，单位为ps-1（请参考AMBER手册）
nstlim = 10000, ＃表示计算的步数。
dt =  0.002,#表示步长，单位为ps，0.002表示2fs。
ntpr = 100,
 ntwx = 500, 
ntwr = 1000
/
Keep protein fixed with  weak restraints
10.0
RES 1  9
END
END

额外参数： tautp = 0.1：热浴时间常数，缺省为1.0。小的时间常数可以得到较好的耦联。 vlimit = 20.0：保持分子动力学稳定性速度极限。20.0为缺省值，当动力学模拟中原子速度大于极限值时，程序将其速度降低到极限值以下。 comp = 44.6： 溶剂可压缩单位。 ntc = 2：Shake算法使用标志位。1表示不实用使用，2表示氢键将被计算，3表示所有键都将被计算在内。 tol = #.#####：坐标位置重新设置的几何位置相对容忍度。 我们将使用一个较小的作用力，10kcal/mol。在分子动力学中，当ntr＝1时，作用力只需要5-10kcal/mol（我们需要引用一个坐标文件做分子动力学过程的比较，我们需要使用"-ref"参数）。太大的作用力同时使用Shake算法和2fs步长将使整个系统变得不稳定，因为大的作用力使系统中的原子产生大频率的振动，模拟过程并步需要。

运行命令如下：

sander  –O –i md1.in –o md1.out –p oxy.top –c oxy_min2.rst –r oxy_md1.rst –x  oxy_md1.mdcrd –ref oxy_min2.rst –inf  md1.info&

④NPT模拟，配置文件md1.in
oxytocin:  250ps MD
&cntrl
imin = 0, irest = 1, ntx = 7,
ntb = 2, pres0 = 1.0,  ntp = 1,
taup = 2.0,
cut = 10, ntr = 0,
ntc = 2, ntf = 2,
tempi =  300.0, temp0 = 300.0,
ntt = 3, gamma_ln = 1.0,
nstlim = 125000, dt =  0.002,
ntpr = 100, ntwx = 500, ntwr =  1000
/

上面一些参数解释如下:
ntb=2：表示分子动力学过程的压力常数。
ntp＝1：表示系统动力学过程各向同性。
taup =2.0：压力缓解时间，单位为ps。
pres＝1：引用1个单位的压强。
使用以下命令进行MD：

sander –O –i md2.in  –o md2.out –p oxy.top –c oxy_md1.rst – r oxy_md2.rst –x oxy_md2.mdcrd –ref  oxy_md1.rst –inf md2.info  &

4、数据分析：process_mdout.perl md1.out md2.out

(1)可用xmgrace查看，如温度，压力，密度等

5、转换为pdb格式查看MD后构型：

(1)ambpdb -p oxy.top -c min1.rst > min1.pdb

(2)ambpdb -p oxy.top -c min2.rst > min2.pdb

(3)ambpdb -p oxy.top -c md1.rst > md1.pdb

(4)ambpdb -p oxy.top -c md2.rst > md2.pdb

注：ambpdb转换为pdb格式，S-S没有连接信息，需手动添加，转换为mol2格式所有键都被链接起来（异常）