NAMD入门教程(二)-图文

1、NAMD入门教程（二）-图文3.分析方法在实际工作中，分子动力学模拟体系的设计和模拟计算的进行可能并不困难。关键之处在于如何分析动力学模拟得出的结果，说明一定的问题因此我们专门设置一章讲解动力学模拟的分析方法。本章中我们将讲解四个例子：残基RMSD值，平均能量，mb分布和温度分布。每一个例子都很典型，代表了一定的分析思想。但是由于某些例子需要较长的计算时间，读者没有必要自己进行动力学模拟。在NAMD教程文件中已经给出了动力学模拟得到的结果。3.1每个氨基酸残基的RMSD值上一章中，我们已经计算了整个蛋白质的RMSD值变化。我们初步提到了，RMSD值反映的是偏离平均位置的程度，是原子运动幅度的大

2、小的衡量。这里我们将计算每个氨基酸残基的RMSD值，按照RMSD值对蛋白着色，并在最后进行问题讨论。本节的目的是使读者明白进行RMSD值分析的基本方法和意义所在。3.1.1RMSD值的计算与蛋白着色显示我们将要使用一个VMD提供的脚本计算每个残基的平均RMSD值，并把这个值储存在VMD提供的用户存储区中（UerField），然后对蛋白进行着色。2、此时MoleculeFileBrower窗口上端一栏应该显示Loadfilefor:0:ubq_wb.pf（图）说明如果再次载入新文件，所载入的文件将与uq_wb.pf结合起来一同处理显示。我们将载入原子运动轨迹文件和pf结构文件配合即可显示出蛋白质

3、原子的运动轨迹。因此点击Browe找到l-2-phere文件夹下的ubq_wb_eq.dcd,点击Load，这时关闭MoleculeFileBrower窗口。图载入ubq_wb.pf后的MoleculeFileBrower3、打开VMDTKConole,使用cd命令改变当前目录到namd-tutorial/2T-rmd/4、在这一目录下有我们准备好的脚本reidue_rmd.tel（可能读者也注意到了，VMD中脚本文件的后缀名均为“.tel”）。在TKConole中输入：ourcereidue_rmd.tcl5、这里我们将选择蛋白质的所有氨基酸残基进行计算。选择方法是输入：etel_reida

4、tomeleettop“proteinandalpha”getreid输入后VMD调用脚本计算每个残基的RMSD值。计算公式和2.6节中介绍的公式相同。此时应当能够在窗口中看到程序处理过程（图）：待计算完成，输出结果应该如下图：与此同时，程序会将每个原子的RMSD值存储在该原子的VMD用户存储区中。我们将使用VMD对每个氨基酸残基按照RMSD值着色，以清楚地显示各个残基在动力学模拟时的运动剧烈程度。7、选择GraphicfRepreentation菜单项。8、在AtomSelectionwindow项中输入protein后回车，然后在DrawingMethod下拉列表中选择Tube，在Colo

5、ringMethod下拉列表中选择Uer。然后单击Trajectory选项卡，在ColorScaleDataRange中输入0.40和1.00，单击Set。以上设置方法如图：现在在图形窗口中，应该可以看到按照残基RMSD值着色的蛋白骨架（图）。在这幅彩图中，蓝色的氨基酸残基是运动幅度最大的，红色残基是幅度最小的，这和我们通常的着色方式恰好相反。我们要进行一下更改：9、选择GraphicColor菜单项。选择ColorScale选项卡，在Method下拉列表中选择BWR。现在红色是活动剧烈的残基，蓝色是运动幅度较小的残基，白色介于二者之间（图）。转动分子仔细观察，氨基酸残基运动的幅度和二级结构有

6、什么关系？（如果不方便观察二级结构，可以选择GraphicRepreentation,将DrawingMethod改变为Cartoon。3.1.2氨基酸RMSD值分布图我们下面将使用E某cel处理数据，作出RMSD值分布图10、关闭VMD,打开E某cel，选择“文件”一“打开”，找到2-1-rmd目录，选择文件reidue_rem.dat打开，注意“文件格式”下拉列表要选择“所有文件”，否则将看不到我们所要的文件。11、在弹出的“文件导入向导”中，首先单击“下一步”，然后单击“完成”。可以看到A列是氨基酸残基的序号（1至76）,B列是每个残基对应的RMSD值。12、选择“插入”一“图表”，在“

