NAMD入门教程（二） - 图文

3. 分析方法

在实际工作中，分子动力学模拟体系的设计和模拟计算的进行可能并不困难。关键之处在于如何分析动力学模拟得出的结果，说明一定的问题。因此我们专门设置一章讲解动力学模拟的分析方法。

本章中我们将讲解四个例子：残基RMSD值，平均能量，mb分布和温度分布。每一个例子都很典型，代表了一定的分析思想。但是由于某些例子需要较长的计算时间，读者没有必要自己进行动力学模拟。在NAMD教程文件中已经给出了动力学模拟得到的结果。

3.1 每个氨基酸残基的RMSD值

上一章中，我们已经计算了整个蛋白质的RMSD值变化。我们初步提到了，RMSD值反映的是偏离平均位置的程度，是原子运动幅度的大小的衡量。这里我们将计算每个氨基酸残基的RMSD值，按照RMSD值对蛋白着色，并在最后进行问题讨论。本节的目的是使读者明白进行RMSD值分析的基本方法和意义所在。

3.1.1 RMSD值的计算与蛋白着色显示

我们将要使用一个VMD提供的脚本计算每个残基的平均RMSD值，并把这个值储存在VMD提供的用户存储区中（User Field），然后对蛋白进行着色。

1、首先，打开VMD，选择File→New Molecule菜单项，在弹出的文件浏览中找到 common目录下的文件usq_wb.psf，点击按钮Load。这样我们就载入了蛋白质的结构信息。（注意图形窗口中是空白的，并没有显示出蛋白结构，因为psf文件不存储各个原子的坐标）

2、此时Molecule File Browser窗口上端一栏应该显示Load files for: 0: ubq_wb.psf （图）说明如果再次载入新文件，所载入的文件将与usq_wb.psf结合起来一同处理显示。我们将载入原子运动轨迹文件，和psf结构文件配合即可显示出蛋白质原子的运动轨迹。因此点击Browse，找到1-2-sphere文件夹下的ubq_wb_eq.dcd，点击Load，这时关闭Molecule File Browser窗口。

图 载入ubq_wb.psf后的Molecule File Browser 3、打开VMD TKConsole，使用cd命令改变当前目录到namd-tutorial/2-1-rmsd/ 4、在这一目录下有我们准备好的脚本residue_rmsd.tcl（可能读者也注意到了，VMD中脚本文件的后缀名均为“.tcl”）。在TKConsole中输入：

source residue_rmsd.tcl

这一命令并不进行RMSD值的求算，而是载入了该脚本以供用户调用其中的命令。求算RMSD值所用的命令是：rmsd_residue_over_time。其格式为：

rmsd_residue_over_time 目标分子 所选残基

5、这里我们将选择蛋白质的所有氨基酸残基进行计算。选择方法是输入： set sel_resid [[atomselect top “protein and alpha”] get resid]

这一命令的作用是：定义一个变量sel_resid， 储存氨基酸残基的序号列表，方便下一步调用。在TKConsole中可以看到列出的氨基酸残基序号列表（图），由1至76依次编号。

6、对蛋白的所有残基（1号至76号）进行RMSD值计算。输入命令： rmsd_residue_over_time top $sel_resid

输入后VMD调用脚本计算每个残基的RMSD值。计算公式和2.6节中介绍的公式相同。此时应当能够在窗口中看到程序处理过程（图）：

待计算完成，输出结果应该如下图：

与此同时，程序会将每个原子的RMSD值存储在该原子的VMD用户存储区中。我们将使用VMD对每个氨基酸残基按照RMSD值着色，以清楚地显示各个残基在动力学模拟时的运动剧烈程度。

7、选择Graphics→Representations菜单项。

8、在Atom Selection window一项中输入protein后回车，然后在Drawing Method下拉列表中选择Tube，在Coloring Method下拉列表中选择User。然后单击Trajectory 选项卡，在Color Scale Data Range中输入0.40 和1.00，单击Set。

以上设置方法如图：

现在在图形窗口中，应该可以看到按照残基RMSD值着色的蛋白骨架（图）。

在这幅彩图中，蓝色的氨基酸残基是运动幅度最大的，红色残基是幅度最小的，这和我们通常的着色方式恰好相反。我们要进行一下更改：

9、选择Graphics→Colors菜单项。选择Color Scale 选项卡，在Method下拉列表中选择BWR。现在红色是活动剧烈的残基，蓝色是运动幅度较小的残基，白色介于二者之间（图）。转动分子仔细观察，氨基酸残基运动的幅度和二级结构有什么关系？（如果不方便观察二级结构，可以选择Graphics→Representation，将Drawing Method改变为Cartoon。

3.1.2 氨基酸RMSD值分布图

我们下面将使用Excel处理数据，作出RMSD值分布图

