转好文章(amber 动力学)

amber 动力学过程：

1、选择力场 tleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99

2、导入晶体结构 model=loadpdb "sbp_lin.pdb"

保存crd和top文件 saveamberparm model polyAT_vac.top polyAT_vac.crd

此时注意电荷是否平衡：

如果缺正电荷 addions model Na+ 0 负离子就用Cl-

选择水箱 solvateoct model TIP3PBOX 8.0

3、保存crd和top文件 saveamberparm model model_wat.top model_wat.crd

4、退出tleap quit

5、保存新的pdb ambpdb -p model_wat.top <> model2.pdb

6、溶剂环境能量最优化。这一步保持溶质（蛋白）不变，去除溶剂中能量不正常的范德华相互作用。

该步骤的配置文件min1.in如下：

---------------------------------------------------------------------

oxytocin: initial minimisation solvent + ions ##说明信息######

&cntrl ##模拟参数起始

imin = 1, ##任务是优化，0 是分子动力学

cut = 10 ##非键相互作用的截断值为10 挨

ntb = 1, ##周期边界条件 0 不采用；1 定容 ；2 定压

maxcyc = 4000, ##优化步数

ntr = 1, ##优化时需要一些约束原子 -ref

ncyc = 2000, ##前2000最陡下降，后面步骤共轭梯度

/

Hold the protein fixed ##约束说明

500.0 ##作用在肽键上的力 kcal/mol

RES 1 9 ##限制的残基序号 同restrain=’:1-9’

END

END

------------------------------------------------------------------------------

任务命令：如果

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

7、对蛋白进行优化，min2.in文件将min1.in中的限制原子修改，限制水的位置。

也可以考虑利用restrainmask=’:1-9@CA,N,C’约束蛋白主链上的原子。

sander -O -i min2.in -p model_wat.top -c min1.rst -o min2.out -r min2.rst&

8、整体的优化,去掉限制条件

sander -O -i min3.in -p model_wat.top -c min2.rst -o min3.out -r min3.rst &

9、有限制的分子动力学

第一步分子动力学保持蛋白分子位置不变，但是不是完全固定每个原子，同时缓解蛋白分子周围的水分子，是溶剂环境能量优化。在这个步骤中，我们将主要目的是对特定的原子使用作用力使其能量优化。

Eq1.in 如下：

fix protein ,relax H2O

&cntrl

nstlim=25000, dt=0.002, ntx=1, irest=0, ntpr=500, ntwr=500,ntwx=500,

tempi=0.0,temp0=300,ntt=3, gamma_ln=1.0,

ntb=1, ntp=0,

nrespa=1,

cut = 10,

ntc=2,ntf=2,

NTR=1,

/

fix protein and HEM

10

RES 1 284

END

END

-----------------------------

nstlim = #：＃表示计算的步数。dt = 0.002：表示步长，单位为ps，0.002表示2fs。ntx=1 irest=0 默认 ntb = 1：表示分子动力学过程保持体积固定。

imin = 0：表示模拟过程为分子动力学，不是能量最优化。

temp0 = 300：表示最后系统到达并保持的温度，单位为K。

tempi = 100：系统开始时的温度。 ntc=2,ntf=2 忽略氢键

gamma_ln = 1：表示当ntt＝3时的碰撞频率，单位为ps-1（请参考AMBER手册）

ntt = 3：温度转变控制，3表示使用兰格氏动力学。

sander -O -i eq1.in -p model_wat.top -c min3.rst -o eq1.out -r eq1.rst -ref min3.rst -x eq1.mdcrd

11整系统分子动力学模拟: eq2

-------------------------------------------------

f2:500ps MD

&cntrl

imin = 0, irest=1, ntx=5,

ntb=2, pres0 = 1.0, ntp=1,

taup = 2.0, ntc=2, ntf=2,

cut = 10, ntr = 0,

ntt = 3, gamma_ln = 1.0,

tempi = 300.0 , temp0 = 300.0

