matlab 重映射,AmberMD仿真实战初步教程(根据俺实际仿真step by step记录)

最新推荐文章于 2023-10-06 16:16:36 发布

浪客宏非

最新推荐文章于 2023-10-06 16:16:36 发布

阅读量905

点赞数

文章标签： matlab 重映射

该教程详细介绍了如何使用AmberMD进行DNA模拟，从创建pdb结构，到力场参数化，再到最小化和分子动力学（MD）模拟。在真空中、隐式溶剂和显式溶剂环境下进行了MD模拟，并分析了能量、RMSD等关键参数，展示了在不同条件下的DNA行为。最后，通过VMD软件可视化了MD轨迹。

摘要由CSDN通过智能技术生成

说明，本实战教程是AmberMD的初步tutorial，俺根据官网，一步步实战，全部是好使的，才记录下来，与诸君共飨。

http://ambermd.org/tutorials/basic/tutorial1/index.htm

Amber学习：模拟DNA的片段 (有些翻译表达，不尽如人意，见谅。)

Part 1 Introduction

该学习手册首先演示了如何构建一个十倍体pdb结构，并利用这个结构构建sander运行必要的输入文件。MD仿真内核基于AMBER14。

MD仿真需要若干文件，主要包括：

 prmtop – 该文件包含分子的拓扑结构，和必要的力场参数。

 inpcrd (或一个上次运行产生的 restrt文件x.rst) –该文件包含原子坐标，或者也可能包含速度和周期性边界盒子尺寸等信息。periodic box dimensions.

 mdin – 该文件是sander的输入文件(*.in)，包含一系列文件名、用于决定仿真类型和仿真选项的控制变量。

该教程第一部分是，咱们先用AMBERTOOL创建prmtop和inpcrd文件，并分别置于真空系统和溶剂系统。然后执行最小化、MD仿真，最终得到仿真后的DNA。

由于使用显式溶剂进行这些模拟可能很昂贵，我们还将使用一些模型，其中隐含地包含溶剂效应。

Part 2 创建双链DNA polyA-polyT

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

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

molecule m;

m = fd_helix( "abdna", "aaaaaaaaaa", "dna" );

