August 30, 2007

New Blog, Old man......

I register a new blog site on Hi.baidu.com: Sensenbobo Hi Baidu.

You can find something help about basic MD and other thing there, and here will still be my home. ^_^

由 sen 发表于 04:40 PM | 回复 (1)

August 25, 2007

GROMACS程序DEMO例程

####################
### 概述 ###
####################
-----------------------------------------------------------------
-----------------------------------------------------------------
该例程来自Gromacs程序/share/tutor/目录下。整个例程大概只需要十分钟就可以完成，非常适合初学者学习。

该例程是一个完整的分子动力学模拟过程，涵盖了Gromacs程序基本的使用方法。模拟内容是一个水环境下的小肽链。模拟唯一要求是该小肽链的PDB文件。

文档翻译如果有错，请你给我发信：sen@nkbbs.org。相关内容请参阅Gromacs文档，或者给Gromacs开发组询问。

Gromacs官方网站：http://www.gromacs.org
Email: gromacs@gromacs.org
-----------------------------------------------------------------
-----------------------------------------------------------------

#########################
### 环境变量设置 ###
#########################
-----------------------------------------------------------------
-----------------------------------------------------------------
在以下的例程中，所有命令都直接运行，没有添加绝对路径。所以，必须将Gromacs安装路径的bin文件夹加入到系统PATH变量中。如果不加入PATH变量，那么运行时要加入命令的绝对路径。
-----------------------------------------------------------------
-----------------------------------------------------------------

################
### PDB2GMX ###
################
-----------------------------------------------------------------
-----------------------------------------------------------------
在进行如何分子动力学模拟之前，必须建立分子的拓扑文件。Gromacs的分子拓扑文件是用pdb2gmx命令生成，文件后缀名为 .top。攑db2gmx的输入唯一文件是分子的PDB文件，可以从www.pdb.org寻找，文件后缀名 .pdb。”

不是所有的PDB文件都含有氢原子，pdb2gmx将添加所有缺失的氢原子。在pdb2gmx输出的gromacs结构文件中，包含了蛋白质结构的每一个原子，并定义了分子结构的尺寸大小。文件后缀名为 .gro。”
-----------------------------------------------------------------
-----------------------------------------------------------------
假设分子PDB文件名为MOL.pdb，命令运行格式：

pdb2gmx -f MOL.pdb -o MOL.gro -p MOL.top > output.pdb2gmx

命令得到三个文件：MOL.gro是gromacs结构文件；MOL.top分子拓扑文件；output.pdb2gmx是pdb2gmx命令输出文件。

###############
### GENBOX ###
###############
-----------------------------------------------------------------
-----------------------------------------------------------------
真空条件下进行分子模拟输出结果误差较大，所以模拟之前，必须给分子添加水环境。Gromacs中使用genbox命令完成。

Genbox命令读入gromacs的结构文件，并读入水盒子尺寸的大小，输出文件包含分子文件和水盒子。同时genbox更改原来的拓扑文件，使其包含水分子。在使用genbox之前，要使用editconf命令定义水盒子大小。

-----------------------------------------------------------------
-----------------------------------------------------------------
命令格式：

editconf -f MOL.gro -o MOL.gro -d 0.5 > output.genbox

genbox -cp MOL.gro -cs -o MOL_b4em.gro -p MOL.top >> output.genbox


##############
### GROMPP ###
##############
-----------------------------------------------------------------
-----------------------------------------------------------------
原则上讲，到此为止已经可以开始进行模拟计算了。但是，添加在准备分子系统过程中，系统有很多原子距离太近，局部能量太高。这些相互距离太近的原子多是由于genbox程序产生的，溶剂分子与蛋白分子之间存在不稳定的高能量。如果现在开始模拟计算，而不进行能量最优化，系统将可能很不稳定。

去除这些局部高能量的办法是对系统进行能量最优化。能量最优化过程是改变系统中局部高能量的原子的位置，降低这些点的能量。

在进行能量最优化之前，我们先用GROMACS的预处理程序grompp处理所有输入文件。grompp预处理拓扑文件（.top），坐标文件（.gro）和一个参数文件（.mdp），然后输出一个二进制拓扑文件（.tpr）。这个二进制文件包含所有模拟需要的信息，利用这个文件即可以进行能量最优化和动力学模拟。
-----------------------------------------------------------------
em.mdp文件范例，详细请参考Gromacs手册
-----------------------------------------------------------------
title = MOL
cpp = /usr/bin/cpp
define = -DFLEX_SPC
constraints = none
integrator = steep
dt = 0.002 ; ps !
nsteps = 100
nstlist = 10
ns_type = grid
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
;
; Energy minimizing stuff
;
emtol = 1000.0
emstep = 0.01
---------------------------------------------------------------------------------

