NAMD入门教程(三)

4 受控分子动力学模拟(Steered Molecular Dynamics)所谓受控分子动力学模拟（Steered Molecular Dynamics，SMD），就是指在进行分子动力学模拟时，人为地给分子中的某个或某几个原子施加一个假想的外力，或者人为地固定某个或某几个原子的位置。从而控制整个分子的行为。我们前面进行的动力学模拟都是通过各种参数设置，尽量逼近蛋白质分子在溶液体系中的真实状态，以研究分子的各种性质和行为。但受控分子动力学模拟却要用一个假想的外力干扰控制生物大分子的行为，这是为什么呢？（需要查一下更多的SMD的应用实例，详细了解SMD究竟可以用于研究什么问题，才能做出回答）在下面的例子中，我们将首先固定泛素分子中一个原子的位置，然后用一个假想的外力牵拉另一个原子。球形的泛素分子会因此而被逐渐拉开，最终成为伸展状态的肽链。进行SMD时，我们需要用已经完成能量最小化和能量平衡，达到稳定状态的蛋白结构。如果用含有扭曲、拉伸、变形构象的原始蛋白结构，我们将无法分清蛋白各部分的运动是由于蛋白内部的形变张力引起的，还是由于SMD实验附加的外力造成的。因此，进行SMD之前必须预先进行一次平衡态分子动力学模拟，获得稳定结构。这里，我们使用2.4节球状水体分子动力学模拟输出的恢复文件（.restart）进行SMD。恢复文件输出时，泛素已在水体中完成了能量最小化和能量平衡，达到了稳定状态。4.1 除去水分子为了节省计算时间，我们在进行本次SMD之前将除去体系中所有的水分子。但读者必须要注意：在真正进行SMD实验的时候决不可以将水分子除去！1、打开VMD，选择FileNew Molecule菜单项，载入common目录下的文件ubq_ws.psf。不要关闭窗口，此时窗口“Load file for:”一项应该显示“0:ubq_ws.psf”。再次单击按钮Browser，找到1-2-sphere目录下的文件ubq_ws_eq.restart.coor，载入该文件。关闭Molecule File Browser窗口。现在我们载入了2.4节球状水体动力学模拟时输出的恢复文件(.restart.coor)。2、选择Extension tk菜单项，打开tk，首先用cd命令改变当前目录至common目录下，然后输入：set selprotein atomselect top protein$selprotein writepdb ubq_ww_eq.pdb这样我们便在common文件夹下新建了一个pdb文件ubq_ww_eq.pdb，储存已经达到平衡态的蛋白质分子，但没有水分子。3、删除当前分子，但不要关闭VMD。4.2 恒速牵拉图 恒速牵拉的图示。图中SMD原子是蓝色的，假原子（dummy atom）是红色的，二者之间有一根假想的弹簧相连。假原子以恒定的速度运动，弹簧发生形变后作用于SMD原子，于是SMD原子开始运动。在这一节我们将使用恒速牵拉进行SMD动力学模拟。进行恒速牵拉时，我们需要首先设定一个假想的原子作为施力原子，这个实际并不存在的原子就叫做假原子（dummy atom）。假原子通过一根假想的弹簧与系统中真实存在的某个原子（称为SMD原子）相连，牵拉时，假原子以恒定的速度运动，因此称这种牵拉方式为恒速牵拉。实际上被牵拉原子（SMD原子）的速度不是恒定的。SMD原子所受拉力由弹簧的形变和弹性系数决定，满足胡克定律：。又有：其中是SMD原子所受的力，是假想弹簧的形变量。它等于假原子的位移减去SMD原子的位移，其中是SMD原子的初始位置，是SMD原子的当前位置。因此，SMD原子的受力可以计算出来，然后根据牛顿定律计算它的运动情况。4.2.1 设定SMD原子和固定原子（Fixed atom）NAMD使用pdb文件中的B因子一栏区别哪些原子是被固定的，哪些原子未被固定：如果某个原子在该栏中的对应值是1（非0值即可），那么该原子被固定，在整个动力学模拟过程中，它的空间坐标不会变化；反之如果为0，那么该原子不受影响。知识链接：PDB文件的格式使用写字板打开common目录下我们刚刚制作的ubq_ww_eq.pdb，可以看到文件的内容分成几栏（column）：在上图中，给我们提供有用信息的共有10栏数据，已在图中标注出：（1）表示所指示为原子（2）该原子序列号（3）IUPAC标准格式的原子名称（4）残基名称（5）残基序列号（6）原子的X坐标（7）原子的Y坐标（8）原子的Z坐标（9）位置（Occupancy）（10）温度因子或B因子（beta factor）其中后两栏在动力学模拟的时候被NAMD用于标记两类特殊的原子：位置（Occupancy）栏用于标记SMD原子，B因子一栏用于标记固定原子（fixed atom）。4.2.1 设定固定原子（fixed atom）和SMD原子下面我们将使用VMD创建我们所需的pdb文件。我们在本例中将使用pdb文件中B因子一栏（即温度因子）指定哪个原子将被固定，使用位置（occupancy）一栏指定哪个原子作为SMD原子被牵拉。1、在VMD中选择FileNew Molecule 菜单项，单击Browse按钮找到common目录下的文件ubq.psf，载入它。不要关闭Molecule File Browser窗口，注意Load file for 一项应当显示1:ubq.psf。再次单击Browse按钮，载入common目录下的ubq_ww_eq.pdb。关闭Molecule File Browse窗口，在VMD图形窗口中应当可以看到没有水分子的泛素分子（图）。图 载入除去水分子后的泛素分子2、在VMD tk中输入以下命令，固定第一个氨基酸残基的碳：set allatoms atomselect top all$allatoms set beta 0set fixedatom atomselect top “resid 1 and name CA”$fixedatom set beta 1下面是对命令的解释：¨ set allatoms atomselect top all 新建了一个变量allatoms，它代表体系中的所有原子。¨ $allatoms set beta 0 所有原子在pdb文件中的B栏（即beta）全部设为0。¨ set fixedatom atomselect top “resid 1 and name CA” 新建一个变量fixedatom，代表第一个氨基酸残基的碳。¨ $fixedatom set beta 1 fixedatom所代表的原子第一个氨基酸残基的碳的B栏被设置为1。这样NAMD会在分子动力学模拟时保持该原子固定。同样地，我们还需要设定哪一个原子将被牵拉。前面提到过，被牵拉的原子就叫做SMD原子。对于这一种原子，NAMD使用pdb文件中的“位置”（occupancy）一栏进行标记。同样地，0代表不受影响的原子，标记为1的原子将被牵拉。设置方法仍然是在tk中输入命令进行。$allatoms set occupancy 0set smdatom atomselect top “resid 76 and name CA”$smdatom set occupancy 1以下是对上述命令的解释：¨ $allatoms set occupancy 0 当前pdb文件中所有原子的“occupancy”一栏设置为0。¨ set smdatom atomselect top “resid 76 and name CA” 新建了一个变量smdatom，代表第76个氨基酸（即最后一个氨基酸）的碳。¨ $smdatom set occupancy 1 变量smdatom所代表的原子的occupancy 一栏设定为1，这样NAMD会在分子动力学模拟时将假想的力作用于该原子上。完成上述设定之后，我们需要让VMD输出我们设置好的文件。如果当前目录不是common，用cd命令改变当前目录到common中，然后输入：$allatoms writepdb ubq_ww_eq.ref 回车后会在当前目录生成文件ubq_ww_eq.ref，这就是我们进行SMD所需的文件，它记录了泛素的结构，以及哪一个原子作为固定原子（fixed atom），哪一个作为牵拉原子。不要关闭VMD，将所有窗口最小化。使用写字板打开common目录下的文件ubq_ww_eq.ref，注意打开时“文件类型”依然要选择“所有文件（*.*）”否则将看不到该文件。在文件开头第六行即可找到固定原子第一个氨基酸甲硫氨酸的碳。这个原子的B因子一栏为1.00，而其它原子都为0.00（图）。拖动滚动条向下，一直到倒数第四行，可以看到我们定义的SMD原子最后一个氨基酸甘氨酸的碳。这个原子的Occupancy一栏是1.00，而其它原子都是0.00（图）。图 固定原子（Fixed atom）的设置图 SMD原子的设置需要说明一点：在上面两张图中出现的原子坐标数值（第7、8、9栏数值）可能和读者看到的不一致。这是因为动力学模拟中，能量平衡（equilibration）一步中各个原子的初始速度是随机选择的。我们使用的pdb文件是各个原子能量平衡结束后的末速度，因为初速度是随机的，末速度也不一定相同。4.2.2 设定拉力的方向现在，关闭写字板。我们已经设定了固定原子（fixed atom）和SMD原子。下面我们需要做的是设定拉力的方向。我们定义拉力的方向为沿固定原子（fixed atom）和SMD原子的连线方向（即由固定原子到SMD原子的矢量的方向）。在VMD tk中输入：set smdpos lindex $smdatom get x y z 0set fixedpos lindex $fixedatom get x y z 0vecnorm vecsub $smdpos $fixedpos输入以上命令后，tk中给出三个数值，这三个值就是固定原子（fixed atom）到SMD原子的矢量坐标，分别为x，y，z坐标。记下这三个值（图），我们后面还会用到。然后关闭VMD即可。图 假想的拉力矢量的坐标4.2.3 配置文件现在我们已经得到了文件：ubq_ww_eq.ref，这个文件储存了原子的位置，以及哪个原子是固定原子，哪个原子是SMD原子。下一步我们将制作NAMD配置文件。注意在下面的步骤中一定要避免输入错误，否则很可能导致动力学模拟无法正常进行。我们将对一个样本配置文件进行修改，制作我们所需要的配置文件。在实际工作中，读者也可以修改NAMD教程中提供的样本配置文件，得到自己所需要的文件。1、使用Windows资源管理器获得我们所需的样本配置文件sample.conf。文件在common目录下。将这一文件粘贴到3-1-pullcv目录下,然后改名为ubq_ww_pcv.conf以和其他文件相区分。2、使用写字板打开ubq_ww_pcv.conf。现在我们看到的是一个配置文件的样本（图）。浏览一下可以发现很多参数是用户需要根据实际情况进行设定的。事实上，这个样本配置文件就是NAMD动力学模拟配置文件的基本蓝本，用户可以在此基础上量身定做符合自己需要的配置文件。图 NAMD配置文件样本3、在Job Description 一栏加入:# N- C- Termini Constant Velocity Pulling当然这只是一个注释， 不加也不会影响我们