Amber学习第一天：模拟DNA的片段

Simulating a DNA polyA-polyT Decamer

1.用NAB创建DNA双螺旋并产生可用的文件

(1)创建一个文件nuc.nab,写入：

molecule m;
m = fd_helix( "abdna", "aaaaaaaaaa", "dna" );
putpdb( "nuc.pdb", m, "-wwpdb");

(2)运行这个文件：

nab nuc.nab

./a.out

之后将产生一个nuc.pdb文件，这个文件是双螺旋结构模型。

(3)将这个文件导入leap中：

在终端输入xleap或tleap,再输入力场source leaprc.ff99SB,即$AMBERHOME/exe/xleap -s -f $AMBERHOME/dat/leap/cmd/leaprc.ff99SB

导入pdb文件：> dna1=loadpdb "nuc.pdb" 在xleap窗口中输出：Loading PDB file: ./nuc.pdb  total atoms in file: 638  说明导入正确。

在xleap中看这个结构：> edit dna1

(4)产生.prmtop和.inpcrd文件：

保存：> saveamberparm dna1 ployAT_vac.prmtop ployAT_vac.inpcrd

会输出：

Checking Unit.
WARNING: The unperturbed charge of the unit: -18.000000 is not zero.
 -- ignoring the warning.
Building topology.
Building atom parameters.
Building bond parameters.
Building angle parameters.
Building proper torsion parameters.
Building improper torsion parameters.
 total 110 improper torsions applied
Building H-Bond parameters.
Not Marking per-residue atom chain types.
Marking per-residue atom chain types.
 (no restraints)

和创建两个文件:

polyAT_vac.prmtop：参数／拓扑文件。它是静态的，在模拟的过程中不会发生变化。

polyAT_vac.inpcrd：坐标，不是静态的在模拟中会发生变化。

(5)加中性化离子：

>addions dna1 Na+ 0  (0意味着中性化)

18 Na+ ions required to neutralize.
Adding 18 counter ions to "dna1" using 1A grid
Grid extends from solute vdw + 1.87  to  7.97
Resolution:      1.00 Angstrom.
grid build: 0 sec
 (no solvent present)
Calculating grid charges
charges: 0 sec
Placed Na+ in dna1 at (-0.73, 10.83, 18.51).
Placed Na+ in dna1 at (6.27, -8.17, 18.51).
Placed Na+ in dna1 at (10.27, 4.83, 11.51).
Placed Na+ in dna1 at (-6.73, -8.17, 11.51).
Placed Na+ in dna1 at (-5.73, 3.83, 13.51).
Placed Na+ in dna1 at (-9.73, 3.83, 25.51).
Placed Na+ in dna1 at (6.27, -8.17, 5.51).
Placed Na+ in dna1 at (-5.73, 8.83, 1.51).
Placed Na+ in dna1 at (11.27, 1.83, 25.51).
Placed Na+ in dna1 at (-10.73, -4.17, 28.51).
Placed Na+ in dna1 at (5.27, 3.83, 3.51).
Placed Na+ in dna1 at (2.27, -6.17, 27.51).
Placed Na+ in dna1 at (-10.73, -3.17, 8.51).
Placed Na+ in dna1 at (6.27, 7.83, 15.51).
Placed Na+ in dna1 at (10.27, -3.17, 8.51).
Placed Na+ in dna1 at (-1.73, -12.17, 16.51).
Placed Na+ in dna1 at (-9.73, 3.83, 4.51).
Placed Na+ in dna1 at (-7.73, 8.83, 21.51).
Done adding ions.
产生中性化的prmtop和inpcrd文件：
>saveamberparm dnal ployAT_cio prmtop ployAT_cio.inpcrd
输出：polyAT_cio.prmtop,  polyAT_cio.inpcrd
最终的目的是建立带有抗性离子的可溶的DNA文件，抗性已建立，接下来创建可溶性的文件，也就是加水。
先复制dnal : 
>dna2 = copy dna1
在DNA周围创建矩形的水盒子：
>solvatebox dna1 TIP3PBOX 8.0
>edit dna1
在DNA周围创建八面体的水盒子：
>solvateoct dna2 TIP3PBOX 8.0
>edit dna2
再保存：
saveamberparm dna2 polyAT_wat.prmtop polyAT_wat.inpcr
Output files:  polyAT_wat.prmtop,  polyAT_wat.inpcrd

到现在已经产生运行最小化和分子动力学的文件。下部分将要用到。

2.最小化和分子动力学的运行

(1)最小化  目的：为了出去坏的联系

sander的语法：sander [-O] -i mdin -o mdout -p prmtop -c inpcrd -r restrt [-ref refc] [-x mdcrd] [-v mdvel] [-e mden] [-inf mdinfo]

