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格式所有键都被链接起来(异常)