putpdb( "nuc.pdb", m, "-wwpdb&quot

;

(2)运行这个文件：

nab nuc.nab

./a.out

之后将产生一个nuc.pdb文件，这个文件是双螺旋结构模型。 &--如果已经有pdb文件，直接从第三步开始 --songjun

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

在终端输入tleap,

再输入力场 source leaprc.ff99SB,

在官网 http://ambermd.org/tutorials/basic/tutorial1/section2.htm

中有关于力场的说明，在早期版本(amber10 or 11)中推荐用ff99SB，但在amber12中推荐用ff12SB，在amber14中推荐ff14SB(它是ff99SB的优化版)

看amber版本用： echo $AMBERHOME

source leaprc.ff14SB

导入pdb文件：> dna1= loadpdb "nuc.pdb"

在tleap窗口中输出：Loading PDB file: ./nuc.pdb total atoms in file: 638 说明导入正确。

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

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

会输出：

Checking Unit.

WARNING: The unperturbed charge of the unit: -18.000000 is not zero.

-- ignoring the warning.

Building topology.

它将创建两个文件:

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

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

(5)加中性化离子：

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

将输出

addions dna1 Na+ 0

18 Na+ ions required to neutralize.

Adding 18 counter ions to "dna1" using 1A grid

Grid extends from solute vdw + 3.65 to 9.75

Resolution: 1.00 Angstrom.

grid build: 0 sec

(no solvent present)

Calculating grid charges

charges: 1 sec

Placed Na+ in dna1 at (6.44, 3.95, 17.79).

Placed Na+ in dna1 at (5.44, -5.05, 10.79).

Placed Na+ in dna1 at (-10.56, 5.95, 13.79).

Placed Na+ in dna1 at (-10.56, -6.05, 19.79).

Placed Na+ in dna1 at (-1.56, 11.95, 9.79).

Placed Na+ in dna1 at (-10.56, -4.05, 6.79).

Done adding ions.

产生中性化的prmtop和inpcrd文件：

> saveamberparm dna1 polyAT_cio.prmtop polyAT_cio.inpcrd

输出：polyAT_cio.prmtop, polyAT_cio.inpcrd

最终的目的是建立带有抗性离子的可溶的DNA文件，抗性已建立，接下来创建可溶性的文件，也就是加水。

> dna2 = copy dna1

> solvatebox dna1 TIP3PBOX 10.0

solvatebox dna1 TIP3PBOX 12.0 矩形盒子，尺寸12A

将输出

> solvatebox dna1 TIP3PBOX 8.0

Solute vdw bounding box: 27.738 26.738 40.099

Total bounding box for atom centers: 43.738 42.738 56.099

Solvent unit box: 18.774 18.774 18.774

Total vdw box size: 46.743 45.963 58.910 angstroms.

Volume: 126564.801 A^3

Total mass 56296.124 amu, Density 0.739 g/cc

Added 2767 residues.

> solvateoct dna2 TIP3PBOX 10.0 八面体盒子，尺寸10A

solvatebox dna2 TIP3PBOX 12.0

saveamberparm dna2 polyAT_wat.prmtop polyAT_wat.inpcrd

Output files: polyAT_wat.prmtop, polyAT_wat.inpcrd

savepdb DNA1 XXX.pdb &保存pdb的命令

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

&--------- 至此需要退出tleap, 用命令quit ！ ----------&

Part 3 最小化和分子动力学运行

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

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

 []中是可选项

 可选项为指定时，系统默认采用默认值。

 -O 覆盖原有的输出文件；

 -i 输入文件的名字, 默认是mdin文件(*.in)；

 -o 输出文件名,默认 mdout (*.out)；

 -p 参数/拓扑文件,默认 prmtop；

 -c 坐标文件,默认 inpcrd；

 -r 本次MD或最小化的坐标设置文件，默认restrt (*.rst)；

 -ref 位置约束的坐标参考文件, 是可选项,默认 refc (*.rst或inpcrd)；

 -x MD的轨迹文件 (针对MD仿真，若是最小化，无此项)，默认mdcrd；

 -v MD的速度文件，可选项(if running MD)，默认mdvel ；

 -e MD的能量总结文件，可选项 (if running MD)，默认 mden.

 -inf 最小化时每次写入能量信息的总结文件，在check仿真过程时很有用，默认 mdinfo文件。

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

polyAT_vac_init_min.in中写入：

polyA-polyT 10-mer: initial minimisation prior to MD

&cntrl

imin = 1,

maxcyc = 500,

ncyc = 250,

ntb = 0,

igb = 0,

cut = 12

接下来咱们第一次采用 sander命令来运行 for the first time

用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

Input files: polyAT_vac_init_min.in, polyAT_vac.inpcrd, polyAT_vac.prmtop

Output files: polyAT_vac_init_min.out, polyAT_vac_init_min.rst

这一步运行比较快，大约需要几秒的时间。

这一步生产的.rst文件是坐标文件，很关键。

此外，打开polyAT_vac_init_min.out能够发现能量在第一步和最后一步的变化：

NSTEP ENERGY RMS GMAX NAME NUMBER

1 -1.4123E+03 1.9444E+01 9.4276E+01 C2' 634

NSTEP ENERGY RMS GMAX NAME NUMBER

500 -2.3331E+03 6.1397E-01 3.8981E+00 C6 208

 利用坐标文件(.rst)创建 PDB文件

> ambpdb -p polyAT_vac.prmtop < polyAT_vac_init_min.rst > polyAT_vac_init_min.pdb

生成新的pdb文件polyAT_vac_init_min.pdb

 真空中运行 MD

我们将运行来两个：一个截断值为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秒

cut = 12.0 //截断值为12埃

空两行

创建polyAT_vac_md1_nocut.in文件，写入：

10-mer DNA MD in-vacuo, no 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命令，

第一个：

> 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

这个命令至少需要10分钟左右。

第二个：

> 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

这个也需要10分钟左右。

可用top命令看sander是否运行完毕。(上述命令结尾也可加符号 & )

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

输入文件包括: polyAT_vac_md1_12Acut.in, polyAT_vac_md1_nocut.in, polyAT_vac_init_min.rst, polyAT_vac.prmtop

输出文件包括: 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

 结果分析

(1) 从mdout文件中提取能量：截断值12埃的文件:

mkdir polyAT_vac_md1_12Acut

cd polyAT_vac_md1_12Acut

>process_mdout.perl ../polyAT_vac_md1_12Acut.out

同样处理nocut的文件

mkdir polyAT_vac_md1_nocut

cd polyAT_vac_md1_nocut

>process_mdout.perl ../polyAT_vac_md1_nocut.out

(2)接下来，画轨迹：

xmgrace ./polyAT_vac_md1_12Acut/summary.EPTOT ./polyAT_vac_md1_nocut/summary.EPTOT

xmgrace命令有问题，是因为没有安装相关xmgrace组件。它实际上是一款画曲线的软件，咱们可以copy出两个summary.EPTOT文件，用matlab来画曲线，也是一样。

注：summary.EPTOT中的数据摘录如下：

0.000 -2267.9419

0.100 -1570.4415

0.200 -1424.1312

0.300 -1513.0011

0.400 -1247.5700

0.500 -1101.3187

0.600 -1003.6084

0.700 -1183.8510

0.800 -1137.9767

0.900 -915.5758

1.000 -656.8817

1.100 -933.2193

1.200 -990.0922

1.300 -885.8356

1.400 -916.7535

问题解决，接着往下走。

(3)计算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_nocut.mdcrd

rms first mass out polyAT_vac_md1_nocut.rms time 0.1

运行 cpptraj 命令，如下:

$AMBERHOME/bin/cpptraj -p polyAT_vac.prmtop -i polyAT_vac_md1_12Acut.calc_rms

$AMBERHOME/bin/cpptraj -p polyAT_vac.prmtop -i polyAT_vac_md1_nocut.calc_rms

运行结果大致如下：

CPPTRAJ: Trajectory Analysis. V15.00

___ ___ ___ ___

| \/ | \/ | \/ |

_|_/\_|_/\_|_/\_|_

Reading 'polyAT_vac.prmtop' as Amber Topology

INPUT: Reading Input from file polyAT_vac_md1_12Acut.calc_rms

[trajin polyAT_vac_md1_12Acut.mdcrd]

Reading 'polyAT_vac_md1_12Acut.mdcrd' as Amber Trajectory

[rms first mass out polyAT_vac_md1_12Acut.rms time 0.1]

RMSD: (*), reference is first frame (*), with fitting, mass-weighted.

---------- RUN BEGIN -------------------------------------------------

PARAMETER FILES:

0: 'polyAT_vac.prmtop', 638 atoms, 20 res, box: None, 2 mol, 1000 frames

INPUT TRAJECTORIES:

0: 'polyAT_vac_md1_12Acut.mdcrd' is an AMBER trajectory, Parm polyAT_vac.prmtop (reading 1000 of 1000)

Coordinate processing will occur on 1000 frames.

TIME: Run Initialization took 0.0001 seconds.

BEGIN TRAJECTORY PROCESSING:

.....................................................

ACTION SETUP FOR PARM 'polyAT_vac.prmtop' (1 actions):

0: [rms first mass out polyAT_vac_md1_12Acut.rms time 0.1]

Target mask:

(638)

Reference mask:

(638)

----- polyAT_vac_md1_12Acut.mdcrd (1-1000, 1) -----

0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% Complete.

Read 1000 frames and processed 1000 frames.

TIME: Trajectory processing: 0.1576 s

TIME: Avg. throughput= 6345.0166 frames / second.

ACTION OUTPUT:

DATASETS:

1 data set:

RMSD_00000 "RMSD_00000" (double, rms), size is 1000

DATAFILES:

polyAT_vac_md1_12Acut.rms (Standard Data File): RMSD_00000

---------- RUN END ---------------------------------------------------

TIME: Total execution time: 0.1940 seconds.

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

To cite CPPTRAJ use:

Daniel R. Roe and Thomas E. Cheatham, III, "PTRAJ and CPPTRAJ: Software for

Processing and Analysis of Molecular Dynamics Trajectory Data". J. Chem.

Theory Comput., 2013, 9 (7), pp 3084-3095.

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

然后可以绘制rmsd随时间变化曲线了：

>xmgrace polyAT_vac_md1_12Acut.rms polyAT_vac_md1_nocut.rms

同样可以用matlab来画。官方网站的Tutorial是个好东西，错误比较少，很权威。

VMD使用方法-windows

用VMD可视化观测simualtion轨迹

http://ambermd.org/tutorials/basic/tutorial1/section3.htm

(1)下载VMD软件windows版本到自己的电脑。

https://www.ks.uiuc.edu/Research/vmd/

这是伊利诺伊大学香槟分校的产品。

(2)安装VMD,并启动 File -> New Molecule

首先装载polyAT_vac.prmtop 文件，

然后再装载轨迹文件 polyAT_vac_md1_12Acut.mdcrd ，该文件一般较大。

装载ok后，可以看到如下可视化界面：

(论坛中如何插入图片，俺还没学会，所以这里的图就不显示了。)

Part 4 隐式溶剂中最小化并运行MD

上述Section一直在真空中(Vacuo)模拟MD，接下来在隐式溶剂(implicit solvent)中运行MD

(1) 在MD模拟前，释放系统Relaxing the System Prior to MD

这一步，还使用原来的prmtop and inpcrd 文件，(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 命令，最小化

$AMBERHOME/bin/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

这一步大约需要30秒。

打开(polyAT_gb_init_min.out)文件可以发现，能量收敛

NSTEP ENERGY RMS GMAX NAME NUMBER

1 -2.2247E+03 1.9139E+01 9.3100E+01 C2' 634

NSTEP ENERGY RMS GMAX NAME NUMBER

500 -3.3168E+03 8.7962E-01 6.8734E+00 N1 618

生成新的PDB文件：

Again we can create PDB files for the start (polyAT_vac.inpcrd) and final structures (polyAT_gb_init_min.rst). [polyAT_gb_init_min.pdb].

$AMBERHOME/bin/ambpdb -p polyAT_vac.prmtop -c polyAT_gb_init_min.rst > polyAT_gb_init_min.pdb

原始结构(green),和最小化之后的(blue)

savepdb DNA1 XXX.pdb &保存pdb的命令

(2)在归一化的 Born solvation溶剂中运行MD simulation

上一节，我们在真空中运行vacuo，结果不能令人满意，那是因为它不是真实的生物体环境。这一节，我们在Born solvation溶剂再次运行MD。

首先，需要两个类似的输入文件mdin

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

polyAT_gb_md1_nocut.in : no long range cutoff, Generalized Born

polyAT_gb_md1_12Acut.in

10-mer DNA MD Generalise Born, 12A 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

空格

polyAT_gb_md1_nocut.in

10-mer DNA MD Generalise Born, no 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

接下来的仿真，需要用到上面生产的最小化结构文件： polyAT_gb_init_min.rst

MD仿真命令如下：

$AMBERHOME/bin/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

$AMBERHOME/bin/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

这个MD运行时间比较长，每个都需要1-5个多小时。

输入文件是: polyAT_gb_md1_12Acut.in, polyAT_gb_md1_nocut.in, polyAT_gb_init_min.rst, polyAT_vac.prmtop

输出文件: 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

注：mdcrd文件大约15MB大小。

好不容易，运行MD结束，可以进行结果分析了。

(3)结果分析

由于用于结果分析的文件较多，需要创建两个文件夹。

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

注意：

xmgrace ./polyAT_gb_md1_12Acut/summary.EPTOT ./polyAT_gb_md1_nocut/summary.EPTOT

这个命令是用来画势能曲线的，在很多服务器上没有安装这个组件。

因此，我们可以用Matlab或excel来根据summary.EPTOT的数据进行绘制曲线，性质是一样的。势能曲线大致如下：

在接下来，我们分析均方根误差RMSd

首先创建cpptraj 输入文件

polyAT_gb_md1_12Acut.calc_rms文件：

trajin polyAT_gb_md1_12Acut.mdcrd

rms first mass out polyAT_gb_md1_12Acut.rms time 0.1

polyAT_gb_md1_nocut.calc_rms文件：

trajin polyAT_gb_md1_nocut.mdcrd

rms first mass out polyAT_gb_md1_nocut.rms time 0.1

运行cpptraj命令：

$AMBERHOME/bin/cpptraj -p polyAT_vac.prmtop -i polyAT_gb_md1_12Acut.calc_rms

$AMBERHOME/bin/cpptraj -p polyAT_vac.prmtop -i polyAT_gb_md1_nocut.calc_rms

得到两个输出文件，为均方差文件：

polyAT_gb_md1_12Acut.rms

polyAT_gb_md1_nocut.rms

同样可以用Matlab或excel绘制均方差曲线，大致如下

附：最后我们可以在自己的电脑上运行VMD，看DNA仿真运行轨迹。

首先装载文件polyAT_vac.prmtop 并 then select AMBER7 Parm.

Load后，再加载12A的轨迹文件polyAT_gb_md1_12Acut.mdcrd

然后即可播放video进行轨迹观测了。如下图：

Part 5 显式溶剂中最小化并运行MD (explicit solvent)

5.1)平衡并仿真溶剂化的10倍体polyAT

对原有真空中结构进行平衡的原因，及周期性边界条件简介

(1)初步最小化

仿真第一步，建立输入文件mdin

polyAT_wat_min1.in

polyA-polyT 10-mer: initial minimisation solvent + ions

&cntrl

imin = 1,

maxcyc = 4000,

ncyc = 2000,

ntb = 1,

ntr = 1,

cut = 10.0

Hold the DNA fixed

500.0

RES 1 20

END

如果是60bp，则改为RES 1 120，约束原子是1到120(60对)

仿真第二步，sander命令运行最小化

(注意，我们在命令行(-ref)上有一个额外的选项，这其实是指定了需要约束的原子结构，在这种情况下，我们使用的是inpcrd文件中的初始结构。我们将使用在本教程开始时创建的溶剂化prmtop和inpcrd文件。)

$AMBERHOME/bin/sander -O -i polyAT_wat_min1.in -o polyAT_wat_min1.out -p polyAT_wat.prmtop -c polyAT_wat.inpcrd -r polyAT_wat_min1.rst -ref polyAT_wat.inpcrd

输入文件是: polyAT_wat_min1.in, polyAT_wat.prmtop, polyAT_wat.inpcrd

输出文件是: polyAT_wat_min1.out, polyAT_wat_min1.rst

这个最小化过程大约需要3分钟。

如果在这一步结束后，我们查看.out文件，可以发现，van der Waals (VDWAALS and 1-4 VDW)很大，静电能electrostatic energy (EEL)下降很多，说明，水和DNA有不良接触，优化也需要进一步进行。

(2)整个系统最小化

首先需要输入文件mdin

polyAT_wat_min2.in

polyA-polyT 10-mer: initial minimisation whole system

&cntrl

imin = 1,

maxcyc = 4000,

ncyc = 2000,

ntb = 1,

ntr = 0,

cut = 10.0

然后，用sander命令执行最小化：

(注意，这个命令使用了上个步骤的restraint文件.rst，它含有上次最小化后的最终结构，因此，不需要初始的inpcrd结构文件了。)

$AMBERHOME/bin/sander -O -i polyAT_wat_min2.in -o polyAT_wat_min2.out -p polyAT_wat.prmtop -c polyAT_wat_min1.rst -r polyAT_wat_min2.rst

pmemd.cuda -O -i polyAT_wat_min2.in -o polyAT_wat_min2.out -p polyAT_wat.prmtop -c polyAT_wat_min1.rst -r polyAT_wat_min2.rst

输入文件: polyAT_wat_min2.in, polyAT_wat.prmtop, polyAT_wat_min1.rst (上一步的结构文件)

输出文件: polyAT_wat_min2.out, polyAT_wat_min2.rst

这个最小化，大约需要6-20分钟。

这一步最小化后，我们需要看看pdb结构依然稳定，然后再进行MD仿真。

$AMBERHOME/bin/ambpdb -p polyAT_wat.prmtop < polyAT_wat.inpcrd > polyAT_wat.pdb

$AMBERHOME/bin/ambpdb -p polyAT_wat.prmtop < polyAT_wat_min2.rst > polyAT_wat_min2.pdb

生成两个pdb文件：(polyAT_wat.pdb, polyAT_wat_min2.pdb)，咱们可以用VMD观察。

在本机windows上启动VMD

vmd

File -> New Molecule (load the first pdb polyAT_wat.pdb)

(3) 基于溶剂带加热和带约束的MD仿真

这是一个带弱位置约束的20ps的MD仿真。

输入文件polyAT_wat_md1.in (后面我有改写的in文件)

polyA-polyT 10-mer: 20ps MD with res on DNA

&cntrl

imin = 0,

irest = 0,

ntx = 1,

ntb = 1,

cut = 10.0,

ntr = 1,

ntc = 2,

ntf = 2,

tempi = 0.0,

temp0 = 300.0,

ntt = 3,

gamma_ln = 1.0,

nstlim = 10000, dt = 0.002,

ntpr = 100, ntwx = 100, ntwr = 1000

Keep DNA fixed with weak restraints

10.0

RES 1 20

END

Sander命令如下：

$AMBERHOME/bin/sander -O -i polyAT_wat_md1.in -o polyAT_wat_md1.out -p polyAT_wat.prmtop -c polyAT_wat_min2.rst -r polyAT_wat_md1.rst -x polyAT_wat_md1.mdcrd -ref polyAT_wat_min2.rst

pmemd.cuda -O -i polyAT_wat_md1.in -o polyAT_wat_md1.out -p polyAT_wat.prmtop -c polyAT_wat_min2.rst -r polyAT_wat_md1.rst -x polyAT_wat_md1.mdcrd -ref polyAT_wat_min2.rst

采用的是上次第二步最小化产生的约束文件xx_min2.rst

输入文件: polyAT_wat_md1.in, polyAT_wat.prmtop, polyAT_wat_min2.rst

输出文件: polyAT_wat_md1.out, polyAT_wat_md1.rst, polyAT_wat_md1.mdcrd

这个MD过程大约需要30分钟。

(4) 对整个系统运行平衡化的MD 仿真

首先，输入文件设计运行时间是100ps，以便有充足时间relax。

polyAT_wat_md2.in

polyA-polyT 10-mer: 100ps MD

&cntrl

imin = 0, irest = 1, ntx = 7,

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

taup = 2.0,

cut = 10.0, ntr = 0,

ntc = 2, ntf = 2,

tempi = 300.0, temp0 = 300.0,

ntt = 3, gamma_ln = 1.0,

nstlim = 50000, dt = 0.002,

ntpr = 100, ntwx = 100, ntwr = 1000

这样，上次是100帧轨迹文件，这次是500帧，加起来所以是600帧轨迹文件(见下图)。

MD的sander命令如下：

$AMBERHOME/bin/sander -O -i polyAT_wat_md2.in -o polyAT_wat_md2.out -p polyAT_wat.prmtop -c polyAT_wat_md1.rst -r polyAT_wat_md2.rst -x polyAT_wat_md2.mdcrd

注：这次使用的restrt 文件是刚刚MD仿真(不是最小化仿真输出的rst文件)输出的rst文件，因为它含有最终的MD 框架、速度和溶液盒子信息。velocities and box information.

输入文件: polyAT_wat_md2.in, polyAT_wat.prmtop, polyAT_wat_md1.rst

输出文件: polyAT_wat_md2.out, polyAT_wat_md2.rst, polyAT_wat_md2.mdcrd

这里的xx.mdcrd文件很大，约45MB。运行时间大约2小时。如果是20碱基对，运行超过12小时。

MD后台运行方法：

nohup + 命令 + &

再按 ctrl + d 或exit 退出shell就可以了。

5.2) 分析结果以确定稳定性