10、关闭VMD，打开Excel，选择“文件”→“打开”，找到2-1-rmsd目录，选择文件residue_rems.dat打开，注意“文件格式”下拉列表要选择“所有文件”，否则将看不到我们所要的文件。

11、在弹出的“文件导入向导”中，首先单击“下一步”，然后单击“完成”。可以看到A列是氨基酸残基的序号（1至76），B列是每个残基对应的RMSD值。

12、选择“插入”→“图表”，在“图表类型”中选择“XY 散点图”，在“子图表类型”中选择左下角的“折线散点图”，点击“完成”。

得到的图表（图）就是泛素的氨基酸残基RMSD值变化分布。

1.61.41.210.80.60.40.20020406080系列1

当然，必须记住我们不是为了作图而作图。我们作图的目的是分析RMSD值的变化，发现问题和规律。仔细观察你得到的RMSD值图，可以看到一些有意思的问题。

以下是我们提出的两个问题，供读者思考：

1、我们将蛋白质的二级结构大体分成三类：α螺旋，β片层和无规则卷曲（loop）。使用前面讲过的DeepView可以得到泛素二级结构的氨基酸残基分布：

1-7 片层；11-17 片层；23-33 螺旋；40-45 片层；48-50 片层；57-59 螺旋；65-72 片层，其余部分认为是无规则卷曲（loop）。

在Excel中计算一下哪一种二级结构的平均RMSD值较高？哪一种二级结构的平均RMSD值较低？注：二级结构的平均RMSD值 = 该种二级结构的总RMSD值 / 该二级结构氨基酸残基数。思考一下这种规律背后的原因是什么？

参考答案： β片层：0.537； α螺旋：0.562； loop：0.744

2、从我们作出的图中可以看到，C末端最后4个氨基酸的RMSD值明显较高。 在VMD中观察C末端loop的构象，有什么特色？

前面提到，泛素的功能是通过标记特定蛋白质，使得被标记的蛋白特异性降解（参见知识链接：泛素——死亡之吻 ）。你认为泛素最有可能是通过哪一端和被降解蛋白相连？C末端还是N末端？如果已知泛素和标记蛋白的半胱氨酸相连，推测一下泛素末端的哪一个氨基酸负责连接？

参考答案：C末端；Arg

思考过上面两个问题后，可能读者已经初步理解了RMSD值分析的意义。我们从X射线衍射得到的蛋白结构仅仅是一个静态的“死”模型，这一结构并不能直接告诉我们各个氨基酸残基的运动特性。但是进行动力学模拟时，我们手中的蛋白质模型就由“死”变“活”，运动起来了。而通过RMSD值分析，我们可以初步了解各个氨基酸残基的空间位阻大小，运动的方式，以及不同的二级结构的运动特点等等信息。这些信息往往同蛋白质的功能相联系（比如上面第二个问题）。总而言之，RMSD值分析可以为我们阐明结构和功能的关系提供思路和突破口。

3.2 能量分析

上一节中我们所分析的RMSD值反映的是体系中不同原子的运动情况。在这一节中，我们将对体系能量进行分析。下面我们将计算整个体系中不同形式的能量（如动能，键能，静电势能等）在动力学模拟时的平均值。

我们将计算上一章完成的立方水体分子动力学模拟过程中的能量平均值。计算时我们又一次需要使用一个脚本文件：namdstats.tcl。这个脚本文件在进行能量平均值的计算时是十分有用的。

1、打开VMD，打开TKConsole，使用cd命令改换当前目录至 2-3-energies/ 2、在TKConsole中输入： source namdstats.tcl

data_avg ../1-3-box/ubq_wb_eq.log 101 last

第一行命令将载入脚本文件namdstats.tcl, 第二行命令将提取日志文件中记录能量的那一部分并进行计算。101 last表示计算动力学模拟第101步至最后一步过程中的平均能量值。

3、结果会直接输出在TKConsole中（图）。

图 计算能量平均值的输出结果 上图中所使用的缩写：BOND键能，ANGLE键角能 DIHED二面角能 ELECT静电能 VDW范德华力能 KINETIC动能 TEMP温度。在这里仅需要读者明白这些平均能量计算的方法，熟悉namdstats.tcl脚本的使用即可。

本节到此就算完成了。但是关于脚本文件namdstats.tcl我们还需要多了解一些，因为这一脚本在进行动力学模拟结果分析时是非常常用的。

脚本namdstats.tcl定义了两个命令：data_avg和data_time。

? data_avg 的功能是计算特定时间段内整个体系各种能量的平均值。调用方法是： data_avg <日志文件> <起始> <终止>