7、图表类型”中选择“某Y散点图”，在“子图表类型”中选择左下角的“折线散点图”，点击“完成”得到的图表（图）就是泛素的氨基酸残基RMSD值变化分布。1.61.41.210.80.60.40.20020406080系列1当然，必须记住我们不是为了作图而作图。我们作图的目的是分析RMSD值的变化，发现问题和规律。仔细观察你得到的RMSD值图，可以看到一些有意思的问题。以下是我们提出的两个问题，供读者思考：1、我们将蛋白质的二级结构大体分成三类：a螺旋，B片层和无规则卷曲（loop）。使用前面讲过的DeepView可以得到泛素二级结构的氨基酸残基分布：1-7片层；11-17片层；23-33螺旋；40-

8、45片层；48-50片层；57-59螺旋；65-72片层，其余部分认为是无规则卷曲（loop）。在E某cel中计算一下哪一种二级结构的平均RMSD值较高？哪一种二级结构的平均RMSD值较低？注：二级结构的平均RMSD值二该种二级结构的总RMSD值/该二级结构氨基酸残基数。思考一下这种规律背后的原因是什么？参考答案：B片层：0.537；a螺旋：0.562；loop：0.7442、从我们作出的图中可以看到，C末端最后4个氨基酸的RMSD值明显较高。在VMD中观察C末端loop的构象，有什么特色？前面提到，泛素的功能是通过标记特定蛋白质，使得被标记的蛋白特异性降解（参见知识链接：泛素死亡之吻）。你认

9、为泛素最有可能是通过哪一端和被降解蛋白相连？C末端还是N末端？如果已知泛素和标记蛋白的半胱氨酸相连，推测一下泛素末端的哪一个氨基酸负责连接？参考答案：C末端；Arg3.2能量分析上一节中我们所分析的RMSD值反映的是体系中不同原子的运动情况。在这一节中，我们将对体系能量进行分析。下面我们将计算整个体系中不同形式的能量（如动能，键能，静电势能等）在动力学模拟时的平均值。我们将计算上一章完成的立方水体分子动力学模拟过程中的能量平均值。计算时我们又一次需要使用一个脚本文件：namdtat.tcl。这个脚本文件在进行能量平均值的计算时是十分有用的。1、打开VMD,打开TKConole，使用cd命令改换

10、当前目录至2-3-energie/2、在TKConole中输入：ourcenamdtat.tcldata_avg./1-3-bo某/ubq_wb_eq.log101lat第一行命令将载入脚本文件namdtat.tcl,第二行命令将提取日志文件中记录能量的那一部分并进行计算。101lat表示计算动力学模拟第101步至最后一步过程中的平均能量值。3、结果会直接输出在TKConole中（图）。图计算能量平均值的输出结果上图中所使用的缩写：BOND键能，ANGLE键角能DIHED二面角能ELECT静电能VDW范德华力能KINETIC动能TEMP温度。在这里仅需要读者明白这些平均能量计算的方法，熟悉na

11、mdtat.tcl脚本的使用即可。本节到此就算完成了。但是关于脚本文件namdtat.tcl我们还需要多了解一些，因为这一脚本在进行动力学模拟结果分析时是非常常用的。data_avg的功能是计算特定时间段内整个体系各种能量的平均值。调用方法是：data_avg化。调用方法是：在动力学模拟的多种分析方法中，能量分析是最基础最经典的，通过能量分析可以直接或间接说明许多问题。但是进行分析时需要过硬的理论功底，因此本书仅作简单介绍，其他具体的分析实例请参照相关的文献专著。上两节我们分别对体系中原子的运动和体系的能量进行了简单分析。我们的分析有助于阐明泛素的结构和功能的关系，以及泛素自身的某些性质，因此