这一步有许多参数需要分析鉴定，以便确认平衡的质量。包括：

 Potential, Kinetic and Total energy

 Temperature

 Pressure

 Volume

 Density

 RMSd  势能、动能、总能量

 温度

 压力

 体积

 密度

 方差

(1) 分析输出文件

在amber命令行中新建文件夹：

mkdir analysis

cd analysis

process_mdout.perl ../polyAT_wat_md1.out ../polyAT_wat_md2.out

在analysis文件中将输出许多summary文件，比如summary.EPTOT，summary.EKTOT， summary.ETOT等。

可以用matlab 或excel 画出这些文件(都可用ultraedit打开)中的数据曲线，如下图所示。

(2) 分析轨迹

分析DNA的骨架，输入文件是：

polyAT_wat_calc_backbone_rms.in

trajin polyAT_wat_md1.mdcrd

trajin polyAT_wat_md2.mdcrd

rms first out polyAT_wat_backbone.rms @P,O3',O5',C3',C4',C5' time 0.2

计算骨架的RMSd：

$AMBERHOME/bin/cpptraj -p polyAT_wat.prmtop -i polyAT_wat_calc_backbone_rms.in

输出文件：polyAT_wat_backbone.rms

同样可以用matlab画出rms的走势图。

由于前20 ps是有约束的，因此RMSd稳定，而后在去除约束后，RMSd迅速走高。