其中，<日志文件>表示用于计算的”.log” 文件的位置。<起始> 和<终止>分别代表计算。第一个timestep 可以用first代替，最后一个timestep可以用last代替。因此我们在例子中输入的101 last，表示的是从第101步计算到最后一个timestep。如果省略<起始>和<终止>，仅仅输入data_avg <日志文件>，VMD将默认对于动力学模拟的整个过程进行计算。

? data_time 的功能是输出特定时间段内整个体系某个特定能量随时间的连续变

化。调用方法是：

data_time <能量> <日志文件> <起始> <终止>

NAMD将会从日志文件中提取<能量>这一参数设定的能量值从起始timestep到终止timestep的连续变化过程，并输出到<能量>.dat这一文件中。<起始>和<终止>的设定方法同上。

在动力学模拟的多种分析方法中，能量分析是最基础最经典的，通过能量分析可以直接或间接说明许多问题。 但是进行分析时需要过硬的理论功底，因此本书仅作简单介绍，其他具体的分析实例请参照相关的文献专著。

上两节我们分别对体系中原子的运动和体系的能量进行了简单分析。我们的分析有助于阐明泛素的结构和功能的关系，以及泛素自身的某些性质，因此我们得出的结果具有明确的生物学意义。

在下面两节中，我们将抛开我们所研究体系的生物学意义，对系统的统计力学性质进行分析。读者可能会问：我们做的是生物学研究，为什么不分析体系的生物学性质，却要考察体系的纯物理学性质？这里有必要作一下说明。

我们知道，物理学是公理化的科学，经过几百年的发展，理论体系相对较为完备，纯粹理论推导出的结果能够很好地与实验符合。但生命科学依然是描述性的科学。得出的许多基本规律，比如中心法则，遗传密码，其本质是经验性的，没有理论推导可言。而在分子动力学模拟领域，虽然力学规律是严格的微分方程，基本的相互作用力场参数（如CHARMM，Amber等）却是半经验式的，依靠大量实验得出的，仅在统计意义上成立，而不能保证动力学模拟的结果一定正确。

除了力场参数引入的误差，我们在进行模拟时还必须人为地做出一些假设。读者可以再仔细阅读一下2.3.1一节，配置文件中我们的设定：我们假定范德华力和静电力只在某一范围内起作用，假定计算时某些原子的相互作用被忽略等等，都使得体系越来越不精确。而在真正的研究工作开展时，我们可能需要在体系中放入离子以模拟真实的细胞质基质，需要使用格点计算静电力场分布??这些人为引入的假设和近似能否符合生命体系中的真实过程？从基础理论的推导我们得不出答案。

那么，如何验证我们得出的结果的正确性？一个方法当然是通过实际实验检验，但实验可能需要大量时间，并且要受条件的限制。因此，另一个方法是：使用物理学基本规律进行证伪。前面说到了，物理学规律是建立在严密的理论基础之上的，其结论是经历过实验的考验的。因此我们对所得结果进行分析，如果符合物理学基本规律，并不能证明结果正确——因为还不能说明是否符合生命体系中的真实过程；但是如果不符合，那么一定能说明结果是错误的。因此，物理学基本规律可以作为一个只能证伪的有力判据。

具体到本例，我们将检验分子动力学模拟得出的结果是否符合统计物理中的两个结论：mb分布和温度分布。如果不符合，说明我们的结果不正确，原因可能是体系有问题，力场参数不适合，边界条件设置不当等等。

同时，完成下面两节后，读者应该掌握一些重要的数据处理思想和方法。

3.3 麦克斯韦-波尔兹曼（Maxwell-Boltzmann ）能量分布

下面我们将通过对分子动力学模拟所得结果的分析，验证动力学模拟的结果是否符合Maxwell-Boltzmann能量分布。本节内容较多，也较为综合，并将详细讲述Excel回归分析，拟合数据的基本方法。

知识链接：麦克斯韦-波尔兹曼（Maxwell-Boltzmann ）能量分布 1、诞生背景 在18世纪后期，分子运动论的发展使人们认识到：气体是由大量作无规则热运动的分子组成的。物理学家自然希望能够对大量分子的运动使用牛顿力学进行分析。但是对于大量微粒组成的体系，用牛顿力学的微分方程去逐一推算显然是不现实的。 为了解决这一难题，19世纪的物理学巨匠麦克斯韦首先将数理统计原理引入对分子运动的研究，这样一来，极大的分子数目不但不再是研究的阻碍，反而正好可以使统计平均值有效。麦克斯韦得出了著名的麦克斯韦分子速度分布率。后来，奥地利物理学家波尔兹曼发展了这一理论，得出了波尔兹曼分布。通常我们称为麦克斯韦-波尔兹曼分布。 2、基本内容 概括地说，麦克斯韦-波尔兹曼能量分布（Maxwell-Boltzmann distribution）所描述的是近独立体系（比如一团理想气体）达到平衡态之后，气体分子的能量分布情况，举个通俗的例子，对于给定温度的一团理想气体，动能是100J的分子有多少个? 是200J的分子有多少个？mb分布就是定量描述能量分布的函数。考虑到篇幅限制，这里不经推导直接给出

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