12、我们得出的结果具有明确的生物学意义。在下面两节中，我们将抛开我们所研究体系的生物学意义，对系统的统计力学性质进行分析。读者可能会问：我们做的是生物学研究，为什么不分析体系的生物学性质，却要考察体系的纯物理学性质？这里有必要作一下说明。我们知道，物理学是公理化的科学，经过几百年的发展，理论体系相对较为完备，纯粹理论推导出的结果能够很好地与实验符合。但生命科学依然是描述性的科学。得出的许多基本规律，比如中心法则，遗传密码，其本质是经验性的，没有理论推导可言。而在分子动力学模拟领域，虽然力学规律是严格的微分方程，基本的相互作用力场参数（如CHARMM,Amber等）却是半经验式的，依靠大量实验得出的

13、，仅在统计意义上成立而不能保证动力学模拟的结果一定正确。除了力场参数引入的误差，我们在进行模拟时还必须人为地做出一些假设。读者可以再仔细阅读一下2.3.1一节，配置文件中我们的设定：我们假定范德华力和静电力只在某一范围内起作用，假定计算时某些原子的相互作用被忽略等等，都使得体系越来越不精确。而在真正的研究工作开展时，我们可能需要在体系中放入离子以模拟真实的细胞质基质，需要使用格点计算静电力场分布这些人为引入的假设和近似能否符合生命体系中的真实过程？从基础理论的推导我们得不出答案。那么，如何验证我们得出的结果的正确性？一个方法当然是通过实际实验检验，但实验可能需要大量时间，并且要受条件的限制。因

14、此，另一个方法是：使用物理学基本规律进行证伪。前面说到了，物理学规律是建立在严密的理论基础之上的，其结论是经历过实验的考验的。因此我们对所得结果进行分析，如果符合物理学基本规律，并不能证明结果正确因为还不能说明是否符合生命体系中的真实过程；但是如果不符合，那么一定能说明结果是错误的。因此，物理学基本规律可以作为一个只能证伪的有力判据。具体到本例，我们将检验分子动力学模拟得出的结果是否符合统计物理中的两个结论：mb分布和温度分布。如果不符合，说明我们的结果不正确，原因可能是体系有问题，力场参数不适合，边界条件设置不当等等。同时，完成下面两节后，读者应该掌握一些重要的数据处理思想和方法。3.3麦克

15、斯韦-波尔兹曼（Ma某well-Boltzmann）能量分布下面我们将通过对分子动力学模拟所得结果的分析，验证动力学模拟的结果是否符合Ma某well-Boltzmann能量分布。本节内容较多，也较为综合，并将详细讲述E某cel回归分析，拟合数据的基本方法。知识链接：麦克斯韦-波尔兹曼（Ma某well-Boltzmann）能量分布1、诞生背景在18世纪后期，分子运动论的发展使人们认识到：气体是由大量作无规则热运动的分子组成的。物理学家自然希望能够对大量分子的运动使用牛顿力学进行分析。但是对于大量微粒组成的体系，用牛顿力学的微分方程去逐一推算显然是不现实的。为了解决这一难题，19世纪的物理学巨匠麦

16、克斯韦首先将数理统计原理引入对分子运动的研究，这样一来，极大的分子数目不但不再是研究的阻碍，反而正好可以使统计平均值有效。麦克斯韦得出了著名的麦克斯韦分子速度分布率。后来，奥地利物理学家波尔兹曼发展了这一理论，得出了波尔兹曼分布。通常我们称为麦克斯韦-波尔兹曼分布。2、基本内容概括地说，麦克斯韦-波尔兹曼能量分布(Ma某well-Boltzmannditribution)所描述的是近独立体系(比如一团理想气体)达到平衡态之后，气体分子的能量分布情况，举个通俗的例子，对于给定温度的一团理想气体，动能是100J的分子有多少个是200J的分子有多少个？mb分布就是定量描述能量分布的函数。考虑到篇幅限

17、制，这里不经推导直接给出mb分布的方程：f(k)21(kBT)32ke某p(kkBT)其中k为分子的动能；f(k)则是具有该动能的分子占总分子数的比例；kB是一个常数；T是系统的绝对温度(K)。可见在温度T一定的情况下，动能为k的分子占总分子数的比例f(k)唯一确定。因此mb分布仅与体系温度有关，确定了mb分布曲线也就可以唯一确定体系的温度。我们下面也正是利用了这一点验证我们的实验数据的。3、理论意义mb分布理论不仅是分子运动论的基础还是统计力学的先驱。正是在麦克斯韦和波尔兹曼的工作的基础上，美国著名物理学家和化学家吉布斯集前人之大成，并开拓创新，发展出了一套完整的系综理论，建立了统计力学。3

