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的軌跡：