计算动能分布自然首先要得到各个原子的动能，而动能的计算需要得到某一瞬间各个原子的速率。在2.6节我们提到，NAMD动力学模拟的输出文件包括“.vel”文件，该文件记录原子的瞬时速率。我们这里使用的就是立方水体中泛素分子动力学模拟时原子的速率文件。

虽然我们在2.5节已经进行了立方水体的动力学模拟，但是在模拟时我们设定了Rigid Bonds（见2.3.1节中知识链接：Rigid Bonds ）,步长是2fs，这样的近似设定使得体系不能精确反映各个原子的运动速率。因此我们取消Rigid Bonds，步长1fs重新进行动力学模拟，并将结果输出文件提供给读者。

1、首先，我们仍旧要载入泛素的结构信息。打开VMD，选择File→New Molecule菜单项，找到common目录下的文件ubq_wb.psf，载入该文件。不要关闭Molecule File Browser窗口，窗口第一项应该显示Load file for: 0:ubq_wb.psf。

2、下面载入泛素中各原子的分子速率信息。再次单击Browser按钮，找到2-2-maxwell目录下，载入ubq_wb_eq_1fs.restart.vel文件。在Determine file type：下拉列表中选择NAMD Binary Coordinates，然后再点击Load。现在可以关闭Molecule File Browser。

图 Molecule File Browser 载入速率文件时的设置 可以在图形窗口中看到一个海胆一样浑身是刺的怪物！这是因为VMD把各个原子的速度当作原子的坐标处理（瞬时速度是以x,y,z三个方向的分速度记录的，记录格式和座标一样）。不过这并不影响我们下一步的操作。我们将使用VMD得到只存储各个原子质量和瞬时速度的文件，便于我们计算动能。

3、在VMD TKConsole窗口中，用cd命令改变目录到namd-tutorial/2-2-maxwell。

图 载入速率文件后图形窗口的显示 4、然后输入：

set all [atomselect top all]

这样我们建立了一个变量all，储存系统中的所有原子。（系统返回：atomselect0，表示这一选择的序号是0）

5、下面在当前目录下新建一个文件energy.dat : set fil [open energy.dat w] 6、下面计算每个原子的动能：?k?12mv,并把计算结果写入文件energy.dat。 2foreach m [$all get mass] v [$all get {x y z}] {

puts $fil [expr 0.5 * $m * [vecdot $v $v] ] }

(注意：连续输入既可，不需要每一行回车) 7、关闭文件： close $fil 然后退出VMD

这样我们就获得了文件energy.dat，存储所有的原子及其相应的动能。

3.3.2 绘制实验数据柱状图

下面我们将首先使用Excel处理原始数据，做出柱状图，作为我们下一步拟合的准备工作。

1、打开Excel，选择“文件”→“打开”，找到2-2-maxwell 目录下的energy.dat 文件并打开它，注意“文件类型”设置为“所有文件”

2、此后会弹出文本导入向导，直接点击“完成”，这样在A列显示的就是各个原子的动能值。

下面我们将把这些动能值处理成柱状图，并和Maxwell-Boltzmann分布拟合，最后可以得出系统的平均温度。

1、我们首先要做的是：统计所有原子的动能值的分布频数。打个比喻，我们现在得到了一个班全体同学身高的列表，身高的分布频数就是统计有多少人的身高在[150cm, 155cm], [155cm,165cm],[165cm, 170cm]...区间之内。显然，为了将动能分布等分成区间，我们首先要知道这些动能值中的最大值和最小值。

在B1中输入 =max（A1：A7051）回车。max是Excel提供的求指定范围内最大值的函数。我们指定的范围是A1到A7051，因为向下拖动滚动条可知一共有7051个原子的动能值。同样地，在B2中输入 =min（A1：A7051）回车，可得到这一范围内的最小值。结果如图

因为最大值是6.5左右，最小值在0附近，我们将0到7以0.1为间距等分成70个区间，可以包含所有的动能值。

将B1和B2单元格中的公式删除，然后在B1中输入0，B2中输入 =B1+0.1 回车，然后再次选中B2，单击B2右下角的小黑点后，不要松开鼠标左键，一直向下拖动到B71后松开左键。可以看到单元格已经被自动填充上了我们想要的区间（图）

图 利用max和min函数得到最大最小值 2、下面我们将使用专门的分析工具统计动能值在这70个区间中的分布频数 单击“工具”→“加载宏”，在弹出的对话框中选择“分析工具库”，点击“确定”。再次打开“工具”菜单，应该出现“数据分析”这一项（图）。