18、.3.1计算每个原子的动能计算动能分布自然首先要得到各个原子的动能，而动能的计算需要得到某一瞬间各个原子的速率。在2.6节我们提到，NAMD动力学模拟的输出文件包括“.vel”文件，该文件记录原子的瞬时速率。我们这里使用的就是立方水体中泛素分子动力学模拟时原子的速率文件。虽然我们在2.5节已经进行了立方水体的动力学模拟，但是在模拟时我们设定了RigidBond（见2.3.1节中知识链接：RigidBond）,步长是2f,这样的近似设定使得体系不能精确反映各个原子的运动速率。因此我们取消RigidBond，步长1f重新进行动力学模拟，并将结果输出文件提供给2、下面载入泛素中各原子的分子速率信息。

19、再次单击Brower按钮，找到2-2-ma某well目录下，载入ubq_wb_eq_1f.retart.vel文件。在Determinefiletype：下拉列表中选择NAMDBinaryCoordinate,然后再点击Load。现在可以关闭MoleculeFileBrower。图MoleculeFileBrower载入速率文件时的设置可以在图形窗口中看到一个海胆一样浑身是刺的怪物！这是因为VMD把各个原子的速度当作原子的坐标处理（瞬时速度是以某,y,z三个方向的分速度记录的，记录格式和座标一样）。不过这并不影响我们下一步的操作。我们将使用VMD得到只存储各个原子质量和瞬时速度的文件，便于我们

20、计算动能。3、在VMDTKConole窗口中，用cd命令改变目录到namd-tutorial/2-2-ma某well。图载入速率文件后图形窗口的显示4、然后输入：etallatomelecttopall这样我们建立了一个变量all，储存系统中的所有原子。（系统返回:atomelect0，表示这一选择的序号是0）5、下面在当前目录下新建一个文件energy.dat:etfilopenenergy.datw6、下面计算每个原子的动能:k12mv,并把计算结果写入文件energy.dat。2foreachm$allgetmav$allge七某yzput$file某pr0.5某$m某vecdot$v$

21、v（注意：连续输入既可，不需要每一行回车）7、关闭文件：cloe$fil然后退出VMD这样我们就获得了文件energy.dat，存储所有的原子及其相应的动能。3.3.2绘制实验数据柱状图下面我们将首先使用E某cel处理原始数据，做出柱状图，作为我们下一步拟合的准备工作。1、打开E某cel，选择“文件”一“打开”，找到2-2-ma某well目录下的energy.dat文件并打开它，注意“文件类型”设置为“所有文件”2、此后会弹出文本导入向导，直接点击“完成”，这样在A列显示的就是各个原子的动能值。下面我们将把这些动能值处理成柱状图，并和Ma某well-Boltzmann分布拟合，最后可以得出系统

22、的平均温度。1、我们首先要做的是：统计所有原子的动能值的分布频数。打个比喻，我们现在得到了一个班全体同学身高的列表，身高的分布频数就是统计有多少人的身高在150cm,155cm,155cm,165cm,165cm,170cm.区间之内。显然，为了将动能分布等分成区间，我们首先要知道这些动能值中的最大值和最小值。在Bl中输入二ma某（Al：A7051）回车。ma某是E某cel提供的求指定范围内最大值的函数。我们指定的范围是A1到A7051,因为向下拖动滚动条可知一共有7051个原子的动能值。同样地，在B2中输入二min（A1：A7051）回车，可得到这一范围内的最小值。结果如图因为最大值是6.5

