分子动力学模拟I.doc

Gromacs中文教程淮海一粟 分子动力学（MD）模拟分为三步：首先，要准备好模拟系统；然后，对准备好的系统进行模拟；最后，对模拟结果进行分析。虽然第二步是最耗费计算资源的，有时候需要计算几个月，但是最耗费体力的步骤在于模拟系统准备和结果分析。本教程涉及模拟系统准备、模拟和结果分析。一、数据格式处理 准备好模拟系统是MD最重要的步骤之一。MD模拟原子尺度的动力学过程，可用于理解实验现象、验证理论假说，或者为一个待验证的新假说提供基础。然而，对于上述各种情形，都需要根据实际情况对模拟过程进行设计；这意味着模拟的时候必须十分小心。丢失的残基、原子和非标准基团 本教程模拟的是蛋白质。首先需要找到蛋白质序列并选择其起始结构，见前述；然后就要检查这个结构是否包含所有的残基和原子，这些残基和原子有时候也是模拟所必需的。本教程假定不存在缺失，故略去。 另一个需要注意的问题是结构文件中可能包含非标准残基，被修饰过的残基或者配体，这些基团还没有力场参数。如果有这些基团，要么被除去，要么就需要补充力场参数，这牵涉到MD的高级技巧。本教程假定所有的蛋白质不含这类残基。结构质量 对结构文件进行检查以了解结构文件的质量是一个很好的练习。例如，晶体结构解析过程中，对于谷氨酰胺和天冬酰胺有可能产生不正确的构象；对于组氨酸的质子化状态和侧链构象的解析也可能有问题。为了得到正确的结构，可以利用一些程序和服务器(如 WHATIF)。本教程假定所用的结构没有问题，我们只进行数据格式处理。二、结构转换和拓扑化 一个分子可以由各个原子的坐标、键接情况与非键相互作用来确定。由于.pdb结构文件只含有原子坐标，我们首先必须建立拓扑文件，该文件描述了原子类型、电荷、成键情况等信息。拓扑文件对应着一种力场，选择何种力场对于拓扑文件的建立是一个值得仔细考虑的问题。这里我们用的是GROMOS96 53a6连接原子力场，该力场对于氨基酸侧链的自由能预测较好，并且与NMR试验结果较吻合。 拓扑数据要与结构吻合，这点很重要；这意味着结构数据也需要在所用的力场中进行转换。 为了转换结构并建立拓扑数据，可以用pdb2gmx程序，该程序对特定结构块（如氨基酸）建立拓扑数据。它需要使用结构块库，如果结构里面有库中没有的结构单元，转换过程就会失败。输入以下命令来进行结构转换，选择GROMOS 53a6力场和SPC水分子模型. 选项 ignh表示起始文件中的氢原子会先被除掉，然后在相应的力场中重建。pdb2gmx -f protein.pdb -o protein.gro -p protein.top -ignh You can always generate a pdb file from a gro file using the following command:editconf -f XXXXX.gro -o XXXXX.pdb 仔细阅读屏幕提示，并检查对于组氨酸质子化所做的选择，注意蛋白质总电荷数；也要浏览输入文件（protein.pdb, pdb格式）和输出文件（protein.gro, gromacs格式）。注意两者的区别；最后浏览拓扑文件及其结构。三、结构的能量最小化（真空中） 现在，一个适合于所选力场的正确结构文件已经建立好，同时还得到了相应的拓扑数据。然而，由于结构转换中，牵涉到氢原子的删除和重建，这可能会带来一些相互作用力的变化，比如由于两个原子的位置靠得太近引起的作用力。因此我们必须对结构进行一次能量优化，使之松弛一点。这其实就是一个模拟步骤，包括两个过程：第一步，结构和拓扑数据被合成在一起，同时还包含了一些控制参数。这个过程对产生一个文件，该文件作为输入文件，用于第二步mdrun程序的动力学模拟。 把这个过程分为两步的好处是，输入文件可以传送到专门的（高性能）计算机网络或者超级计算机，在那里完成模拟计算。然而，一般只在正式动力学模拟中才会这么做。 为了执行第一步，需要用到一个.mdp文件，该文件（文件名minim.mdp）已经被做好并放在你们的目录下面。该文件包含一系列能量最小化的控制参数，请查看之。注意，文件以integrator行开始，表示指定的数学模型。本例中，被指定采用最速下降法计算5000步。现在，用以下命令合并结构和拓扑文件，以及一些控制参数：grompp -v -f minim.mdp -c protein.gro -p protein.top -o protein-EM-vacuum.tpr grompp程序将会标明此体系存在非零电荷，并将写入有关系统和控制参数的其他信息。该程序还会产生一个额外的输出文件(mdout.mdp)，该文件包含所有的控制参数设置. 下一步，运行mdrun：mdrun -v -deffnm protein-EM-vacuum -c protein-EM-vacuum.gro 由于该蛋白只含有1500个原子，能量优化非常快。选项 -v 将把每步的势能和最大势能力打印出来，以便于跟踪能量优化过程。-deffnm选项设定了输出文件的默认文件名，而不需要对每个输出文件都进行命名，以减少选项设置。现在，结构松弛下来了，我们该设定周期性边界并添加溶剂。Which method was used for the energy 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 energy minimization by loading them into PyMol四、周期性边界条件设定 在添加溶剂之前，需要选定模拟体系的大致外形。大多数常见的模拟都是在周期性边界条件(PBC)下进行的，这意味着要定义一个外形单元框，该框可用于空间填充堆积。这样，就可以对一个无限、周期性系统进行模拟，以避免由于边框造成的边界效应。建立PBC时，仅有少数几个通用的形状。我们将采用rhombic dodecahedron，因为其对应于最优堆积空间，因此是自由旋转原子的最佳选择。为了不产生周期性影像的直接接触，我们设定了蛋白质距离边框的最小距离为1.0，这样的话，两个相邻的堆积单元的距离就会大于2.0 nm。周期性边界条件用editconf命令设定： editconf -f protein-EM-vacuum.gro -o protein-PBC.gro -bt dodecahedron -d 1.2 看看editconf的输出，注意体积的变化。同时，也看看文件protein-PBC.gro的最后一行。在gromacs格式文件 (.gro)中，最后一行表示溶剂盒子的形状。我们常常用三斜晶系矩阵表示法，其中前三个数字表示对角线元素 (xx, yy, zz)，最后六个表示非对角线元素(xy, xz, yx, yz, zx, zy)。What is the new unit cell volume? 五、溶剂条件设定 溶剂盒子设定好了之后，就可添加溶剂了。溶剂模型有好几种，每种模型多少都与一种力场相对应。GROMOS96力场常用于SPC水分子模型。溶剂添加不需要写入拓扑数据，但是也可以在拓扑数据中包含溶剂模型。向盒子中添加溶剂用以下命令：genbox -cp protein-PBC.gro -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 correspond to? How many atoms are in the system now? 用以前学过的命令，将 .gro文件转换为 .pdb文件。用Pymol程序读入溶剂化的结构文件 (protein-water.pdb)。Pymol中，可用“show cell”命令显示盒子形状。show cell 这将画出一个框架线来显示盒子的三维空间边界。这个三斜晶图形对应菱形十二面体，但可能不太明显。另一个值得注意的事情是，并非所有的溶剂都在盒子内，但以砖形排布。非常类似地，蛋白质也有部分伸出水外面。 Why is it not a problem to have the protein sticking out 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程序的输入文件。我们设定了NA+/CL-(-pname NA+ -nname CL-)的浓度为0.15 M (-conc 0.15)，并指定特定离子的加入以是体系电中性 (-neutral)。Genion程序将询问用离子取代何种分子，选择SOL。genion -s protein-water.tpr -o protein-solvated.gro -conc 0.15 -neutral -pname NA+ -nname CL- 注意添加离子的数量，并注意一个额外的氯离子被加入，以中和体系电荷。通过将水分子置换为离子，体系的拓扑数据文件（protein.top）将不正确了（你可以修改拓扑文件，减少溶剂分子数目；同时加入一行命令，说明添加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?七、溶剂化体系的能量最小化 到此为止，模拟系统已经准备好了。在进行正式模拟前，还需要使体系再次松弛。因为溶剂的添加以及离子的置换，可能产生不利的相互作用力，例如，原子重叠、同种电荷的接近，等，因此需要对体系再次能量优化，其步骤与前面相同：先运行grompp，再运行mdrun：grompp -v -f minim.mdp -c protein-solvated.gro -p protein.top -o protein-EM-solvated.tpr mdrun -v -deffnm protein-EM-solvated -c protein-EM-solvated.gro How many steps were specified in the parameter file and how many steps did the minimization take? What is the final potential energy of the system? 八、溶剂和氢原子位置的松弛：位置限制性MD 为了缓解体系的随机分布的张力，我们是溶剂与蛋白质相适应，比如，我们使溶剂自由移动，而使得蛋白质上面的非氢原子或多或少地固定在相对位置。这么做的目的是为了 确保溶剂的构象“适合”蛋白质。这一步是正式MD的第一步。这一步的控制参数在pr.mdp文件中，执行时需要被读入。看看这个文件， 注意integrate行和define行。后者用于允许拓扑文件的流动控制（to allow flow control in the topology file）： define = -DPOSRES 将导致定义关键字POSRES。查看拓扑文件的最后一行，看看怎样进行位置限制。采用 grompp和mdrun程序运行模拟：Remark: edit the file pr.mdp and change the value of gen_seed grompp -v -f pr.mdp -c protein-EM-solvated.gro -p protein.top -o protein-PR.tpr mdrun -v -deffnm protein-PR What is the length of the simulation in picoseconds? How is the inclusion of the position restraints handled in the topology file? 你会很有趣地看到总能量、势能和动能的变化。在前一步，粒子没有速度，所以没有动能，因此也没有温度。现在，模拟一开始，原子就被赋予速度因此获得动能。模拟过程的能量信息被保存为不可读的二进制文件格式（.edr文件）。该文件的信息可用gromacs工具g_energy程序提取。 注：在下一步模拟运行时，你能用以下命令并在执行过程中回答有关问题，以节省时间。如果这一步不奏效，让该程序继续在这里运行，而打开另外一个终端，运行下一步的grompp和 mdrun程序，但要注意两个终端的都要在正确的文件夹下运行。g_energy -f protein-PR.edr -o thermodynamics-PR.xvg 你将要在一列能量参数中进行选择，选择好的参数结果将会被输出。输入与温度、势能、动能和总能量相对应的数字，然后输入0（零），回车。 输出文件 (.xvg) 可以用xmgr或xmgrace程序查看，这些程序需要预先安装好。也可以用Python语言脚本xvg2ascii.py查看，它能在终端窗口上显示出图线。xmgrace -nxy thermodynamics-PR.xvg 注：绿色曲线是总能量，黑色是势能，蓝色是温度，红色是动能。你也可以点击曲线并填写“图标”窗口以改变图标。别忘了点击 accept 以接受改变。What happens with the temperature? What happens to the potential/kinetic/total energy and how can this be explained? 九、非限制性MD，打开温度耦合 系统现在应该相对松弛了一些，我们就引入温度并开始对系统进行热浴耦合吧。换句话说，将在打开温度耦合的情况下运行一小段时间，让系统松弛到一个新的状态。下载控制参数文件nvt.mdp，然后检查一下其中的温度耦合参数；最后用grompp和mdrun运行模拟：grompp -v -f nvt.mdp -c protein-PR.gro -p protein.top -o protein-NVT.tpr mdrun -v -deffnm protein-NVT There are two groups which are separately coupled to a heat bath. Which groups are that and what do you think is in each of these groups? 运行结束后，查看一下能量和温度。再次提醒，你可以在查看能量的同时，运行下一步 。数据提取方法同前：g_energy -f protein-NVT.edr -o thermodynamics-NVT.xvg 观察能量图。What happens to the temperature? What happens to the potential/kinetic/total energy and how can this be explained? 十、非限制性MD模拟：打开压力耦合 当引入温度并松弛后，就该引入压力了。导入控制参数文件npt.mdp，短时间运行有压力下的模拟。先查看一下此控制参数文件，找到控制引入压力的句柄。此处也应该注意，因为我们需要重新赋予离子速度，通过gen_vel行实现。在这行下面，设定了速度分布（Maxwell分布）的温度，并通过gen_seed产生初始值。gen_seed现在设定为9999，但是会被改为一个新的值。 MD模拟是最重要的一步，所以每次都以相同的坐标体系、速度和相同的参数开始模拟时，就会导致模拟的结果一模一样，这不是我们所希望的。 当你设定好了控制参数后，运行以下命令进行模拟。Remark: edit the file npt.mdp and change the value of gen_seed grompp -v -f npt.mdp -c protein-NVT.gro -p protein.top -o protein-NPT.tpr mdrun -v -deffnm protein-NPT 运行完成后，再次查看能量和温度，同时注意观察压力变化。数据提取方法同前：g_energy -f protein-NPT.edr -o thermodynamics-NPT.xvg 看看能量曲线、温度和压力的变化。What happens to the temperature? What happens to the potential/kinetic/total energy and how can this be explained? 最后一个问题与从有限数量的粒子系统中提取热动力特性直接有关。粗略地说，热动力学指的是大量粒子（例如，几十亿个而不是几千个）的行为。取大量粒子的平均特性能够使数据波动减少，相反，仅对少量粒子进行平均，波动就较大。 由于我们引入了压力，系统的密度会发生变化，所以从能量文件中提取密度数据，方法如下： g_energy -f protein-NPT.edr -o density-NPT.xvg How does the density behave over time? Why does the density of the system change if pressure coupling is turned on? 十一、正式模拟 到了最后