图 得到分布区间

图 “加载宏”对话框

单击打开“数据分析”。选择“直方图”，在弹出的“直方图”对话框中填入：

输入数据：$A$1:$A$7051 接受区域：$B$1:$B$71

单击“确定”后会自动新建一个新的工作表，A列为“接受”，B列为“频率”。我们可以保存一下以防工作丢失。不要单击

，选择“文件”→“另存为”，“保存类型”选择第一项“Microsoft

图 “直方图” 对话框 图 “工具”中新增的“数据分析”菜单项 Office Excel 工作簿”，重新命名为rmsd，单击“保存”。

3、保存之后，我们开始对数据进行作图。我们希望获得的是标准化的频率分布图。所谓标准化（normalization）就是使得频率分布柱状图的面积之和为一。为了达到这一点，需要先对数据进行如下处理：

在C2单元格中输入公式 =B2*0.1 回车，然后按照上面用过的方法单击并拖动C2单元格的右下角一直到第72行，将C2到C72填充该公式。

4、在D2单元格输入公式 =sum（C2：C72）回车。sum是Excel提供的求和函数。这样在D2中求得了C2至C72之中所有数据的和。

5、在E2单元格中输入公式 =B2/D$2回车。同样地，单击并拖动E2的右下角一直到第72行。这样数据处理过程就进行完毕了。完成后的工作表顶端应当如图：

6、下面开始作图。选择“插入”→“图表”，在弹出的对话框中选择图表类型：“柱形图”，子图表类型：左上角第一项“簇状柱形图”，单击“下一步”，选择“系列”选项卡（图）

图 完成后的工作表顶端 7、然后不断单击按钮“删除”，删除掉所有的已有系列。再按“添加”新增一个系列，单击“值”这一项最右边的按钮（图）。这个按钮的作用是允许用户开始选定特定的区域。因此我们以后称这个按钮

为“区域选择按钮”

图 系列选项卡 8、此时图表向导对话框消失，并出现一浮动条（图）。此时鼠标变成十字指针，用鼠标选定E2到E72单元格，单击区域按

钮回到图表向导（图）。

图 单击区域选择按钮 图：图表源数据选择 图：再次单击按钮回到图表向导

9、然后以同样的方法设置“分类X轴标志”的区域：单击最右侧的区域选择按钮

图 完成以上过程后的图表向导对话框 ，

用鼠标选定A2到A72单元格，再次单击按钮回到图表向导对话框。完成后对话框应该如图：

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分布理论曲线

上一小节，我们使用Excel绘制出了实验数据的柱状图。下面我们将根据mb分布方程绘制理论曲线，以便直观的看到拟合前后发生的变化。

1、在F2单元格输入以下公式：

= (2 / SQRT(pi() * G$2^3)) * SQRT(A2) * EXP(-A2/G$2)

回车后，单元格内显示错误，这是因为我们还需要在G2单元格中输入公式的一项参数。输入0.5（注意：这个值是随意的，仅仅是作为一个初值。读者可以输入0.3，0.7等等），作为初始参数。然后选中F2单元格，单击按住单元格右下角的黑点向下拖动，一直到F72单元格松开鼠标，将这些单元格内填充上同样的公式。这样我们就在F列中得到了理论计算求得的Maxwell-Boltzmann分布。

在拟合之前，我们先直观地看一下理论求得的Maxwell分布和我们的实验值（E列）的偏差有多大。

3、在刚刚生成的柱状图上先单击左键，再单击右键，这时弹出的右键菜单第一项应当是“数据系列格式”（图）。如果不是，说明刚才单击左键的位置不对。注意一定要左键单

击有柱形图存在的区域，不要单击空白区域。

4、右键菜单中选择“数据系列格式”，找到最后一个选项卡“选项”，更改“分类间距”为0（图）。确定后，可以看到柱形图中一个个条形之间已经没有间距（图）。 图：更改数据系列格 式所弹出的右键菜单

6.8

图 更改分类间距为0 0.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单元格，然后再次单击按钮

，

，回到源数据对话框，单击确定。这时图表中

淡蓝色的条纹被新加入的深红色条纹间隔开。如果看不清楚，可以单击图表，拖动四周的控

6.46.80246制句柄将图表放大（图），可以看到红蓝相间的柱状图。

图 插入系列2后的柱状图 图 放大后的柱状图局部 6、下面，首先用左键单击柱状图的红色条形（注意不要单击蓝色条形），然后右键，在弹出的菜单中选择“图表类型”（图），在“图表类型”中选择“折线图”，“子图表类型”中选择左上角第一项，点击“确定”，这样我们便得到了Maxwell-Boltzmann能量分布理论曲线。（注意得到的是折线）进行一下修饰后，读者得到的图表应当和下面的类似：