Grompp命令格式：
grompp -f em -c MOL_b4em -p MOL -o MOL_em >&! output.grompp_em

################
### MDRUN EM ###
################
-----------------------------------------------------------------
-----------------------------------------------------------------
有了二进制拓扑文件，现在可以开始进行能量最优化了。进行能量最优化的程序是mdrun，在
Gromacs程序中，所有模拟都是用mdrun程序进行。

在进行能量最优化过程中，注意查看mdrun程序的输出文件。在输出文件在中，从左到右第一个数字是模拟步数，第二个数字是计算步长，第三个数字是系统能量。如果模拟顺利利，可以看到系统能量从一个很高的值迅速降低，最后稳定在一个大负值。
-----------------------------------------------------------------
-----------------------------------------------------------------
命令格式：
mdrun -nice 4 -s MOL_em -o MOL_em -c MOL_b4pr -v >& ! output.mdrun_em


#################
### GROMPP PR ###
#################
-----------------------------------------------------------------
能量最优化之后，首先保持蛋白质不动，对蛋白质周围水环境进行动力学模拟，该过程称为位置限制性分子动力学（position restrained MD）

位置限制分子动力学保持蛋白质位置不变，对溶剂分子进行平衡计算，可以使溶剂分子填补空间空洞。这个可能存在的空洞是genbox程序产生的。

首先对输入文件进行预处理得到二进制拓扑文件，这些输入文件就是能量最优化得到的输出文件。然后再写一个参数配置文件和索引文件。

默认情况下，程序对模拟系统是分部分的。我们利用两个部分进行位置限制性模拟：Protein和
SOL部分，分别表示蛋白质和溶剂。这这个过程中，即保持protein位置不变。

参数文件（.mdp）包含了位置限制性模拟的所有参数，包括步长，步数，温度等等。同时参数文件告诉GROMACS模拟的类型，如能量最优化、位置限制性或者是分子动力学模拟。
-----------------------------------------------------------------
pr参数文件
-----------------------------------------------------------------
title = MOL position restraining
cpp = /usr/bin/cpp
define = -DPOSRES
constraints = all-bonds
integrator = md
dt = 0.002 ; ps !
nsteps = 500 ; total 1.0 ps.
nstcomm = 1
nstxout = 10
nstvout = 1000
nstfout = 0
nstlog = 10
nstenergy = 10
nstlist = 10
ns_type = grid
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
; Berendsen temperature coupling is on in two groups
Tcoupl = berendsen
tau_t = 0.1 0.1
tc-grps = protein sol
ref_t = 300 300
; Pressure coupling is not on
Pcoupl = no
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Generate velocites is on at 300 K.
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529

----------------------------------------------------------------------------------------
命令格式

grompp -f pr -c MOL_b4pr -r MOL_b4pr -p MOL -o MOL_pr >& ! output.grompp_pr

################
### MDRUN PR ###
################
-----------------------------------------------------------------
-----------------------------------------------------------------
现在进行1ps的位置限制性分子动力学模拟（the Position restrained Molecular Dynamics simulation）在这里进行1ps只是为了节省例程演示实践时间，在实际研究中，1ps模拟将太短，不符合实际情况。

-----------------------------------------------------------------
-----------------------------------------------------------------
命令格式
mdrun -nice 4 -s MOL_pr -o MOL_pr -c MOL_b4md -v >& ! output.mdrun_pr

#################
### GROMPP MD ###
#################
-----------------------------------------------------------------
-----------------------------------------------------------------
到此为止，分子系统已经可以进行正式的分子动力学模拟了。再一次利用grompp程序将输入文件转化成二进制拓扑文件（.tpb/.tpr文件后缀）。
分子动力学模拟配置文件（参数详细意义请参考Gromacs文档）
-----------------------------------------------------------------
title = MOL MD
cpp = /usr/bin/cpp
constraints = all-bonds
integrator = md
dt = 0.002 ; ps !
nsteps = 5000 ; total 5 ps.
nstcomm = 1
nstxout = 50
nstvout = 0
nstfout = 0
nstlist = 10
ns_type = grid
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
; Berendsen temperature coupling is on in two groups
Tcoupl = berendsen
tau_t = 0.1 0.1
tc-grps = protein sol
ref_t = 300 300
; Pressure coupling is not on
Pcoupl = no
tau_p = 0.5
compressibility = 4.5e-5
ref_p = 1.0
; Generate velocites is on at 300 K.
gen_vel = yes
gen_temp = 300.0
gen_seed = 173529
----------------------------------------------------------------