(3) 可视化MD结果

用本机的VMD软件

 Step 1： load the prmtop file (polyAT_wat.prmtop) ，注意要选择： "AMBER7 Parm"

 Step 2： load the two trajectory files (polyAT_wat_md1.mdcrd, polyAT_wat_md2.mdcrd), 注意：这次选项是 AMBER Coordinates with Periodic Box as file type since our simulations are now periodic boundary runs with the box information stored in the trajectory file.

0ps 120ps

可以发现，在MD后程，因为去除了约束，水跑了。

如果只选中“all not water”可发现(方法如前面的VMD操作所述)，骨架也受到威胁。

(4) 把轨迹重映射到初始盒子

输入文件是：

polyAT_wat_reimage.ptraj (此乃文件名)

trajin polyAT_wat_md1.mdcrd

trajin polyAT_wat_md2.mdcrd

trajout polyAT_wat_md_reimaged.mdcrd

center :1-20

image familiar

用amber命令：

$AMBERHOME/bin/cpptraj -p polyAT_wat.prmtop -i polyAT_wat_reimage.ptraj

生产新的轨迹文件polyAT_wat_md_reimaged.mdcrd(比较大，至少100MB多)

这会将md2轨迹文件附加到md1文件，更改轨迹，使残差1到20的几何中心放置在原点，然后将初始盒子外部的所有坐标映射到盒子box内部。熟悉的关键字确保cpptraj将盒子写成熟悉的截断八面体而不是默认的三斜成像。最后，转换后的轨迹将被写入polyAT_wat_md_reimaged.mdcrd。