1.210.80.60.40.20实验能量分布柱状图理论能量分布曲线0240.40.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.98657x10-3kcal/( mol· K),由此可求得T=252K，这和我们的实验温度310K相差甚远，难怪曲线和实验值不符合！

那么，我们是不是只要把T=310K带入理论方程，就可以得到和实验值符合的很好的理论曲线？不一定。我们在前面对数据进行了繁杂的处理，到了这一步，可能读者已经不清楚我们前面这么多步操作的目的究竟是什么。所以有必要稍作停顿，总结一下我们的思路。

概括地说，我们前面首先使用VMD计算分子动能，并在Excel中处理得出柱状图，这个图是实验得出的能量分布；然后我们又利用mb分布的理论公式，得出了理论能量分布曲线。

我们的目的是为了验证动力学模拟的实验结果是否满足mb能量分布理论，换句话说：在同样的温度下，实验得出的能量分布柱状图是否符合理论曲线？为了验证，我们可以在作理论曲线时代入T=310K，看看所做出的曲线是否和柱状图的形状相符合。但问题是，这样缺少一个定量的标准，仅凭人眼目测很难判断二者符合的程度。

所以我们不如反其道而行之：先将理论曲线向实验值靠近（就是所谓“拟合”）。前面提到了，理论曲线的形状仅与方程中T取值有关。我们通过拟合，就可以得到使理论曲线最符合实验柱状图的一个T值（用T*表示）。下面的推理就很明显了：如果分子动力学模拟

实验和理论相符，那么模拟时设定的温度应当和拟合后的理论温度T* 一致。如果二者有差异，那么差值越大，实验越不符合理论。这样，拟合后求得的理论温度值T* 与分子动力学模拟实验时设定的温度值310K的差就可以作为一个定量标准，衡量实验值和mb能量分布理论值的符合程度。

请读者务必仔细理解上述思路。数据处理的具体手段和方法是次要的，思路才是本次练习的精髓。

3.3.4 数据拟合

知识链接：拟合（fit）与最小二乘法（最小平方法） 拟合（fit）是数理统计中回归分析（regression analysis）的重要组成部分。简而言之，拟合就是通过实验数据，求解理论函数的未知系数。比如我们知道分光光度法测定时符合朗伯-比尔定律：吸光度A和溶液浓度c成线性关系，但是比例系数ε未知。我们只要将得到的大量实验数据作图得到一条直线，就可以求出未知系数ε。这一过程实际上就是线性拟合。 然而，由于实验值总不可能完全位于理论曲线上，用手工绘图的方法拟合误差较大。我们需要采用具有严格数学依据的方法——最小二乘法。 最小二乘法的基本原理是：假设实验数据是Ei（i=0,1,2,…,k），对应的理论值是Ti (i=0,1,2,…,k)，我们使用极值求解的办法，使实验数据和理论值的偏差平方和有最小值，此时确定的方程系数就可以使理论曲线最符合实验值。 1、在上一小节，我们在F列求得了理论值，E列是对应的实验值，我们首先要求一下二者的偏差平方和，即：

?(E?T)iii?1k2?(Ei?Fi)i?2722

在H2单元格输入求方差的公式 =（E2-F2）^2回车。选中H2，拖动右下角的黑点，填充单元格H3至H72。下面我们在I2输入公式 =sum(H2:H72)回车进行求和，就得到了偏差平方和

2、根据最小二乘法，我们希望得到一个最适G2值，使得偏差的平方和最小，理论值和实验值最接近。怎样才能实现这一点？其实Excel已经为我们提供了自动求解的工具。选择“工具”→“规划求解”，弹出“规划求解参数”对话框（图）。

图 规划求解参数对话框 首先设置目标单元格，输入$I$2，然后再下面选择“最小值”，意思是通过求解，我们希望得到单元格I2的最小值。最后在“可变单元格”一栏中输入$G$2，意思是在求解过程中，允许G2单元格中的值变化，以找到I2的最小值。设置完毕后，单击“求解”会弹出“规划求解结果”对话框，单击“确定”即可得到I2的最小值和最适的G2值，完成拟合。

再看一下柱状图，发现实验柱状图已经和理论曲线相当接近了（图）。

0.90.80.70.60.50.40.30.20.10系列1系列20240.40.81.21.62.42.83.23.64.44.85.25.666.4图 拟合后的能量分布柱状图和理论曲线

拟合完成之后，数据处理就全部完成了。

6.8 3、通过拟合我们求得了G2 = 6.0347957 x 10-1, 由于G2 = kBT，（kB

=1.98657x10-3kcal/( mol*K) ）可得 T = G2 / kB = 303.8K, 和实验设定值310K还是有较大差距的。这是为什么呢？说明我们的实验体系有问题？