• Arguments in []'s are optional
• If an argument is not specified, the default name will be used.
• -O    overwrite all output files (the default behavior is to quit if any output files already exist)
• -i      the name of the input file (which describes the simulation options), mdin by default.
• -o     the name of the output file, mdout by default.
• -p     the parameter/topology file, prmtop by default.
• -c     the set of initial coordinates for this run, inpcrd by default.
• -r     the final set of coordinates from this MD or minimization run, restrt by default.
• -ref  reference coordinates for positional restraints, if this option is specified in the input file, refc by default.
• -x    the molecular dynamics trajectory file (if running MD), mdcrd by default.
• -v    the molecular dynamics velocities file (if running MD), mdvel by default.
• -e    a summary file of the energies (if running MD), mden by default.
• -inf  a summary file written every time energy information is printed in the output file for the current step of the minimization of MD, useful for checking on the progress of a simulation, mdinfo by default.

能量优化由sander模块完成，运行sander至少需要三个输入文件，分别是分子的拓扑文件，坐标文件以及sander的控制文件。现在分子的拓扑文件和坐标文件已经完成，需要建立mdin输入文件 : 

polyAT_vac_init_min.in中写入：

polyA-polyT 10-mer: initial minimization prior to MD
  &cntrl
  imin   = 1,
  maxcyc = 500,
  ncyc   = 250,
  ntb    = 0,
  igb    = 0,
  cut    = 12
 /

如：一些变量的意义：

Initial minimisation of our structures
&cntrl
imin=1, maxcyc=4000, ncyc=2000,
cut=10, ntb=1,ntr=1,
restraint_wt=0.5
restraintmask=':1-283'
/

文件首行是说明，说明这项任务的基本情况； &cntrl与/之间的部分是模拟的参数

其中imin=1表示任务是能量优化，maxcyc=4000表示能量优化共进行4000步，ncyc=2000表示在整个能量优化的4000步中，前2000步采用最陡下降法，在2000步之后转换为共轭梯度法，如果模拟的时候不希望进行方法的转换，可以再加入另一个关键词NTMIN，如果NTMIN=0则全程使用共轭梯度法，NTMIN=2则全程使用最陡下降法，此外还有=3和=4的选项，分别是xmin法和lmod法，具体情况可以看手册。

第二行的cut=10表示非键相互作用的截断值，单位是埃， ntb=1表示使用周期边界条件，这个选项要和前面生成的拓扑文件坐标文件相匹配，如果前面加溶剂时候用的是盒子水，就设置ntb=1，如果加的是层水，那就应该选择ntb=0；ntr=1表示在能量优化的过程中要约束一些原子，约束的是哪些原子呢？后面有。

第三行和第四行都是关于约束原子的信息，restraint_wt=0.5限定了约束的力常数，在这里约束原子就是把原子用一根弹簧拉在固定的位置上，一旦原子偏离固定的位置，系统就会给他施加一个回复力，偏离的越远，回复力越大，回复力就是由这个力常数决定的，单位是Kcal/(mol*A)。restraintmask=':1-283'表示约束的是从1到283号残基，在这个分子中，1-283号残基是蛋白质上的氨基酸残基，从284号开始，就都是水了，所以这个命令的意思就是，约束整个蛋白质，自由地优化溶剂分子。因为溶剂分子是前面的tleap自动给加上的，所以一定要额外多关注一些。

(注：imin = 0 没有最小化能量优化（仅仅做分子动力学，默认）imin＝1 进行最小化能量优化（不做分子动力学）imin＝5 读轨迹并分析)

用sander运行最小化：

在工作目录下终端输入：sander -O -i polyAT_vac_init_min.in -o polyAT_vac_init_min.out -c polyAT_vac.inpcrd -p polyAT_vac.prmtop -r polyAT_vac_init_min.rst

输入文件：polyAT_vac_init_min.in, polyAT_vac.inpcrd, polyAT_vac.prmtop
输出文件:  polyAT_vac_init_min.out, polyAT_vac_init_min.rst

创建最小化后的PDB文件：

在终端输入：ambpdb -p polyAT_vac.prmtop < polyAT_vac_init_min.rst > polyAT_vac_init_min.pdb

之后就能在xleap看这个分子的结构。

(2)在真空(气相中)运行分子动力学

我们将运行来两个：一个截断值为12埃，一个为0

创建polyAT_vac_md1_12Acut.in文件写入：   

10-mer DNA MD in-vacuo, 12 angstrom cut off         
 &cntrl
  imin = 0, ntb = 0,                                                              //imin = 0,关闭最小化 ；ntb = 0,没有周期
  igb = 0, ntpr = 100, ntwx = 100,                                    //igb = 0,没有溶剂存在；后面的每一百步输出一个轨迹坐标
  ntt = 3, gamma_ln = 1.0,                                               // ntt = 3,恒温器(Langevin dynamics);
  tempi = 300.0, temp0 = 300.0                                        //温度为300K
  nstlim = 100000, dt = 0.001,                                        //运行100000步，每步为0.001秒，所以总长为100000＊0.001＝100
  cut = 12.0                                                                         //截断值为12埃
 /

