《分子动力学模拟I》由会员分享，可在线阅读，更多相关《分子动力学模拟I（13页珍藏版）》请在金锄头文库上搜索。
1、Gromacs 中文教程淮海一粟分子动力学（MD）模拟分为三步：首先，要准备好模拟系统；然后，对准备好的系 统进行模拟；最后，对模拟结果进行分析。虽然第二步是最耗费计算资源的，有时候需要 计算几个月，但是最耗费体力的步骤在于模拟系统准备和结果分析。本教程涉及模拟系统 准备、模拟和结果分析。一、数据格式处理准备好模拟系统是MD最重要的步骤之一。MD模拟原子尺度的动力学过程，可用 于理解实验现象、验证理论假说，或者为一个待验证的新假说提供基础。然而，对于上述 各种情形，都需要根据实际情况对模拟过程进行设计；这意味着模拟的时候必须十分小丢失的残基、原子和非标准基团本教程模拟的是蛋白质。首先需要找到蛋
2、白质序列并选择其起始结构，见前述；然 后就要检查这个结构是否包含所有的残基和原子，这些残基和原子有时候也是模拟所必需 的。本教程假定不存在缺失，故略去。另一个需要注意的问题是结构文件中可能包含非标准残基，被修饰过的残基或者配 体，这些基团还没有力场参数。如果有这些基团，要么被除去，要么就需要补充力场参 数，这牵涉到MD的高级技巧。本教程假定所有的蛋白质不含这类残基。结构质量对结构文件进行检查以了解结构文件的质量是一个很好的练习。例如，晶体结构解 析过程中，对于谷氨酰胺和天冬酰胺有可能产生不正确的构象；对于组氨酸的质子化状态 和侧链构象的解析也可能有问题。为了得到正确的结构，可以利用一些程序和服
3、务器（如 WHATIF）。本教程假定所用的结构没有问题，我们只进行数据格式处理。二、结构转换和拓扑化个分子可以由各个原子的坐标、键接情况与非键相互作用来确定。由于.pdb结构 文件只含有原子坐标，我们首先必须建立拓扑文件，该文件描述了原子类型、电荷、成键 情况等信息。拓扑文件对应着一种力场，选择何种力场对于拓扑文件的建立是一个值得仔细考虑的问题。这里我们用的是GROMOS96 53a6连接原子力场，该力场对于氨基酸侧 链的自由能预测较好，并且与NMR试验结果较吻合。拓扑数据要与结构吻合，这点很重要；这意味着结构数据也需要在所用的力场中进 行转换。为了转换结构并建立拓扑数据，可以用pdb2gmx
4、程序，该程序对特定结构块（如氨 基酸）建立拓扑数据。它需要使用结构块库，如果结构里面有库中没有的结构单元，转换 过程就会失败。输入以下命令来进行结构转换，选择GROMOS 53a6力场和SPC水分子 模型选项-gnh表示起始文件中的氢原子会先被除掉，然后在相应的力场中重建。仔细阅读屏幕提示，并检查对于组氨酸质子化所做的选择，注意蛋白质总电荷数； 也要浏览输入文件（protein.pdb, pdb格式）和输出文件（protein.gro, gromacs格式）。 注意两者的区别；最后浏览拓扑文件及其结构。三、结构的能量最小化（真空中）Cnntr d IP QramfiTe ri EMVqcju
5、mCnntr d IP QramfiTe ri EMVqcju m现在，一个适合于所选力场的正确结构文件已经建立好，同时还得到了相应的拓扑 数据。然而，由于结构转换中，牵涉到氢原子的删除和重建，这可能会带来一些相互作用 力的变化，比如由于两个原子的位置靠得太近引起的作用力。因此我们必须对结构进行一 次能量优化，使之松弛一点。这其实就是一个模拟步骤，包括两个过程：第一步，结构和 拓扑数据被合成在一起，同时还包含了一些控制参数。这个过程对产生一个文件，该文件 作为输入文件，用于第二步mdrun程序的动力学模拟。把这个过程分为两步的好处是， 输入文件可以传送到专门的（高性能）计算机网络或者超级计算机
6、，在那里完成模拟计 算。然而，一般只在正式动力学模拟中才会这么做。为了执行第一步，需要用到一个.mdp文件，该文件（文件名minim.mdp）已经被做 好并放在你们的目录下面。该文件包含一系列能量最小化的控制参数，请查看之。注意， 文件以integ rato r行开始，表示指定的数学模型。本例中，被指定采用最速下降法计算 5000步。现在，用以下命令合并结构和拓扑文件，以及一些控制参数：grompp -v -f mini m mdp -c protein gro -p protein.top -o protein-EM-vacuum.tprgrompp程序将会标明此体系存在非零电荷，并将写入有
7、关系统和控制参数的其他信 息。该程序还会产生一个额外的输出文件（mdout.mdp），该文件包含所有的控制参数设置. 下一步，运行mdrun :mdrun -v -deffnm protein-EM-vacuum -c protein-EM-vacuum gro由于该蛋白只含有1500个原子，能量优化非常快。选项-v将把每步的势能和最大 势能力打印出来，以便于跟踪能量优化过程。-deffnm选项设定了输出文件的默认文件 名，而不需要对每个输出文件都进行命名，以减少选项设置。现在，结构松弛下来了，我 们该设定周期性边界并添加溶剂。Which method was used for the ene
8、rgy minimization?How many steps were specified in the parameter file and how many steps did the minimization take?What could cause the energy minimization to stop before the specified number of steps was reached?What is the final potential energy of the system?Compare the structures before and after
9、 energy minimization by loading them into PyMol四、周期性边界条件设定在添加溶剂之前，需要选定模拟体系的大致外形。大多数常见的模拟都是在周期性 边界条件(PBC)下进行的，这意味着要定义一个外形单元框，该框可用于空间填充堆积。 这样，就可以对一个无限、周期性系统进行模拟，以避免由于边框造成的边界效应。建立 PBC时，仅有少数几个通用的形状。我们将采用rhombic dodecahedron，因为其对应于 最优堆积空间，因此是自由旋转原子的最佳选择。为了不产生周期性影像的直接接触，我 们设定了蛋白质距离边框的最小距离为1.0,这样的话，两个相邻的堆积
10、单元的距离就会 大于2.0 nm。周期性边界条件用editconf命令设定：看看editconf的输出，注意体积的变化。同时，也看看文件protein-PBC.gro的最后 一行。在gromacs格式文件(.gro)中，最后一行表示溶剂盒子的形状。我们常常用三斜晶系矩阵表示法，其中前三个数字表示对角线元素(XX, yy, zz),最后六个表示非对角线元素 (xy, xz, yx, yz, zx, zy)oWhat is the new unit cell volume?五、溶剂条件设定S true tur* eW i t h VVa! erTop ohgylnc ludingYVnter溶剂盒
11、子设定好了之后，就可添加溶剂了。溶剂模型有好几种，每种模型多少都与 一种力场相对应。GROMOS96力场常用于SPC水分子模型。溶剂添加不需要写入拓扑数 据，但是也可以在拓扑数据中包含溶剂模型。向盒子中添加溶剂用以下命令：genbox -cp protein-PBCgro -cs spc216 gro -p protein.top -o protein-water.gro看看protein.top的结尾溶剂添加的情况，添加了多少溶剂？How many water molecules are added to the system and what volume does that corres
12、pond to?How many atoms are in the system now?用以前学过的命令，将.gro文件转换为.pdb文件。用Pymol程序读入溶剂化的结 构文件(protein-water.pdb)。Pymol中，可用show cel”命令显示盒子形状。show cell这将画出一个框架线来显示盒子的三维空间边界。这个三斜晶图形对应菱形十二面 体，但可能不太明显。另一个值得注意的事情是，并非所有的溶剂都在盒子内，但以砖形 排布。非常类似地，蛋白质也有部分伸出水外面。Why is it not a problem to have the protein sticking ou
13、t of the water?六添加相反的离子以平衡体系电荷目前我们已经有了溶剂化的蛋白质了，但是体系的净电荷数不为零。为了使体系电 中性，我们需要添加一定数量的相反离子。添加一定离子的步骤，是一个很好的练习，程 序genion能完成这些任务，但是它需要一个输入文件，例如，一个包含了结构和拓扑数 据的文件。就像前面一样，我们可以用grompp来产生这样的文件：grompp -v -f minim mdp -c protein-water.gro -p protein.top -o protein-water.tpr文件protein-water.tpr既可作为genion程序的输入文件。我们设
14、定了 NA+/CL-(- pname NA+ -nname CL-)的浓度为0.15 M (-cone 0.15)，并指定特定离子的加入以是体系 电中性(-neutral)。Genion程序将询问用离子取代何种分子，选择SOL。genion -s protein-water.tpr -o protein-solvatedgro -conc 015 -neutral -pname NA+ -nname CL-注意添加离子的数量，并注意一个额外的氯离子被加入，以中和体系电荷。通过将水 分子置换为离子，体系的拓扑数据文件(protein.top)将不正确了(你可以修改拓扑文 件，减少溶剂分子数目；同
15、时加入一行命令，说明添加Na离子的数量并写入另一行命令 说明Cl离子的数量)。Edit the topology file and decrease the number of solvent molecules Also add a line specifying the number of NA ions and a line specifying the amount of CL.How many sodium and chloride ions are added to the system?七、溶剂化体系的能量最小化Cq ntro I Pnramet ersEArt So I ven
16、tSy st e mConf i gum tio n5ystEmTDpologySIMULATION十efti到此为止，模拟系统已经准备好了。在进行正式模拟前，还需要使体系再次松弛。 因为溶剂的添加以及离子的置换，可能产生不利的相互作用力，例如，原子重叠、同种电 荷的接近，等，因此需要对体系再次能量优化，其步骤与前面相同：先运行grompp,再 运行mdrun :grompp -v -f minim.mdp -c protein-solvatedgro -p protein.top -o protein-EM- solvated.tprmdrun -v -deffnm protein-EM-solvated -c protein-EM-solvated