其实，我们的实验值和理论值有较大差距，主要原因首先是样本太小：mb分布是描述宏观体系分子热运动的理论，应用到仅有76个氨基酸残基，7051个原子的的泛素上，随机误差就会很大；其次，理论上应当采用平衡系统在一段时间内各个原子的平均速度进行计算，但我们仅仅取动力学模拟的最后一帧中各原子的瞬时速率进行研究，故也引入了随机误差，导致最后的结果和mb分布方程并不十分符合。在这里，考虑到种种误差因素，得到的值在310?10K范围内，可以认为实验结果和mb分布相符合。

3.4 温度分布(Temperature Distribution)

上一节中，我们研究的是体系进行动力学模拟时某一瞬间所有原子的动能分布，是体系的在某一瞬间的统计性质；而在这一节中，我们将研究体系在整个动力学模拟时间中温度值的波动情况，是体系在整个时间段内的统计性质。

知识链接：温度分布 使用写字板打开common目录下的ubq_ws_eq.log文件，这是我们第一次模拟所得到的输出文件。 向下拖动右侧的滚动条一直到文件最后，向右拖动下方的滚动条一直到文件最右，图中所标注出的体系的温度。可以看到，该温度并不是我们所设定的310K，而在310K上下波动。这是为什么呢？ 我们知道，对于这样一个只有7000多个原子的微观体系，温度就是体系中各原子的平均动能，事实上NAMD在输出温度时就是这样计算的。很显然，由于体系中动能-势能是不断相互转化的，体系平均动能不可能是一个定值，而会在310K上下波动。理论上来说波动是一个随机过程，应该符合正态分布。我们下面所要做的就是验证动力学模拟过程中，温度这一随机变量的波动是否符合正态分布。

3.4.1 动力学模拟

本次动力学模拟的配置文件已经提供给读者了。文件是2-4-temp /ubq_nve.conf。这次的配置文件和前两次又有一些不同：

? 本次动力学模拟将以上一章的球形水体中泛素动力学模拟的恢复文件（restart

file）作为起始文件。

? 初始温度可以恢复文件中记录分子速度的文件（.restart.vel）进行定义。这一初

始温度相当于313K。

? 动力学模拟时的步长为2fs，Rigid Bonds 设定为on。模拟进行1ns，即50万步。 模拟所需时间可能较长，因此我们不推荐读者进行动力学模拟。如果读者希望自己进行模拟，则在terminal中使用cd命令改变目录到namd-tutorial /NAMD目录下，然后输入：

namd2 ../2-4-temp/ubq-nve.conf > ../2-4-temp/ubq-nve.log 3.4.2 作出温度分布柱状图

如果没有进行动力学模拟，则打开“我的电脑”，找到目录namd-tutorial /2-4-temp/example-output，里面的文件是提供给读者的动力学模拟结果。将其中的文件全部拷贝到外层目录 namd-tutorial /2-4-temp中。

1、下面我们将首先使用VMD进行数据处理。打开VMD，选择Extension→TKConsole打开TKConsole，使用cd 命令改变当前目录至 2-4-temp，然后输入下列命令，再次载入脚本文件 namdstats.tcl：

source ../2-3-energies/namdstats.tcl

2、然后我们将把分子动力学模拟的整个时间段内，体系温度随时间的连续变化储存到TEMP.dat中：

data_time TEMP ubq-nve.log

计算可能需要花费一定时间完成。完成后关闭VMD

3、得到TEMP.dat之后，我们可以用Excel进行作图分析。打开Excel，单击“文件”→“打开”，“文件类型”选择“所有文件”，然后找到2-4-temp目录下我们刚刚生成的TEMP.dat文件，打开它

4、这时会弹出文本导入向导。单击“下一步”，在“分割符号”一栏中选中“Tab 键”，然后单击“完成”步数和温度值将分别显示在A列和B列（图）。