命令格式：

grompp -f md -c MOL_b4md -p MOL -o MOL_md >& ! output.grompp_md

################
### MDRUN MD ###
################
-----------------------------------------------------------------
-----------------------------------------------------------------
现在进行分子动力学模拟计算。详细观察模拟步数的增加。（整体模拟步数为2500，时间长度为5ps。再一次强调，模拟时间太短，仅适合教学演示。）
-----------------------------------------------------------------
-----------------------------------------------------------------
mdrun -nice 4 -s MOL_md -o MOL_md -c MOL_after_md -v >& ! output.mdrun_md

############
### NGMX ###
############
-----------------------------------------------------------------
-----------------------------------------------------------------
最后，可以使用Gromacs的ngmx程序观看动力学模拟轨迹（文件后缀.trj）。轨迹文件包含了分子动力学过程的坐标，原子速度和受力情况等等。可供分析使用。Ngmx的详细使用方法，请参考Gromacs文档。由于分子动力学模拟分析方法各异，这里不再讨论。

本文来自gromacs程序教学例程，仅供参考。^_^

由 sen 发表于 11:09 AM | 回复 (3)

August 08, 2007

分子动力学模拟概论

五、经验力场

对生物分子的理论研究主要是从原子水平探究其结构、功能与动力学之间的关系。这样的问题涉及大量原子，目前还无法直接使用量子力学解决。但是如果使用经验力场，问题则可以很好的解决。所谓经验力场，简单的说就是使用一系列势能方程描述原子之间的相互作用。使用经验力场解决生物分子中各个原子之间的相互作用，可以节约大量计算时间，但是是牺牲一定准确性为代价的。

经验力场发展到今天，已经很好的在计算精确度和计算效率之间找到一个平衡。这些经验力场一般经过一定使用时间后还不停与实验结果进行校准与修改。经验力场对物理实验的再现能力主要包括：X光衍射和NMR得到的生物结构信息、光谱学方法和无弹性中子散射得到的动力学信息、统计热力学信息等等。经验力场方程的各个参数的校订是一件费时费力的工作，需要不停的大量优化。在过去二十年中，经验力场研究已经形成一个独立的领域，不同的研究小组还不停的加入到这个领域中去。现在比较著名的经验力场有AMBER力场、CHARMM力场、GROMOS力场和OPLS/AMBER力场等。这些力场的持续发展重视基础研究同时更重视实际应用。

上文提到，使用经验力场有一定的限制性。经验力场最明显也最重要的一个限制性就是模拟的生物分子不能发生重大构象变化，如化学键的断裂和重建。为了突破这样的限制性，新型的经验力场已经开始引进量子力学机制。这一类力场还需要一定的完善时间，如有兴趣可以参阅其他资料。

经验力场之间一般大同小异，并一般公开使用，如NAMD软件使用CHARMM力场，也可以使用AMBER力场；DOCK实用部分AMBER力场最为评分函数。在此简单介绍CHARMM22力场，该力场可以用于模拟蛋白(MacKerell et al 1998)、核酸(MacKerell et al. 1995)、流体(Schlenkrich et al. 1996)和碳水化合物(Ha et al. 1988)等。

References
A. D. MacKerell Jr., J. Wiórkiewicz-Kuczera, and M. Karplus (1995). An All-Atom Empirical Energy Function for the Simulation of Nucleic Acids, J. Am. Chem. Soc. 117, 11946-11975.

A. D. MacKerell, Jr., D. Bashford, M. Bellott, R. L. Dunbrack Jr., J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, III, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin, M. Karplus (1998). All-atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. 102, 3586-3616.

M. Schlenkrich, J. Brickmann, A. D. MacKerell Jr., and M. Karplus (1996). An Empirical Potential Energy Function for Phospholipids: Criteria for Parameter Optimization and Applications, in Biological Membranes: A Molecular Perspective from Computation and Experiment, K. M. Merz, Jr. and B. Roux, editors (Birkhauser), 31-81.

S.N. Ha, A. Giammona, M. Field, J.W. Brady, A revised potential-energy surface for molecular mechanics studies of carbohydrates, Carbohydr. Res. (1988), 180(2), 207-221.

