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

相关文章推荐

想出国的勺子 · 行业大模型标准体系及能力架构研究报告（附下载 ...· 1 年前 ·

傻傻的课本 · 如何在正则表达式中对嵌套捕获组进行编号？-腾 ...· 2 年前 ·

不要命的芒果 · 将标题添加到pandas hist ...· 3 年前 ·

爽快的莲藕 · r - ggpubr: change ...· 3 年前 ·

豁达的匕首 · NSToolbar.SetConfigura ...· 3 年前 ·

之后将产生一个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：

建立两个输入文件：

polyAT_gb_md1_12Acut.in : 12.0 angstrom long range cutoff, Generalized Born

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
process_mdout.perl “../ polyAT_gb_md1_12Acut.out”
cd ../polyAT_gb_md1_nocut
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的轨迹：