现在，咱们可以在VMD中打开此轨迹文件，会发现我们的飞水滴更完整。在vmd中进行表示、去除水并仅观察DNA等。观察在关闭约束后的前20ps后DNA分子的动态如何变化。此外，我们还可以观察到，现在的DNA链比在真空模拟中更稳定。

5.3) Summary

这个教程主要包含了初始PDB文件的创建，以及pdb文件修改(比如加入中性原子，加水盒子等)，以适应LEaP加载(tleap有，xleap没有)。此外，还包含如何中和该DNA结构，如何溶解在水中，如何创建MD输入文件(xx.in)，以及如何采用sander命令运行MD。

需要注意的是，咱们得学会输入文件中每一行参数的具体含义和用法，然后是在真空、隐性溶剂、显性溶剂中运行MD仿真。在该tutorial中，还有最小化、平衡等概念需要掌握。最后，简单的结果分析和轨迹文件可视化也是咱们需要了解的。尤其是VMD软件功能强大，不仅可视化，还可以进行谱密度计算和具有信号处理signal processing功能。

注：

正因为VMD的功能很强大，特别是在谱密度计算和signal processing方面，俺还没有完全搞明白，故借此机会，求助各位大神，如果有用过或知道怎么用，可发帖与大家共享。

Georgecumt 于 Library of Lancaster University, U.K. 2018.11.12.