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)