5、和mb分布分析时相同，我们需要首先将这一组温度数据等分成几个区间，然后统计各个区间的频数。为了区间划分，我们必须首先知道数据的最大值和最小值。所以在C1单元格内输入：=max(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、下面我们进行频数统计。选择“工具”→“数据分析”（如果没有这一项，则请选择“加载宏”，在弹出的对话框中选择“分析工具库”，勾选前面的对号，然后点击“确定”）弹出的对话框选择“直方图”。

8、确定后弹出直方图对话框，在“输入区域”一栏输入：$B$1:$B$10001 在“接收区域”一栏输入：$C$1: $C$121。确定后得到了新的工作表，包括“接收”和“频率”两列。

9、下面我们就可以作图了。选择菜单项“插入”→“图表”，图表类型选择“柱形图”，子图表类型选择左上角第一项“簇状柱形图”，单击“下一步”，然后直接单击“完成”即可得到温度的频率分布图（图）。

图 载入TEMP.dat 之后的工作簿局部 频率400350300250200150100500频率300307314321301.7303.5305.2308.7310.5312.2315.7317.5319.2322.7324.5326.2328 图 未经处理的频率分布图 3.4.2 与正态分布拟合

下面我们要将正态分布对实验值进行拟合。过程和mb分布的拟合方法大同小异，读者如果仔细做完3.3节的实例，应当能顺利完成本例的操作。

1、首先，在C2输入理论公式： =D$2 * EXP(-((A2-D$3)^2/D$4)

回车后，显示错误。这是因为正态分布方程中的三个参数我们还没有输入。在D2，D3，D4单元格中分别输入400，315，100作为初始值。和2.3.1一节中一样，这三个值是随意设定的。

输入初始参数之后，单击C2单元格，拖动填充C3至C122单元格，这样在C列生成的就是正态分布的理论值。下面我们要在图表中作出正态分布曲线，方便我们下一步拟合。

2、在刚刚生成的柱状图上先单击左键，再单击右键，这时弹出的右键菜单第一项应当是“数据系列格式”（图）。如果不是，说明刚才单击左键的位置不对。注意一定要左键单击有柱形图存在的区域，不要单击空白区域。 图：更改数据系列格 式所弹出的右键菜单 3、右键菜单中选择“数据系列格式”，找到最后一个选项卡“选项”，更改“分类间距”为0。确定后，可以看到柱形图中一个个条形之间已经没有间距。

4、在图表的空白处单击右键，选择“源数据”，在弹出的源

数据对话框中，单击“添加”新添加一个系列“系列2”。这个系列的“值”这一项按照刚才讲解的方法，单击右边的区域选择按钮

，选中C2到C122单元格，然后单击按钮

，回到源

329.7

数据对话框，单击确定。这时的图表如下图：

450400350300250200150100500频率系列2300307314321301.7303.5305.2308.7310.5312.2315.7317.5319.2322.7324.5326.2328 329.7 6、新加入的“系列2”就是理论正态分布。和2.3.1例中一样，我们希望得到理论值得分布曲线而不是分布柱形图，因此需要更改图表类型。首先左键单击系列2柱形图部分（图），然后点右键，在弹出菜单中选择“图表类型”（图），选择“折线图”中的左上角第一项。单击“确定”，系列2就以紫红色曲线表示出来。稍稍修饰后，可以得到如下图表：

450400350300250200150100500实验温度分布理论温度分布300305310315320325302.5307.5312.5317.5322.5301.25303.75306.25308.75311.25313.75316.25318.75321.25323.75326.25327.5图 实验温度分布柱形图和理论温度分布曲线 7、可以看到，未进行拟合之前，理论曲线和实验值有很大差距。下面我们将用3.2.1例中用过的方法，通过最小二乘法拟合。在E2单元格中输入公式=（C2-B2）^2回车，然后再次选中E2，拖动右下角填充E3至E122。在F2单元格中输入公式 =sum(E2:E122)回车，F2中的值就是误差平方和。

8、我们仍旧使用规划求解去求误差平方和的最小值。选择“工具”→“规划求解”，然后如下设置：

设置目标单元格：$F$2 等于： 最小值

可变单元格：$D$2:$D$4

单击“求解”，Excel将变动D2、D3、D4三个参数，求解目标单元格的最小值，此时曲线已经和实验值拟合得非常好（图）。

328.75330 400350300250200150100500实验温度分布理论温度分布300305310315320325301.25303.75306.25308.75311.25313.75316.25318.75321.25323.75326.25 图 拟合后的温度分布柱状图和理论曲线 事实上，我们刚才在D2，D3，D4三个单元格中填入的三个值分别是归一化常数（normalization constant），体系温度T（即正态分布中的平均值?），以及分布的方差?。因此我们也求得了体系的正态分布理论平均温度T：313.92K。我们下面使用namdstats.tcl脚本中data_avg命令求解体系的实验平均温度，看一看二者是否符合？

打开VMD，进入TKConsole，首先用cd命令改变当前目录至2-4-temp，然后输入： source ../2-3-energies/namdstats.tcl 载入脚本之后，输入： data_avg ubq-nve.log

耐心等待计算完成后，得到实际平均温度313.75K，可见理论值和实际值还是相当符合的。另外，如果不想使用VMD完成，也可以在Excel中使用average函数完成。

完成上面两个例子后，我们就完成了对于分子动力学模拟实验结果的证伪过程。上面两个例子说明我们的实验结果较好地符合了物理学基本规律。但是实验结果能否符合生命过程的真实情况，还需要通过真实的实验进行检验。

2

328.75302.5307.5312.5317.5322.5327.5330

本文来源：https://www.bwwdw.com/article/t08t.html