创建polyAT_vac_md1_nocut.in文件写入：To run without a cutoff we simply set CUT to be larger than the extent of the system (e.g. CUT = 999).

10-mer DNA MD in-vacuo, infinite cut off
 &cntrl
  imin = 0, ntb = 0,
  igb = 0, ntpr = 100, ntwx = 100,
  ntt = 3, gamma_ln = 1.0,
  tempi = 300.0, temp0 = 300.0
  nstlim = 100000, dt = 0.001,
  cut = 999
 /

 在工作目录下终端运行：sander -O -i polyAT_vac_md1_12Acut.in -o polyAT_vac_md1_12Acut.out -c polyAT_vac_init_min.rst -p polyAT_vac.prmtop -r polyAT_vac_md1_12Acut.rst -x polyAT_vac_md1_12Acut.mdcrd

                                                sander -O -i polyAT_vac_md1_nocut.in -o polyAT_vac_md1_nocut.out -c polyAT_vac_init_min.rst -p polyAT_vac.prmtop -r polyAT_vac_md1_nocut.rst -x polyAT_vac_md1_nocut.mdcrd

看运行的进程：输入：tail -f polyAT_vac_md1_12Acut.out

Input files: polyAT_vac_md1_12Acut.in, polyAT_vac_md1_nocut.in, polyAT_vac_init_min.rst(最小化的文件), polyAT_vac.prmtop(开始没有加水的文件)

Output files: polyAT_vac_md1_12Acut.out, polyAT_vac_md1_nocut.out, polyAT_vac_md1_12Acut.rst, polyAT_vac_md1_nocut.rst, polyAT_vac_md1_12Acut.mdcrd, polyAT_vac_md1_nocut.mdcrd

(3)分析结果

a.提取能量：截断值12埃的文件:

mkdir polyAT_vac_md1_12Acut
cd polyAT_vac_md1_12Acut

perl process_mdout.perl “../polyAT_vac_md1_12Acut.out”    ／／运行perl程序会自动生成，没有截断值也一样运行

xmgr ./polyAT_vac_md1_12Acut/summary.EPTOT ./polyAT_vac_md1_nocut/summary.EPTOT

如果在redhat上没有安装xmgr的话，首先在http://rpmfind.net/linux/rpm2html/search.php?query=libXp.so.6下载Red Hat Linux 5.0 for i386   的XFree86-libs-3.3.1-14.i386.rpm包   安装：rpm -ivh XFree86 ＊.rpm

                                                         之后在ftp://plasma-gate.weizmann.ac.il/pub/xmgr4/RPMS/i386/README.html下载 Binary RPM for libc.so.6 的xmgr-semistatic-4.1.2-12.i386.rpm.hurricane 安装：rpm -ivh xmgr-se * .rpm.hurricane

这就OK 了，能运行 xmgr ./polyAT_vac_md1_12Acut/summary.EPTOT ./polyAT_vac_md1_nocut/summary.EPTOT ，就可以看两个文件总电荷的图了，从而可以区别了。

b.计算Rmsd值随时间的变化：

首先建两个文件：polyAT_vac_md1_12Acut.calc_rms 和 polyAT_vac_md1_nocut.calc_rms

分别写入：trajin polyAT_vac_md1_12Acut.mdcrd
                 rms first mass out polyAT_vac_md1_12Acut.rms time 0.1

                 trajin polyAT_vac_md1_12Acut.mdcrd
                 rms first mass out polyAT_vac_md1_nocut.rms time 0.1

运行：ptraj polyAT_vac.prmtop < polyAT_vac_md1_12Acut.calc_rms

           ptraj polyAT_vac.prmtop < polyAT_vac_md1_nocut.calc_rms

会产生两个文件：polyAT_vac_md1_12Acut.rms ，polyAT_vac_md1_nocut.rms

之后运行：xmgr polyAT_vac_md1_12Acut.rms polyAT_vac_md1_nocut.rms

横坐标为ps,纵坐标为埃值。12埃截断值的图为一条直线，很平稳；而没有截断值的为逐步上升的趋势。

c.用VMD来看轨迹

安装VMD；下载vmd-1.8.7.bin.LINUX.opengl.tar.gz
解压：tar -zxvf vmd-1.8.7.bin.LINUX.opengl.tar.gz

          cd vmd-1.8.7

           ./configure  LINUX OPENGL FLTK TK ACTC CUDA IMD LIBSBALL XINPUT LIBTACHYON VRPN NETCDF TCL PTHREADS SILENT (在文件configure.options)

         cd src/

        make install

安装好了：启动文件在/usr/local/bin下的vmd