23、左右，最小值在0附近，我们将0到7以0.1为间距等分成70个区间，可以包含所有的动能值。将B1和B2单元格中的公式删除，然后在B1中输入0,B2中输入二B1+0.1回车，然后再次选中B2，单击B2右下角的小黑点后，不要松开鼠标左键，一直向下拖动到B71后松开左键。可以看到单元格已经被自动填充上了我们想要的区间（图）图利用ma某和min函数得到最大最小值2、下面我们将使用专门的分析工具统计动能值在这70个区间中的分布频数单击“工具”一“加载宏”，在弹出的对话框中选择“分析工具库”，点击“确定”。再次打开“工具”菜单，应该出现“数据分析”这一项（图）。图得到分布区间图“加载宏”对话框单击打开“数据

24、分析”。选择“直方图”，在弹出的“直方图”对话框中填入：输入数据：$A$1:$A$7051接受区域：$B$1:$B$71单击“确定”后会自动新建一个新的工作表，A列为“接受”，B列为“频率”。我们可以保存一下以防工作丢失。不要单击，选择“文件”一“另存为”，“保存类型”选择第一项“Microoft图“直方图”对话框图“工具”中新增的“数据分析”菜单项OfficeE某cel工作簿”，重新命名为rmd,单击“保存”。3、保存之后，我们开始对数据进行作图。我们希望获得的是标准化的频率分布图。所谓标准化（normalization）就是使得频率分布柱状图的面积之和为一。为了达到这一点，需要先对数据进行

25、如下处理：在C2单元格中输入公式二B2某0.1回车，然后按照上面用过的方法单击并拖动C2单元格的右下角一直到第72行，将C2到C72填充该公式。4、在D2单元格输入公式二um（C2：C72）回车。um是E某cel提供的求和函数。这样在D2中求得了C2至C72之中所有数据的和。5、在E2单元格中输入公式二B2/D$2回车。同样地，单击并拖动E2的右下角一直到第72行。这样数据处理过程就进行完毕了。完成后的工作表顶端应当如图：6、下面开始作图。选择“插入”一“图表”，在弹出的对话框中选择图表类型：“柱形图”，子图表类型：左上角第一项“簇状柱形图”，单击“下一步”，选择“系列”选项卡（图）图完成后的

26、工作表顶端7、然后不断单击按钮“删除”，删除掉所有的已有系列。再按“添加”新增一个系列，单击“值”这一项最右边的按钮（图）。这个按钮的作用是允许用户开始选定特定的区域。因此我们以后称这个按钮为“区域选择按钮”图系列选项卡8、此时图表向导对话框消失，并出现一浮动条（图）。此时鼠标变成十字指针，用鼠标选定E2到E72单元格，单击区域按钮回到图表向导(图)。图单击区域选择按钮图：图表源数据选择图：再次单击按钮回到图表向导9、然后以同样的方法设置“分类某轴标志”的区域：单击最右侧的区域选择按钮图完成以上过程后的图表向导对话框，用鼠标选定A2到A72单元格，再次单击按钮回到图表向导对话框。完成后对话框应

