分子动力学模拟I.docx
- 文档编号:407361
- 上传时间:2022-10-09
- 格式:DOCX
- 页数:11
- 大小:101.43KB
分子动力学模拟I.docx
《分子动力学模拟I.docx》由会员分享,可在线阅读,更多相关《分子动力学模拟I.docx(11页珍藏版)》请在冰豆网上搜索。
分子动力学模拟I
Gromacs中文教程
淮海一粟
分子动力学〔MD〕模拟分为三步:
首先,要准备好模拟系统;然后,对准备好的系统进展模拟;最后,对模拟结果进展分析。
虽然第二步是最消耗计算资源的,有时候需要计算几个月,但是最消耗体力的步骤在于模拟系统准备和结果分析。
本教程涉与模拟系统准备、模拟和结果分析。
一、数据格式处理
准备好模拟系统是MD最重要的步骤之一。
MD模拟原子尺度的动力学过程,可用于理解实验现象、验证理论假说,或者为一个待验证的新假说提供根底。
然而,对于上述各种情形,都需要根据实际情况对模拟过程进展设计;这意味着模拟的时候必须十分小心。
丢失的残基、原子和非标准基团
本教程模拟的是蛋白质。
首先需要找到蛋白质序列并选择其起始结构,见前述;然后就要检查这个结构是否包含所有的残基和原子,这些残基和原子有时候也是模拟所必需的。
本教程假定不存在缺失,故略去。
另一个需要注意的问题是结构文件中可能包含非标准残基,被修饰过的残基或者配体,这些基团还没有力场参数。
如果有这些基团,要么被除去,要么就需要补充力场参数,这牵涉到MD的高级技巧。
本教程假定所有的蛋白质不含这类残基。
结构质量
对结构文件进展检查以了解结构文件的质量是一个很好的练习。
例如,晶体结构解析过程中,对于谷氨酰胺和天冬酰胺有可能产生不正确的构象;对于组氨酸的质子化状态和侧链构象的解析也可能有问题。
为了得到正确的结构,可以利用一些程序和服务器(如WHATIF)。
本教程假定所用的结构没有问题,我们只进展数据格式处理。
二、结构转换和拓扑化
一个分子可以由各个原子的坐标、键接情况与非键相互作用来确定。
由于.pdb结构文件只含有原子坐标,我们首先必须建立拓扑文件,该文件描述了原子类型、电荷、成键情况等信息。
拓扑文件对应着一种力场,选择何种力场对于拓扑文件的建立是一个值得仔细考虑的问题。
这里我们用的是GROMOS9653a6连接原子力场,该力场对于氨基酸侧链的自由能预测较好,并且与NMR试验结果较吻合。
拓扑数据要与结构吻合,这点很重要;这意味着结构数据也需要在所用的力场中进展转换。
为了转换结构并建立拓扑数据,可以用pdb2gmx程序,该程序对特定结构块〔如氨基酸〕建立拓扑数据。
它需要使用结构块库,如果结构里面有库中没有的结构单元,转换过程就会失败。
输入以下命令来进展结构转换,选择GROMOS53a6力场和SPC水分子模型.选项–ignh表示起始文件中的氢原子会先被除掉,然后在相应的力场中重建。
仔细阅读屏幕提示,并检查对于组氨酸质子化所做的选择,注意蛋白质总电荷数;也要浏览输入文件〔protein.pdb,pdb格式〕和输出文件〔protein.gro,gromacs格式〕。
注意两者的区别;最后浏览拓扑文件与其结构。
三、结构的能量最小化〔真空中〕
现在,一个适合于所选力场的正确结构文件已经建立好,同时还得到了相应的拓扑数据。
然而,由于结构转换中,牵涉到氢原子的删除和重建,这可能会带来一些相互作用力的变化,比如由于两个原子的位置靠得太近引起的作用力。
因此我们必须对结构进展一次能量优化,使之松弛一点。
这其实就是一个模拟步骤,包括两个过程:
第一步,结构和拓扑数据被合成在一起,同时还包含了一些控制参数。
这个过程对产生一个文件,该文件作为输入文件,用于第二步mdrun程序的动力学模拟。
把这个过程分为两步的好处是,输入文件可以传送到专门的〔高性能〕计算机网络或者超级计算机,在那里完成模拟计算。
然而,一般只在正式动力学模拟中才会这么做。
为了执行第一步,需要用到一个.mdp文件,该文件〔文件名c:
\iknow\docshare\data\cur_work\md.chem.rug.nl\~mdcourse\minim.mdp〕已经被做好并放在你们的目录下面。
该文件包含一系列能量最小化的控制参数,请查看之。
注意,文件以"integrator"行开始,表示指定的数学模型。
本例中,被指定采用最速下降法计算5000步。
现在,用以下命令合并结构和拓扑文件,以与一些控制参数:
grompp程序将会标明此体系存在非零电荷,并将写入有关系统和控制参数的其他信息。
该程序还会产生一个额外的输出文件(mdout.mdp),该文件包含所有的控制参数设置.下一步,运行mdrun:
由于该蛋白只含有1500个原子,能量优化非常快。
选项-v将把每步的势能和最大势能力打印出来,以便于跟踪能量优化过程。
-deffnm选项设定了输出文件的默认文件名,而不需要对每个输出文件都进展命名,以减少选项设置。
现在,结构松弛下来了,我们该设定周期性边界并添加溶剂。
Whichmethodwasusedfortheenergyminimization?
Howmanystepswerespecifiedintheparameterfileandhowmanystepsdidtheminimizationtake?
Whatcouldcausetheenergyminimizationtostopbeforethespecifiednumberofstepswasreached?
Whatisthefinalpotentialenergyofthesystem?
parethestructuresbeforeandafterenergyminimizationbyloadingthemintoPyMol
四、周期性边界条件设定
在添加溶剂之前,需要选定模拟体系的大致外形。
大多数常见的模拟都是在周期性边界条件(PBC)下进展的,这意味着要定义一个外形单元框,该框可用于空间填充堆积。
这样,就可以对一个无限、周期性系统进展模拟,以防止由于边框造成的边界效应。
建立PBC时,仅有少数几个通用的形状。
我们将采用rhombicdodecahedron,因为其对应于最优堆积空间,因此是自由旋转原子的最优选择。
为了不产生周期性影像的直接接触,我们设定了蛋白质距离边框的最小距离为1.0,这样的话,两个相邻的堆积单元的距离就会大于2.0nm。
周期性边界条件用editconf命令设定:
看看editconf的输出,注意体积的变化。
同时,也看看文件的最后一行。
在gromacs格式文件(.gro)中,最后一行表示溶剂盒子的形状。
我们常常用三斜晶系矩阵表示法,其中前三个数字表示对角线元素(xx,yy,zz),最后六个表示非对角线元素(xy,xz,yx,yz,zx,zy)。
Whatisthenewunitcellvolume?
五、溶剂条件设定
溶剂盒子设定好了之后,就可添加溶剂了。
溶剂模型有好几种,每种模型多少都与一种力场相对应。
GROMOS96力场常用于SPC水分子模型。
溶剂添加不需要写入拓扑数据,但是也可以在拓扑数据中包含溶剂模型。
向盒子中添加溶剂用以下命令:
看看的结尾溶剂添加的情况,添加了多少溶剂?
Howmanywatermoleculesareaddedtothesystemandwhatvolumedoesthatcorrespondto?
Howmanyatomsareinthesystemnow?
用以前学过的命令,将.gro文件转换为.pdb文件。
用Pymol程序读入溶剂化的结构文件(protein-water.pdb)。
Pymol中,可用“showcell〞命令显示盒子形状。
这将画出一个框架线来显示盒子的三维空间边界。
这个三斜晶图形对应菱形十二面体,但可能不太明显。
另一个值得注意的事情是,并非所有的溶剂都在盒子内,但以砖形排布。
非常类似地,蛋白质也有局部伸出水外面。
Whyisitnotaproblemtohavetheproteinstickingoutofthewater?
六添加相反的离子以平衡体系电荷
目前我们已经有了溶剂化的蛋白质了,但是体系的净电荷数不为零。
为了使体系电中性,我们需要添加一定数量的相反离子。
添加一定离子的步骤,是一个很好的练习,程序genion能完成这些任务,但是它需要一个输入文件,例如,一个包含了结构和拓扑数据的文件。
就像前面一样,我们可以用grompp来产生这样的文件:
文件既可作为genion程序的输入文件。
我们设定了NA+/CL-(-pnameNA+-nnameCL-)的浓度为0.15M(-conc0.15),并指定特定离子的参加以是体系电中性(-neutral)。
Genion程序将询问用离子取代何种分子,选择'SOL'。
注意添加离子的数量,并注意一个额外的氯离子被参加,以中和体系电荷。
通过将水分子置换为离子,体系的拓扑数据文件〔〕将不正确了〔你可以修改拓扑文件,减少溶剂分子数目;同时参加一行命令,说明添加Na离子的数量并写入另一行命令说明Cl离子的数量〕。
Howmanysodiumandchlorideionsareaddedtothesystem?
七、溶剂化体系的能量最小化
到此为止,模拟系统已经准备好了。
在进展正式模拟前,还需要使体系再次松弛。
因为溶剂的添加以与离子的置换,可能产生不利的相互作用力,例如,原子重叠、同种电荷的接近,等,因此需要对体系再次能量优化,其步骤与前面一样:
先运行grompp,再运行mdrun:
Howmanystepswerespecifiedintheparameterfileandhowmanystepsdidtheminimizationtake?
Whatisthefinalpotentialenergyofthesystem?
八、溶剂和氢原子位置的松弛:
位置限制性MD
为了缓解体系的随机分布的X力,我们是溶剂与蛋白质相适应,比如,我们使溶剂自由移动,而使得蛋白质上面的非氢原子或多或少地固定在相对位置。
这么做的目的是为了确保溶剂的构象“适合〞蛋白质。
这一步是正式MD的第一步。
这一步的控制参数在c:
\iknow\docshare\data\cur_work\md.chem.rug.nl\~mdcourse\pr.mdp文件中,执行时需要被读入。
看看这个文件,注意integrate行和define行。
后者用于允许拓扑文件的流动控制〔toallowflowcontrolinthetopologyfile〕:
define=-DPOSRES将导致定义关键字"POSRES"。
查看拓扑文件的最后一行,看看怎样进展位置限制。
采用grompp和mdrun程序运行模拟:
Whatisthelengthofthesimulationinpicoseconds?
Howistheinclusionofthepositionrestraintshandledinthetopologyfile?
你会很有趣地看到总能量、势能和动能的变化。
在前一步,粒子没有速度,所以没有动能,因此也没有温度。
现在,模拟一开始,原子就被赋予速度因此获得动能。
模拟过程的能量信息被保存为不可读的二进制文件格式〔.edr文件〕。
该文件的信息可用gromacs工具g_energy程序提取。
注:
在下一步模拟运行时,你能用以下命令并在执行过程中回答有关问题,以节省时间。
如果这一步不奏效,让该程序继续在这里运行,而打开另外一个终端,运行下一步的grompp和mdrun程序,但要注意两个终端的都要在正确的文件夹下运行。
你将要在一列能量参数中进展选择,选择好的参数结果将会被输出。
输入与温度、
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 分子 动力学 模拟