所以设置环境变量：gedit /etc/profile  中写入：        

 #path for VMD
export PATH_HOME="/opt/vmd-1.8.7/"
export CLASSPATH="/usr/local/bin"
pathmunge /usr/local/bin

之后在终端输入vmd即可启动。

(4)在隐式的溶液中运行最小化和分子动力学

首先最小化初始的系统：用到先前开始的文件, polyAT_vac.prmtop, polyAT_vac.inpcrd ，并建立polyAT_gb_init_min.in 写入:

polyA-polyT 10-mer: initial minimization prior to MD GB model
 &cntrl
  imin   = 1,
  maxcyc = 500,
  ncyc   = 250,
  ntb    = 0,
  igb    = 1,
  cut    = 12
 /

运行：sander -O -i polyAT_gb_init_min.in -o polyAT_gb_init_min.out -c polyAT_vac.inpcrd -p polyAT_vac.prmtop -r polyAT_gb_init_min.rst

Input files: polyAT_gb_init_min.in, polyAT_vac.inpcrd, polyAT_vac.prmtop
Output files: polyAT_gb_init_min.out, polyAT_gb_init_min.rst

产生PDB文件：ambpdb -p polyAT_vac.prmtop < polyAT_gb_init_min.rst > polyAT_gb_init_min.pdb

在隐式溶液中进行MD：

建立两个输入文件：

1. polyAT_gb_md1_12Acut.in: 12.0 angstrom long range cutoff, Generalized Born
2. polyAT_gb_md1_nocut.in:no long range cutoff, Generalized Born

分别写入：

10-mer DNA MD Generalized Born, 12 angstrom cut off
 &cntrl
  imin = 0, ntb = 0,
  igb = 1, ntpr = 100, ntwx = 100,
  ntt = 3, gamma_ln = 1.0,
  tempi = 300.0, temp0 = 300.0
  nstlim = 100000, dt = 0.001,
  cut = 12.0
 /

10-mer DNA MD Generalized Born, infinite cut off
 &cntrl
  imin = 0, ntb = 0,
  igb = 1, ntpr = 100, ntwx = 100,
  ntt = 3, gamma_ln = 1.0,
  tempi = 300.0, temp0 = 300.0
  nstlim = 100000, dt = 0.001,
  cut = 999

用最小化的文件运行MD：
有截断值：sander -O -i polyAT_gb_md1_12Acut.in -o polyAT_gb_md1_12Acut.out -c polyAT_gb_init_min.rst -p polyAT_vac.prmtop -r polyAT_gb_md1_12Acut.rst -x polyAT_gb_md1_12Acut.mdcrd

没有截断值：sander -O -i polyAT_gb_md1_nocut.in -o polyAT_gb_md1_nocut.out -c polyAT_gb_init_min.rst -p polyAT_vac.prmtop -r polyAT_gb_md1_nocut.rst -x polyAT_gb_md1_nocut.mdcrd 

Input files: polyAT_gb_md1_12Acut.in, polyAT_gb_md1_nocut.in, polyAT_gb_init_min.rst, polyAT_vac.prmtop

Output files: polyAT_gb_md1_12Acut.out, polyAT_gb_md1_nocut.out, polyAT_gb_md1_12Acut.rst, polyAT_gb_md1_nocut.rst, polyAT_gb_md1_12Acut.mdcrd, polyAT_gb_md1_nocut.mdcrd

结果分析：

能量的比较：

mkdir polyAT_gb_md1_12Acut
mkdir polyAT_gb_md1_nocut
cd polyAT_gb_md1_12Acut
perl process_mdout.perl “../polyAT_gb_md1_12Acut.out”
cd ../polyAT_gb_md1_nocut
perl process_mdout.perl “../polyAT_gb_md1_nocut.out”
cd ..
xmgrace ./polyAT_gb_md1_12Acut/summary.EPTOT ./polyAT_gb_md1_nocut/summary.EPTOT

RMSD的比较：

建立俩个ptraj（polyAT_gb_md1_12Acut.calc_rms, polyAT_gb_md1_nocut.calc_rms.）的文件分别写入：

FILE1
trajin polyAT_gb_md1_12Acut.mdcrd
rms first mass out polyAT_gb_md1_12Acut.rms time 0.1

FILE2
trajin polyAT_gb_md1_nocut.mdcrd
rms first mass out polyAT_gb_md1_nocut.rms time 0.1

运行：

       ptraj  polyAT_vac.prmtop  <  polyAT_gb_md1_12Acut.calc_rms

      ptraj  polyAT_vac.prmtop  <  polyAT_gb_md1_nocut.calc_rms

然后拟合RMSD值：

        xmgrace polyAT_gb_md1_12Acut.rms  polyAT_gb_md1_nocut.rms

最后用VMD看这两个DNA的轨迹：