CHARMM经验力场

一个系统的能量E，可以表达为一个关于系统中所有原子笛卡尔坐标R的方程。系统能量值主要由两个部分组成：分子内部化学键、键角、二面角相关的能量值Ebonded；外部非成键相互作用Enon-bonded，如静电和范德华作用等。




其中Ebonded部分主要有三个能量相组成：



分别代表键长畸变、键角畸变和二面角畸变。




上述等式第一项键长畸变能是以弹簧振子弹性势能方程近似表达的，其中b0表示成键原子之间的平衡位置，Kb为力学常数。这两个参数都是根据原子性质不同而不同的。






力学常数常常需要根据红外线频率延展实验或者量子力学计算进行调整，而键长平衡值则多依据高分辨率晶体结构或者微波光谱分析数据进行调整。

第二项键角畸变是指化学键角偏离平衡角度θ0所具有的弹性势能。平衡键角θ0和力学常熟Kθ根据不用化学键角不同而不同。






以上两项能量都是分子结构偏离平衡位置所产生的势能，在一个平衡分子体系中，该两项理想状态下总合应该接近零。

第三项是涉及二面角畸变所具有的势能，包含三个连续化学键连接的四个原子。二面角畸变能是化学键旋转造成的，势能变化具有周期性，一般使用一个余弦方程表示。





除了以上三项，CHARMM力场还包含其它两项与化学键相关的能量项：
i Urey-Bradley作用能：相隔两个化学键原子之间的相互作用；
ii 非正常二面角作用能：用于手性分子以及保持原子构成的平面。
而力学常数：Kb、Kq、Kf主要是对大量小分子复合物研究得来，并同时与量子力学计算结果相比较调整。


ChARMM力场中对非键对相互作用对势能的影响的考虑主要分两个部分：范德华相互作用能和静电相互作用能。在其他经验力场中，有的还独立计算氢键对势能的影响，而CHARMM力场则将氢键包含在范德华相互作用和静电相互作用中。



两个原子之间的范德华相互作用力主要是在两个原子距离小于平衡距离时的排斥力和距离大于平衡距离时表现出来的吸引力。排斥力与吸引力都与原子的电子云分布有关，详细内容可以参考相关资料。当排斥力与吸引力大小相同，互相抵消的时候，两原子的距离为平衡距离，总体势能最小。一个系统中原子之间都处于平衡位置时，系统最为稳定。系统最低能量值和原子之间平衡距离大小取决于不同原子类型。



范德华相互作用是保持生物分子系统稳定最为重要的因素，其计算模型一般为Lennard-Jones 6-12势能模型：



常量A和C主要取决于原子类型，可以通过如气相光散射的实验测量。
原子之间的静电相互作用可以使用库仑电势能表示：



其中D为有效电介质常数，r为两个原子电量中心之间qi和qk之间的距离。静电相互作用根据系统中各个原子的坐标位置可以各个求值，同时可以确定各个原子的受力情况以应用于分子动力学模拟计算。

经验力场是一种能量近似，存在一定不精确的局限性。

经验力场中各个原子的属性是同一的，原子属性用于定义与其它原子成键或者相互影响时表现的特性。但是，比如一个在sp3电子成键状态下成键的C&alpha原子与一个HIS残基环结构上的C原子的性质是不同，这不是经验力场所能描述的。除了对各个原子的原子属性进行统一处理，经验力场还对不同原子进行分类，以尽可能减少原子类型，简化计算，这回导致原子类型识别错误。某一些原子，如脂肪质中的碳原子或者生物分子中的氢原子对周围环境比较不敏感，这类原子的统一属性定义一般不会出现较明显的误差；但是如氧原子和氮原子与周围其它原子影响十分强烈，这一类原子属性需要根据成键环境不同而不同。

为了降低计算资源的要求，经验力场对生物系统中势能计算使用原子势能成对叠加近似，即系统中一个原子与周围的相互作用是该原子与其他各个原子相互作用的叠加，两个原子自间的相互作用忽略其他原子的影响。三个及三个原子以上之间的同时间的相互作用在经验力场中是不予考虑的，所以极化作用并不能精确的表述。这些都有可能使计算结果与实验结果存在差异。

另外重要一点是经验力场势能方程并没有考虑熵效应，所以势能的最小值并不一定与系统稳定状态或者最有可能处在的状态互相对应，这些状态在实验条件下是与系统势能最低态相对应的。

由 sen 发表于 10:23 AM | 回复 (0)