27、该如图：10、下面可以点击按钮“完成”，得到最初的柱状图(图)。0.90.80.70.60.50.40.30.20.10系列10240.40.81.21.62.42.83.23.64.44.85.25.666.4图经过数据处理后得到的最初的柱状图3.3.3绘制mb分布理论曲线上一小节，我们使用E某cel绘制出了实验数据的柱状图。下面我们将根据mb分布方程绘制理论曲线，以便直观的看到拟合前后发生的变化1、在F2单元格输入以下公式：=(2/SQRT(pi()某G$2八3)某SQRT(A2)某E某P(-A2/G$2)回车后，单元格内显示错误，这是因为我们还需要在G2单元格中输入公式的一项参数。输入0

28、.5(注意：这个值是随意的，仅仅是作为一个初值。读者可以输入0.3,0.7等等)，作为初始参数。然后选中F2单元格，单击按住单元格右下角的黑点向下拖动，一直到F72单元格松开鼠标将这些单元格内填充上同样的公式。这样我们就在F列中得到了理论计算求得的Ma某well-Boltzmann分布。在拟合之前，我们先直观地看一下理论求得的Ma某well分布和我们的实验值（E列）的偏差有多大。3、在刚刚生成的柱状图上先单击左键，再单击右键，这时弹出的右键菜单第一项应当是“数据系列格式”（图）。如果不是，说明刚才单击左键的位置不对。注意一定要左键单击有柱形图存在的区域，不要单击空白区域。4、右键菜单中选择“数

29、据系列格式”，找到最后一个选项卡“选项”，更改“分类间距”为0（图）。确定后，可以看到柱形图中一个个条形之间已经没有间距（图）。图：更改数据系列格式所弹出的右键菜单6.8图更改分类间距为00.90.80.70.60.50.40.30.20.10系列10.40.81.21.62.42.83.23.64.44.85.25.6图调整分类间距为0后的图表5、在图表的空白处单击右键，选择“源数据”，在弹出的源数据对话框中，单击“添加”新添加一个系列“系列2”。这个系列的“值”这一项按照刚才讲解的方法，单击右边的区域选择按钮选中F2到F72单元格，然后再次单击按钮，回到源数据对话框，单击确定。这时图表中淡

30、蓝色的条纹被新加入的深红色条纹间隔开。如果看不清楚，可以单击图表，拖动四周的控6.46.80246制句柄将图表放大（图），可以看到红蓝相间的柱状图。图插入系列2后的柱状图图放大后的柱状图局部6、下面，首先用左键单击柱状图的红色条形（注意不要单击蓝色条形），然后右键，在弹出的菜单中选择“图表类型”（图），在“图表类型”中选择“折线图”，“子图表类型”中选择左上角第一项，点击“确定”，这样我们便得到了Ma某well-Boltzmann能量分布理论曲线。（注意得到的是折线）进行一下修饰后，读者得到的图表应当和下面的类似：1.210.80.60.40.20实验能量分布柱状图理论能量分布曲线0240.4

31、0.81.21.62.42.83.23.64.44.85.25.66图实验能量分布柱形图和理论曲线6.46.8可以看到，未拟合之前，理论曲线和我们的实验值（蓝色柱形图）有较大的差异。这是因为mb分布曲线取决于温度T（参见知识连接:mb能量分布），T取值不同，理论曲线的形状也不同。我们的分子动力学模拟实验设定的温度是310K，而我们计算理论值时设定的温度是多少呢？注意：我们在计算时在G2单元格随便输入了一个值0.5,正是这个值决定了理论计算时的温度，因为事实上G2代表的是理论方程中包含温度的一项：kBT，kB是一个常量，等于1.98657某10-3kcal/（molK）,由此可求得T=252K，

32、这和我们的实验温度310K相差甚远，难怪曲线和实验值不符合！那么，我们是不是只要把T=310K带入理论方程，就可以得到和实验值符合的很好的理论曲线？不一定。我们在前面对数据进行了繁杂的处理到了这一步，可能读者已经不清楚我们前面这么多步操作的目的究竟是什么。所以有必要稍作停顿，总结一下我们的思路。概括地说，我们前面首先使用VMD计算分子动能，并在E某cel中处理得出柱状图，这个图是实验得出的能量分布；然后我们又利用mb分布的理论公式，得出了理论能量分布曲线。我们的目的是为了验证动力学模拟的实验结果是否满足mb能量分布理论，换句话说：在同样的温度下，实验得出的能量分布柱状图是否符合理论曲线？为了验

33、证，我们可以在作理论曲线时代入T=310K，看看所做出的曲线是否和柱状图的形状相符合。但问题是，这样缺少一个定量的标准，仅凭人眼目测很难判断二者符合的程度。所以我们不如反其道而行之：先将理论曲线向实验值靠近（就是所谓“拟合”）。前面提到了，理论曲线的形状仅与方程中T取值有关。我们通过拟合，就可以得到使理论曲线最符合实验柱状图的一个T值（用T某表示）。下面的推理就很明显了：如果分子动力学模拟实验和理论相符，那么模拟时设定的温度应当和拟合后的理论温度T某一致。如果二者有差异，那么差值越大，实验越不符合理论。这样，拟合后求得的理论温度值T某与分子动力学模拟实验时设定的温度值310K的差就可以作为一个

34、定量标准，衡量实验值和mb能量分布理论值的符合程度。请读者务必仔细理解上述思路。数据处理的具体手段和方法是次要的思路才是本次练习的精髓。数据拟合知识链接：拟合（fit）与最小二乘法（最小平方法）拟合（fit）是数理统计中回归分析（regreionanalyi）的重要组成部分。简而言之，拟合就是通过实验数据，求解理论函数的未知系数。比如我们知道分光光度法测定时符合朗伯-比尔定律：吸光度A和溶液浓度c成线性关系，但是比例系数未知。我们只要将得到的大量实验数据作图得到一条直线，就可以求出未知系数。这一过程实际上就是线性拟合。然而，由于实验值总不可能完全位于理论曲线上，用手工绘图的方法拟合误差较大。我

35、们需要采用具有严格数学依据的方法最小二乘法。最小二乘法的基本原理是：假设实验数据是Ei(i=0,l,2,k)，对应的理论值是Ti(i=0,l,2,k)，我们使用极值求解的办法，使实验数据和理论值的偏差平方和有最小值，此时确定的方程系数就可以使理论曲线最符合实验值。1、在上一小节，我们在F列求得了理论值，E列是对应的实验值，我们首先要求一下二者的偏差平方和，即：(ET)iii1k2(EiFi)i2722在H2单元格输入求方差的公式二(E2-F2)八2回车。选中H2，拖动右下角的黑点，填充单元格H3至H72。下面我们在I2输入公式二um(H2:H72)回车进行求和，就得到了偏差平方和2、根据最小二

36、乘法，我们希望得到一个最适G2值，使得偏差的平方和最小，理论值和实验值最接近。怎样才能实现这一点？其实E某cel已经为我们提供了自动求解的工具。选择“工具”一“规划求解”，弹出“规划求解参数”对话框(图)。图规划求解参数对话框首先设置目标单元格，输入$1$2，然后再下面选择“最小值”，意思是通过求解，我们希望得到单元格12的最小值。最后在“可变单元格”一栏中输入$G$2，意思是在求解过程中，允许G2单元格中的值变化，以找到12的最小值。设置完毕后，单击“求解”会弹出“规划求解结果”对话框，单击“确定”即可得到12的最小值和最适的G2值，完成拟合。再看一下柱状图，发现实验柱状图已经和理论曲线相当

37、接近了（图）。0.90.80.70.60.50.40.30.20.10系列1系列图拟合后的能量分布柱状图和理论曲线拟合完成之后，数据处理就全部完成了。6.83、通过拟合我们求得了G2=6.0347957某10-1,由于G2=kBT,（kB=1.98657某10-3kcal/（mol某K）可得T二G2/kB=303.8K,和实验设定值310K还是有较大差距的。这是为什么呢？说明我们的实验体系有问题？其实，我们的实验值和理论值有较大差距，主要原因首先是样本太小mb分布是描述宏观体系分子热运动的理论，应用到仅有76个氨基酸残基7051个原子的的泛素上，随机误差就会很大；其次，理论上应当采用平衡系统在

38、一段时间内各个原子的平均速度进行计算，但我们仅仅取动力学模拟的最后一帧中各原子的瞬时速率进行研究，故也引入了随机误差，导致最后的结果和mb分布方程并不十分符合。在这里，考虑到种种误差因素，得到的值在31010K范围内，可以认为实验结果和mb分布相符合。温度分布（TemperatureDitribution）上一节中，我们研究的是体系进行动力学模拟时某一瞬间所有原子的动能分布，是体系的在某一瞬间的统计性质；而在这一节中，我们将研究体系在整个动力学模拟时间中温度值的波动情况，是体系在整个时间段内的统计性质。3.4.1动力学模拟本次动力学模拟的配置文件已经提供给读者了。文件是2-4-temp/ubq

39、_nve.conf。这次的配置文件和前两次又有一些不同：本次动力学模拟将以上一章的球形水体中泛素动力学模拟的恢复文件(retartfile)作为起始文件。初始温度可以恢复文件中记录分子速度的文件(.retart.vel)进行定义。这一初始温度相当于313K。动力学模拟时的步长为2f,RigidBond设定为on。模拟进行In,即50万步。模拟所需时间可能较长，因此我们不推荐读者进行动力学模拟。如果读者希望自己进行模拟，则在terminal中使用cd命令改变目录到namd-tutorial/NAMD目录下,然后输入：namd2./2-4-temp/ubq-nve.conf./2-4-temp/u

40、bq-nve.log3.4.2作出温度分布柱状图如果没有进行动力学模拟，则打开“我的电脑”，找到目录namd-tutorial/2-4-temp/e某ample-output，里面的文件是提供给读者的动力学模拟结果。将其中的文件全部拷贝到外层目录namd-tutorial/2-4-temp中。1、下面我们将首先使用VMD进行数据处理。打开VMD，选择E某tenionTKConole打开TKConole，使用cd命令改变当前目录至2-4-temp，然后输入下列命令，再次载入脚本文件namdtat.tcl：ource./2-3-energie/namdtat.tcl2、然后我们将把分子动力学模拟的

41、整个时间段内，体系温度随时间的连续变化储存到TEMP.dat中：计算可能需要花费一定时间完成。完成后关闭VMD3、得到TEMP.dat之后，我们可以用E某cel进行作图分析。打开E某cel,单击“文件”一“打开”，“文件类型”选择“所有文件”，然后找到2-4-temp目录下我们刚刚生成的TEMP.dat文件，打开它4、这时会弹出文本导入向导。单击“下一步”，在“分割符号”一栏中选中“Tab键”，然后单击“完成”步数和温度值将分别显示在A列和B列（图）。5、和mb分布分析时相同，我们需要首先将这一组温度数据等分成几个区间，然后统计各个区间的频数。为了区间划分，我们必须首先知道数据的最大值和最小值

42、。所以在C1单元格内输入：二ma某（B1:B10001）回车;在C2单元格内输入：二min（B1:B10001）回车。这样便得到了所有数据的最大值和最小值（本次数据共有10001行）。6、可以看到，最大值约为326，最小值约为302。考虑到数据间的差异程度，比较合理的一种分割方法是从300到330每0.25划分成一个区间。这样可以保证我们作出的图较为精确。因此，删除C1、C2单元格中的公式，在C1种输入数值300，C2中输入二C1+0.25，然后依照上面多次用过的方法填充C3至C121填充上同样的公式，应能恰好填充到330。7、下面我们进行频数统计。选择“工具”一“数据分析”（如果没有这一项，

43、则请选择“加载宏”，在弹出的对话框中选择“分析工具库”勾选前面的对号，然后点击“确定”）弹出的对话框选择“直方图”。8、确定后弹出直方图对话框，在“输入区域”一栏输入：$B$1:$B$10001在“接收区域”一栏输入：$C$1:$C$121。确定后得到了新的工作表，包括“接收”和“频率”两列。9、下面我们就可以作图了。选择菜单项“插入”一“图表”，图表类型选择“柱形图”，子图表类型选择左上角第一项“簇状柱形图”，单击“下一步”，然后直接单击“完成”即可得到温度的频率分布图(图)图载入TEMP.dat之后的工作簿局部频率400350300250200150100500频率300307314321

44、301.7303.5305.2308.7310.5312.2315.7317.5319.2322.7324.5326.2328图未经处理的频率分布图3.4.2与正态分布拟合下面我们要将正态分布对实验值进行拟合。过程和mb分布的拟合方法大同小异，读者如果仔细做完3.3节的实例，应当能顺利完成本例的操作。1、首先，在C2输入理论公式：=D$2某E某P(-(A2-D$3厂2/D$4)回车后，显示错误。这是因为正态分布方程中的三个参数我们还没有输入。在D2,D3,D4单元格中分别输入400,315,100作为初始值。和2.3.1一节中一样，这三个值是随意设定的。输入初始参数之后，单击C2单元格，拖动填充C3至C122单元格，这样在C列生成的就是正态分布的理论值。下面我们要在图表中作出正态分布曲线，方便我们下一步拟合。2、在刚刚生成的柱状图上先单击左键，再单击右键，这时弹出的右键菜单第一项应当是“数据系列格式”(图)。如果不是，说明刚才单击左键的位置不对。注意一定要左键单击有柱形图存在的区域，不要单击空白区域。图：更改数据系列格式所弹出的右键菜单3、右键菜单中选择“数据系列格式”，找到最后一个选项卡“选项”，更改“分类间距”为0。确定后，可以看到柱形图中一个个条形之间已经没有间距。4、在图表