nstlim=500000, dt=0.002,

ntpr=500, ntwr=500, ntwx=500

/

----------------------------------------------------------

ntb=2：表示分子动力学过程的压力常数。Pres0＝1：引用1个单位的压强。

ntp＝1：表示系统动力学过程各向同性。taup = 2.0：压力缓解时间，单位为ps。

使用以下命令进行MD：

sander -O -i eq2.in -p model_wat.top -c eq1.rst -o eq2.out -r eq2.rst -x eq2.mdcrd -ref eq1.rst &

结果处理：

在动力学过程中，将会产生 .rst 和.mdrcd 文件（需要的），由于crd文件很多，可以将其压缩成.gz 文件：gzip filename，将产生一个filename.gz 文件

做出已跑动力学的rmd图，判断是否平衡：

ptraj xxx.top <>

xxx.in 内容：

trajin eq2.mdcrd.gz

trajin eq3.mdcrd.gz

trajin eq4.mdcrd.gz

trajin eq5.mdcrd.gz

rms first mass out xin.rms.dat1 :1-284@CA,C,N time 0.1

#产生一个xin.rms.dat1的文件，整体1-284骨架上的C与N原子

rms first mass out xin.rms.dat2 :88-91,172,201-204,230@CA,C,N time 0.1

#产生一个xin.rms.dat2的文件，保守残基骨架上的C与N原子

----------------------------------------------------------------------------

利用xmgrace 作图：

xmgrace xin.rms.dat1 xin.rms.dat2

如果需要取随即的点的构型：

ptraj xxx.top <>

xxx.in内容：

trajin eq7.mdcrd.gz 117 117 #eq7中的117ps 的构型

strip :WAT #冒号前有空格，后没有，注意wat与top的水的大小写一致

strip :Na+

trajout eq7.pdb pdb nobox 产生一个eq7.pdb.117的文件

-------------------------------------------

traj AR.top <<>

> trajin AR.md7.crd 300 500

> center :1-250

> image center familiar

> rms first mass out rms.dat :1-250

> average average.rst restrt

> average average.pdb pdb

> EOF

对于有些小分子利用sybyl构建的一般缺少原子参数，首先将小分子单独存成mol2，用amber的antechamber 做出参数。

1、利用gaussian

* antechamber -i 49.mol2 -fi mol2 -o 49.in -fo gzmat

这样可以生成49.in文件，下载到windows环境，运行gaussian计算这个文件，如果发现计算时间过长或者内存不足计算中断，可以修改文件选择小一些的基组。

you have Gassian output file.out (don’t forget to put some keyword in the Gassian input for the special format used by amber)

* -antechamber –i file.out –fi gout –o filr.prep –fo prepi –c resp (this step prepare the “file.prep” for your drug -c bcc all right )

* -parmchk –i file.prep –f prepi –o file.frcmod )this step prepare the file “file.frcmod”for the drug)

2、 直接用mol2 转

也可以直接把sybyl 的mol2 做参数

* -antechamber –i file.mol2 –fi mol2 –o filr.prep –fo prepi –c resp / -c bcc

* -parmchk –i file.prep –f prepi –o file.frcmod

now you have 2 files: file.prep and file.frcmod for the drug .

注意修改file.prep 中文件名。要与pdb中的配体一致。（vi prep。假设是abc）

对应pdb与prep文件中的原子名称,，原子顺序要一致。

prep 在xleap中打开

using tleap for the protein or drug or/and complex

-xleap –s –f leaprc.ff99

-mod = loadamberparams file.frcmod

-loadamberprep file.prep

edit abc

pdb 可以用sybyl打开。用vi 修改pdb

-source leaprc.gaff

-RL = loadpdb file.pdb

-Solvateoct RL TIP3PBOX 10

-Addions RL Cl- 0 (0:zero for neutral charge,Cl- can be replaced by Na+)

-Saveamberparm RL fileWT.top fileWT.crd (save topology and coord. With water)