预定目录 1. 分子动力学模拟概论 1.1 分子动力学模拟的发展 1.2 分子动力学模拟的基本原理 1.3 分子动力学模拟相关软件
2. 分子动力学入门 2.1 基本设置
2.2 生成蛋白质结构文件（PSF） 2.3 蛋白质的溶质化
2.4 球状水体中泛素（Ubiquitin）的分子动力学模拟 2.5 立方水体中泛素（Ubiquitin）的分子动力学模拟 2.6 简单的结果分析
3. 分析方法
3.1 平衡态分子动力学模拟分析 3.1.1 每个残基的RMSD值
3.1.2 麦克斯韦-波尔兹曼（Maxwell-Boltzmann ）能量分布 3.1.3 能量分析 3.1.4 温度分布 3.1.5 比热分析
3.2 非平衡态分子动力学模拟分析 3.2.1 热扩散 3.2.2 温度回音
4 人工操纵的分子动力学模拟（SMD） 4.1 除去水分子 4.2 恒速拉伸 4.3 恒力拉伸 4.4 结果分析
1. 分子动力学模拟概论
分子动力学模拟（Molecular Dynamics Simulation）是指利用计算机软件，根据牛顿力学的基本原理，模拟大分子的相互作用和运动变化的研究方法。生命科学的研究往往离不开各种仪器，试管和活的有机体，通过实验手段研究生命现象背后的规律。那么，为什么我们要将生命大分子抽象成二进制数据，由计算机软件模拟其行为呢？
首先，从理论基础上讲，我们能够使用计算机模拟生物大分子的行为。生物体系非常复杂，但生物大分子如蛋白质，脂肪，多糖等也是许多原子由化学键连接起来形成的，所有原子的运动规律都符合量子力学方程，在较大尺度上也近似符合牛顿力学方程，它的行为是要受物理学基本规律支配的。因此我们可以将利用纯数学的手段，近似模拟生物大分子的行为.
其次，从研究需要上讲，我们不仅希望从宏观上研究生命大分子溶液体系的行为，还想直接研究单个生物大分子在原子尺度上的行为，而这是目前的实验仪器难以达到的。比如，我们希望直接研究蛋白质从伸展的肽链折叠成球形的具体过程，使用仪器手段只能收集到间接的数据，但使用软件模拟则可以形象直观的模拟出整个折叠过程，可以具体求算每个键能、键角的变化，研究某几个氨基酸残基之间的相互作用，以及对蛋白质折叠的意义。总之，目前的生物学研究需要我们利用计算机模拟生物大分子的行为，以弥补实验手段的限制，希望能自下而上地阐明生物大分子结构和功能的关系。
最后，从实际意义上讲，分子动力学模拟可以用来指导实验，提供思路和理论依据；分子动力学模拟所得结果的正确性也需要回到实验验证。这样，我们可以将分子动力学模拟和实验研究结合成一个整体，从而能够全面地，深入地研究生命现象的本质规律。 1.1 分子动力学模拟的发展 *暂缺相关文献
1.2 分子动力学模拟的基本原理 *暂缺相关文献
1.3 分子动力学模拟相关软件
随着分子动力学模拟技术的飞速发展，逐步形成了一些商品化的软件。应用于生物大分子领域的商品化分子模拟软件主要有InsightⅡ以及Sybyl，分子模拟是其中的一个重要的模块。InsightⅡ中分子动力学模块使用的是由美国哈佛大学Martin Karplus研究小组等开发的CHARMM（Chemistry at Harvard Macromolecular Mechanics），同时它本身也是一个商品化的软件。而Amber（Assisted Model Building with Energy Refinement）则是另一个非常有名分子动力学模拟软件，它是由美国UCSF的Kollman教授的课题组开发的，商业
化程度和易用性要好于CHARMM，当前版本9.0。以上两个研究小组都为其软件开发了相应的力场，并且现在已经成为分子动力学模拟的经典力场。此外免费和部分免费的软件有NAMD，Gromos，Gromacs，DL_POLY，Tinker等。
在上述软件中，我们选择NAMD作为本章的示范软件。NAMD是由美国伊利诺斯大学理论与计算生物物理研究组开发的一套分子动力学模拟软件，适用于计算生物大分子，并行计算效率非常高，可以使用Amber，CHARMM，X-PLOR，GROMACS，OPLS等多种力场，而且可以兼容Amber，CHARMM的文件格式。NAMD支持几乎所有操作系统，而且免费获取，开放源代码。如配合分子可视化、结果分析软件VMD以及格点计算软件BioCoRE则可使用更多、更强大的功能，进行更大规模的计算，可以说集众多优势于一身。
不仅如此，利用NAMD还可以进行极具特色的IMD（Interactive Molecular Dynamics，交互式分子动力学模拟）和SMD（Steered Molecular Dynamics，可控式分子动力学模拟）。在本教程中，我们将首先讲解使用NAMD进行分子动力学模拟的基本流程，然后讲解经典的结果分析方法，最后我们将简单介绍SMD的基本思想和过程。 2. NAMD分子动力学入门 2.1 软件基本设置
NAMD的最新版本是2.6版，可以从http://www.ks.uiuc.edu/Research/namd/ 免费得到（需要进行免费注册）。此外，我们还需要VMD作为分子可视化和辅助分析软件，可以从http://www.ks.uiuc.edu/Research/vmd/ 免费得到，最新版本是VMD1.85 。
NAMD安装方法：事实上NAMD是不需要安装的。请新建文件夹namd-tutorial, 在该目录中新建文件夹NAMD，下载完成NAMD2.6软件包后，将压缩文件解压到文件夹NAMD中，就可以使用。 下文中为了叙述方便，我们将默认读者的NAMD主程序位于../namd-tutorial/NAMD 目录中（安装VMD程序时可以安装到任意目录，不影响教程操作）。
此外，本教程还需要一系列教程文件。所需文件均可以从 http://www.ks.uiuc.edu/Training/Tutorials/ 下载（图）
下载完成教程所用文件后，请把所有内容解压到namd-tutorial目录下，此后的部分我们将默认教程所用文件位于../namd-tutorial 目录中。
完成上述准备之后，请打开Windows资源管理器，namd-tutorial目录的结构应该如下：（如果目录形式不一致，请务必进行调整）
该文件夹中有我们进行动力学模拟所需的所有文件。
最后，还需要交代的是，NAMD不同于我们所熟悉的大多数Windows软件：它不具有图形界面。打个比方说，我们平常使用Word，Excel，Photoshop等有图形界面的软件，好像是面对面聊天；而现在使用不具有图形界面的NAMD就像是书信往来：动力学模拟的所有参数设定都需要用户通过一个文本文件通知NAMD，NAMD进行处理计算，然后再通过许多输出文件输出结果。不借助其他软件，用户无法直接看到NAMD的工作状态。
由于进行动力学模拟的准备和结果的可视化分析，必不可少的软件是VMD，下面的讲解中也将大量用到VMD。我们假定读者已经对VMD的基本操作有一定的了解。VMD的入门教程可参见本章附录。
下面，我们将使用NAMD进行简单的分子动力学模拟，并进行初步的分析。我们将要进行动力学模拟的分子是一个76个氨基酸的小肽：泛素。
知识连接：泛素——“死亡之吻” 泛素是一个由76个氨基酸组成的高度保守的多肽链, 因其广泛分布于各类细胞而得名。泛素共价地结合于底物蛋白质的赖氨酸残基，被泛素标记的蛋白质将被特异性地识别，并在蛋白酶体中迅速降解。泛素因此得名——“死亡之吻”。因为被其标记的蛋白都摆脱不了被降解的厄运。 随着研究的进一步深入，蛋白质降解过程中泛素的枢纽作用越来越得到重视。蛋白质降解异常与许多疾病（恶性肿瘤,神经退行性疾患等）的发生密切相关。而泛素在蛋白质降解中的作用机制如能被阐明，将对解释多种疾病的发生机制和有重要意义。 Hershko、Ciechanover、Rose三名杰出科学家在泛素标记的蛋白质降解方面做出了突出贡献，他们荣获2004年度诺贝尔化学奖。 使用NAMD进行分子动力学模拟之前，我们需要为NAMD准备好各种必须的数据文件，以供NAMD使用。这些文件包括：
? 蛋白质分子的PDB文件。该文件负责储存蛋白质中所有原子的坐标。在后续课程中
我们还会了解到，PDB文件还可以储存原子运动的速度等信息。
? 蛋白质分子的PSF文件。该文件负责储存蛋白质的结构信息。注意PDB文件只记录
原子的空间位置，并不储存蛋白质中原子之间的成键情况。成键情况由PSF文件负责记录。
? 力场参数文件(force field file)。力场参数文件是分子动力学模拟的核心，文件
中的数学方程决定了原子在力场中的受力如何计算。常用的四种力场是CHAEMM, X-PLOR, AMBER 和GROMACS。NAMD可以使用以上任何一种力场进行分子动力学模拟。
? 配置文件（configuration file）配置文件的目的是告知NAMD分子动力学模拟的
各种参数，比如PDB文件和PSF文件的储存位置，结果应当储存在哪里，体系的温度等等
上述四种文件中，PDB文件通常是从蛋白质结构数据库（Protein Data Bank）中获得。力场参数文件也可以从网上下载， 而PSF文件和用户配置文件是用户根据具体要求自己生成的。下面我们将首先制作蛋白质结构文件（PSF）。 2.2 生成蛋白质结构文件（PSF）
1、单击 开始菜单→程序 → VMD，打开VMD窗口
2、在VMD主窗口中，单击 File → New Molecule 打开Molecule File Browser对话框；单击Browse…按钮，在弹出的文件浏览中找到 namd-tutorial/1-1-build 文件夹，在此文件夹中选择1UBQ.pdb，单击Load按钮载入1UBQ.pdb。 提示：关于文件后缀名 如果浏览文件时看不到“.psf”“.pdb”等后缀名，可以在“我的电脑”中选择“工具”→“文件夹选项”，在“查看”选项卡中取消“隐藏已知文件类型的扩展名”。强烈推荐读者取消这一项，因为这还涉及到下文中的许多操作。 载入之后在图形窗口（VMD 1.8.5.OpenGL Display）中应当可以看到下图（图）：
可以看到，所有的氧原子用红色表示，碳原子以天蓝色表示（碳原子所连的键也是天蓝色，所以整个蛋白骨架为天蓝色），硫原子以黄色表示。注意到没有出现氢原子，这是因为此结构是由X射线晶体衍射得来的，而X射线衍射一般得不到氢原子的精确位置。注意：蛋白周围的红点实际上是水分子，由于没有氢，所以仅显示出一个一个的氧原子。
我们只需要蛋白质分子的结构，因此下面我们将首先除去pdb文件中带有的水分子。 4、单击 Extension → TK Console 菜单项，弹出VMD Tk Console窗口。首先用cd命令改变当前目录到namd-tutorial /1-1-build 下。然后输入下列命令：
set ubq [atomselect top protein] $ubq writepdb ubqp.pdb
(每输入一行命令后按回车键，下同。另外，尤其要注意空格的有无和空格的位置，否则空格位置不对可能造成命令执行错误)
提示：VMD TK Console（VMD控制台）中改变当前目录的方法 在Windows命令行模式中和VMD TK Console 中都是用cd命令改换当前目录的。但是注意二者的使用方法不同。这里简单说明VMD TK Console 中改变当前目录的方法，Windows命令行改变目录的方法将在后面说明。 在VMD TKConsole中，改变目录的命令十分简单。无论是改变到哪一个目录，只需要输入： cd 目标目录 比如本例中，假设需要改变目录到 E:/namd-tutorial/1-1-build ，无论当前目录是什么，只需要在VMD TKConsole中输入以下命令即可： cd e:/namd-tutorial/1-1build 输入以上命令之后，VMD已经在1-1-build目录下生成了文件 ubqp.pdb。这一PDB文件仅包含泛素蛋白，不含水分子。
5、在VMD主窗口中单击1UBQ.pdb，选择Molecule→Delete Molecule菜单项删除当前分子。
6、下面我们将生成泛素蛋白的psf文件。注意：VMD组件中实际上提供了一个全自动的psf文件生成器(选择Extensions→Modeling→Automatic PSF Builder菜单项)。但我们将人工制作所需要的psf文件，以让读者明白制作的详细流程。制作时，需要使用VMD提供的psfgen软件包。
7、首先，打开写字板，输入以下内容：
package require psfgen
topology top all27_prot_lipid.inp pdbalias residue HIS HSE pdbalias atom ILE CD1 CD segment U {pdb ubqp.pdb} coordpdb ubqp.pdb U guesscoord writepdb ubq.pdb writepsf ubq.psf
8、输入完成之后，保存文件。注意文件保存在1-1-build目录中，文件名为ubq.pgn，文件类型选择文本文档。然后退出写字板。这样我们便制作了pgn文件，这一文件可以被psfgen软件包所识别，并处理成我们想要的psf文件。我们需要在VMD中使用该文件调用psfgen数据包
下面我们详细介绍一下刚刚输入的每一行命令的意义：
package require psfgen：通知VMD我们将要调用psfgen数据包
topology top all27_prot_lipid.inp：载入拓扑文件 top_all27_prot_lipid.inp pdbalias residue HIS HSE：改变组氨酸残基名，使得残基名称能够和拓扑文件中的一致。在pdb文件中组氨酸残基名是HIS,而在拓扑文件中组氨酸残基名为 HSE, HSD, HSP 三种。分别对应组氨酸的三个不同的带电荷形式。
pdbalias atom ILE CD1 CD：改变异亮氨酸中的原子名。pdb文件中异亮氨酸δ碳的
名称为CD1，而拓扑文件中原子名应该为CD。
segment U {pdb ubqp.pdb}：生成一个集合（segment）U，包含ubqp.pdb中的所有原
子。
coordpdb ubqp.pdb U：从ubqp.pdb中读取坐标，比较各个原子的名称是否对应，然
后旧的集合名被改换成新的名称“U”。
guesscoord：根据拓扑文件推测缺少的原子（氢原子）的空间位置。
writepdb ubq.pdb：生成新的pdb文件，包含所有原子的坐标，包括刚刚推测出的氢
原子。
writepsf ubq.psf：生成psf文件，该文件包含蛋白结构的全部信息。
知识链接：组氨酸的三种离子模式
知识链接：PDB文件中原子的命名方式 9、如果刚刚关闭了VMD，则重新打开，改变目录至1-1-build。然后输入以下命令： source ubq.pgn
这样我们就成功得到了含有氢原子的psf文件。同时，可以看到VMD TKConsole中显示出系统返回的信息。信息显示我们的系统中有1231个原子，631个原子的坐标是推测的（图）。
现在在你的1-1-build文件夹下应当有ubq.pdb和ubq.psf两个文件。到此为止，我们已经成功制作了下一步分子动力学模拟所需的psf文件。 2.2 蛋白质的溶质化
显然在真实情况下，蛋白质不是在真空中存在下面。所以我们需要把蛋白质放入一个水环境中，以更真实的模拟生物体内的环境。我们可以使用两种水体环境进行动力学模拟：
? 球状水体（water sphere）。水体包围蛋白质，四周则是真空，动力学模拟时没
有周期性边界条件（periodic boundary condition）
? 立方水体(water box)。立方水体是正六面体形状的水体（不一定是正立方体）。
使用立方水体需要我们设定周期性边界条件。
2.2.1 生成球状水体（water sphere）
我们将使用一个脚本文件建立球状水体。脚本文件在1-1-build目录下，文件名是wat_sphere.tcl。
1、如果刚刚关闭了VMD，则重新打开，改变目录至1-1-build。然后输入以下命令： source wat_sphere.tcl
输入之后VMD将会调用脚本文件，之后VMD会反馈一系列信息（图），
2、由所给的信息我们可以看出，VMD生成了两个文件：ubq_ws.pdb文件和ubq_ws.psf。 在最后两行还给出了生成的球状水体的质心坐标（center of mass of sphere）和球状水体的半径(radius of sphere)，精确到小数点后第十位。记下这些数字，以后我们还会用到。
这时在图形窗口中却会出现一个立方水体包围的蛋白质分子（图）。不过没有关系，在VMD主窗口中可见分子名为del_water，并不是我们所要的结果。我们的最终结果已经储存在1-1-build中，文件名分别为ubq_ws.pdb和ubq_ws.psf。
3、下面我们将看一下生成的球形水体究竟是什么样子的。在主窗口中单击del_water 分子，选择Molecule → Delete Molecule 菜单项删除该分子；然后选择 File → New Molecule，单击Browse 按钮，在1-1-build目录下找到ubq_ws.pdb文件，单击Load载入该蛋白，可以看到球状水体包围的蛋白（如图）。说明我们已经成功地生成了球状水体包围的泛素分子。
2.2.2 生成立方水体（water box）
下面我们将把泛素放入一个立方体状的水环境中。我们使用的是VMD提供的solvate软件包。该软件包位于VMD的/plugins/noarch/tcl目录下。不过我们不需要自己找到它。只要通知VMD我们将使用该软件包，VMD就会载入它。
1、打开VMD，选择 Extensions→TK Console菜单项，在VMD TKConsole 窗口中输入：
package require solvate
这时VMD就会载入solvate软件包。窗口返回数字：1.2 说明我们所使用的软件包是solvate 1.2。确保当前目录是1-1-build，否则用cd命令改变当前目录至1-1-build，然后输入:
solvate ubq.psf ubq.pdb –t 5 –o ubq_wb
等待运行结束，VMD就调用solvate将ubq.pdb和ubq.psf所储存的蛋白放入一个立方水体中。在图形窗口可以见到一个立方形的水体包围蛋白（如图）。
参数 –t 5 通知程序如何确定立方体的各边长。方法是在每个坐标方向上选择坐标最大的那个原子，然后再延伸5A，即为该方向立方体面的边界。注意：生成的立方水体并不一定是正立方体。各边长取决于坐标最大（距离原点最远）的原子的位置。还有一个参数 –o ubq_wb是为了通知程序生成的文件名。运行结束后我们得到的两个文件就是ubq_wb.psf 和 ubq_wb.pdb。
2、在VMD TkConsole中输入： set everyone [atomselect top all] measure minmax $everyone
这时返回的数值是整个体系中离原点最近的点和最远的点的坐标。我们需要的是整个立方体的中心，可以自己计算也可以用下面的命令：
measure center $everyone
这时返回的三个数值就是体系的中心。记下这三个数值，我们以后还会用到。
返回值如图：
在开始下一节之前，我们要将生成的文件拷贝到common公用目录下以方便访问。在Windows资源管理器中找到1-1-build目录，按Ctrl选择以下六个文件：ubq.pdb, ubq.psf, ubq_ws.pdb, ubq_ws.psf, ubq_wb.pdb, ubq_ws.psf，然后把它们拷贝到namd-tutorial/common目录下。 提示： 此处生成的立方水体事实上过小了。在实际应用时，应当保证水体足够大，以防止蛋白在拉伸运动时超出水环境。也要避免在使用周期边界条件时蛋白和四周各个单元的蛋白镜像相碰撞。周期边界条件的详细知识见下文。 此外，还应当注意的是在实际应用时应当在水环境中放入离子。特别是当蛋白质的净电荷不为0时，更应当设定离子数目以使得整个体系是中性的。在放入离子时应当将它们放在体系中静电势能的最低点，以节省计算时间。因为离子总会向势能最低点自发运动。 2.3 球状水体中泛素（Ubiquitin）的分子动力学模拟
在这一节中，我们将要对球状水体的泛素分子进行最简单的动力学模拟。 首先，我们要进行的分子动力学模拟的目的是什么？我们将泛素放入球状水体中，水体周围是真空，然后NAMD会根据我们设定好的温度值按照Boltzmann-Maxwell分子速率分布给各原子赋予一定的初始速度，接下来就是要根据牛顿力学方程，求解个水分子以及蛋白质中各原子的运动轨迹。我们得到的结果，就是模拟泛素这一小肽在溶液状态下的运动状态。
知识链接：能量最小化和能量平衡（Minimization and Equilibration） 本次动力学模拟实际包括两个过程：能量最小化和能量平衡（Minimization and Equilibration）。能量最小化时，NAMD设定各原子的速度为0，然后不断改变各个原子的相对位置并计算体系总能量，搜索最低势能点，作为分子动力学模拟的初始状态。这一过程是不记录原子运动轨迹的。因为原子的位置改变只是因为NAMD需要搜索最低能量状态，而不是真实的相互作用引起的运动。 能量平衡是让蛋白质和水分子在设定好的环境温度（即原子的速度）下相互作用，达到能量平衡分配，整个体系达到稳定状态（熵达到极大值）。 为什么需要首先进行能量最小化？这是因为我们提供的体系有可能包含极度扭曲，拉伸或压缩变形的键和键角。它们是解析结构或同源模建时引入的错误结构，含有很高的能量。如果不首先进行最小化，直接进行能量平衡，蛋白质会和水分子相互作用，恢复伸展状态，释放掉这些错误结构中的高能量。这一过程是没有意义的——因为它是错误结构引起的反应，并不是蛋白质在溶液中的真实状态。从而就浪费了计算时间。不仅如此，能量释放引起的剧烈运动和相互作用最终可能使得蛋白质的行为不符合溶液中的真实行为。所以有必要在能量平衡之前，首先人为搜索能量最低点，作为分子动力学模拟的初始状态。 一般地，分子动力学模拟包括多个能量最小化和平衡过程。通常我们会首先将蛋白质固定而仅允许水分子运动，进行能量最小化和能量平衡；然后允许蛋白质和水分子同时运动，再次经历能量最小化和能量平衡这一循环。第一步的目的是使水分子达到能量最小，这通常是一个很快的过程。然后再放开蛋白质，使整个系统达到能量最小。这样可以减小计算量，并防止由于一开始蛋白结构很不稳定而结果产生假象。 在上一节中我们已经获得了所需要的ubq_ws.pdb 和ubq_ws.psf 两个文件。对照本章开始提到NAMD所需的四个文件知，还需要有配置文件就可以提交NAMD进行动力学模拟了。（力场参数文件在common文件夹中）。下面我们将首先得到配置文件。 2.3.1 配置文件
前面我们把使用NAMD比作写信，这里的“信”就是指的配置文件。配置文件记录了进行动力学模拟所需的全部参数和设置，NAMD只要得到这一文件就可以按照相应指令进行操作。对于本次动力学模拟，在namd-tutorial/1-2-sphere 目录下可以得到已经预先制作完毕的配置文件。下面我们将要仔细讲述文件的内容。
打开写字板，在菜单中选择 文件→打开，找到1-2-sphere 目录，文件类型选择全部文件（图），然后打开文件ubq_ws_eq.conf。这个文件看起来好像很复杂，但是我们会仔细分析讲解每一部分的含义。
注意：在配置文件中，每一行开头如果是“#”，则本行内容会被当作注释对待，NAMD会忽略其中的内容。因此为了便于区分，我们用#把文件分割成几大部分。如第一部分是：
############################################ ## JOB DESCRIPTION ## ############################################ 意思是这一部分是对所提交工作的描述。
1、大体浏览一下，可以发现整个文件被分成了以下几部分： ? 工作描述（job description） ? 可调参数 （adjustable parameters） ? 动力学模拟参数（simulation parameters） ? 附加参数（extra parameters） ? 执行脚本 （execution script）
它完整的记录了输入的蛋白质结构文件（pdb和psf文件）的位置，输出结果文件的文件名，以及动力学模拟时的环境温度，截止点，步长等各种参数。在进行动力学模拟时只需要提供给NAMD一个配置文件，NAMD就可以找到输入文件，调整好各种参数，按照要求进行动力学模拟之后输出结果。
2、然后，首先我们来看第一部分：Job Description。这一部分每一行开头都有#，因此只包括注释。它描述的是这一配置文件的目的：Minimization and Equilibration of Ubiquitin in a Water Sphere。就相当于一片文章的题目。
3、Adjustable Parameters 这一部分包括5项参数： ? structure: 给出调用的psf文件的位置
? coordinates：给出调用的坐标文件（即pdb文件）的位置
? set temperature 310: 定义一个变量temperature，并赋值310。以后如果要使
用环境温度值310，只需要用$temperature 代替。这和c语言中的预处理命令#define有些类似。
? set outputname ubq_ws_eq：新建一个变量outputname，并赋值ubq_ws_eq。作
用同上。
? firsttimestep:设定动力学模拟时起始timestep的数值。在重新开始进行一个
被中断的动力学模拟时，这一设定是非常有用的。如果前一次动力学模拟结束时的timestep是533，那么这一次的起始值显然应该是534。
4、Simulation Parameters 这一部分包括许多参数，可以分成以下几部分： ? Input
- paraTypeCharmm: 说明参数文件是否是CHARMM力场格式的。设置为on。 - parameters: 从给出的力场参数文件中调用参数（此例中，力场参数文件为../common/par_all27_prot_lipid.inp）
- temperature: 设定环境的起始温度（K）。如上所述，在这里$temperature 相当于310。设定这一数值后，NAMD会根据Maxwell分子速率分布给体系中的分子分配运动速率。
? Force-Field Parameters
- exclude: 说明哪一种原子-原子相互作用可以忽略。这里的设定值是 scaled1-4。成键相邻原子的编号方式见图。scaled1-4就是说如图中的原子1-2，1-3，2-3之间的相互作用被完全忽略，而原子1-4的相互作用被弱化。
- 1-4scaling: 刚刚提到了原子1-4之间的相互作用会被弱化。这个参数就是为了说明弱化的程度。取值在0～1之间，0表示完全忽略，1表示不进行弱化。 - cutoff：设定范德华力和静电力的截止点。如果不设定此值，NAMD会计算整个体系中任意两个原子的范德华力和静电力相互作用，这显然是没有必要的。注意：如果Particle Mesh Ewald Sum 设定为on，cutoff的定义就会改变，在此不详细叙述。
- switching：设定是否使用过渡函数（switching function），使得在截止点处
范德华力和静电力不会突然降低至0，而是平滑的过渡至0。
- switchdist：设定在哪一点静电力和范德华力函数开始使用过渡函数修正（switch function）以使这两个函数可以平滑过渡，在cutoff处降低为0。 - pairlistdist: 这一设定是为了使得计算更快进行。如果不设定这个值，对于体系中的某个原子，NAMD需要遍历搜索整个体系以找出和该原子有相互作用的所有其他原子。设定之后，在计算某个原子的受力时，NAMD将只搜索设定范围之内的原子。设定值的单位是A。注意这个值必须要大于cutoff值。
图是以上概念的图示说明。 ? Integrator Parameters
- timestep：说明所使用的步长数值。分子动力学模拟的基本原理还是求解牛顿力学方程，但是并不能做到连续求解，而只能每隔一段间隔求解一次，最后生成原子的运动轨迹。步长值就是求解的时间间隔。以一个飞秒（fs,femtoseconds）为单位。2.0即为2fs。
- rigidBonds：设定与氢原子相连的哪一种键是刚性的（不会来回振动）。这里设定值是all，说明所有和氢原子相连的键都被认为是不振动的。
知识链接：Rigid Bonds 为什么要设定RigidBonds？这是因为我们设定的步长是2飞秒。在分子动力学模拟时，键的转动，振动，原子的位移等等速度并不相同。而步长数值显然应该由最快的那一种运动的时间尺度决定。在各种运动形式中，键长的伸缩和键角的扭曲是最快的。键长振动一般是每10-100飞秒一次。其中，最快的当然是与氢原子相连的键长的振动，一般是10飞秒一次，而我们的步长是2飞秒，几乎在一个数量级上。因此无法精确描述这样的键长振动。所以需要先设定认为这些键不振动。其实也相当于取键长伸缩振动的平均值作为固定键长。 大分子的功能和行为一般与较慢的分子构象变化和分子运动关系密切，但和快速的原子振动关系较小。所以认定键长振动不存在也是可以接受的，只是对于精确的分子动力学模拟而言应当尽量避免。对于任何的分子动力学模拟，步长应该是体系中最快运动周期的1/10以下。
- nonbondedFreq: 设定每隔多少步长计算一次非成键相互作用（nonbonded interactions）。适当调整这个值可以节约计算时间。
- fullElectFrequency: 设定每个多少步长计算一次总体静电相互作用（full electrostatic interactions）。
- stepspercycle: 前面提到过，每个原子都有一个pair list，即和它有相互作用的所有原子的列表。这个列表显然是动态变化的。列表更新的周期叫做一个循环（cycle）。这个值设定的是每多少步长更新一次列表，完成一次循环。 ? Constant Temperature Control
- langevin: 设定动力学模拟时是否使用Langevin 动力学。这里设定为on。 - langevinTemp: 设定一个温度值，使用Langevin 动力学将原子保持在恒定的该
温度。
- langvinHydrogen: 设定是否对于氢原子也应用langevin动力学。 ? Output
- outputName：每进行一次动力学模拟，NAMD会输出多个文件。这个参数设定这些文件的前缀名（如ubq.pdb，ubq就是前缀名）都为ubq_ws_eq。NAMD输出的文件包括：一个后缀名为“.coor”的文件，储存经过动力学模拟后的所有原子的坐标；一个后缀名为“.vel”的文件，储存系统动力学模拟结束时所有原子的瞬时速度。所以运行结束后我们可以得到两个文件：ubq_ws_eq.coor 和 ubq_ws_eq.vel。 - restartfreq: 在进行分子动力学模拟时，NAMD还会创建恢复文件（restart file），类似于Word的自动保存，使得用户在动力学模拟意外停止的时候可以用恢复文件继续进行模拟。这个参数就是设定每过多长个步长自动保存一次，生成一个恢复文件。恢复文件的后缀名是 “.restart”，表示刚刚生成的恢复文件；以及“.restart.old”，是前一次保存的恢复文件。
- dcdfreq：dcd文件记录的就是每一个原子的运动轨迹。记录方法是，NAMD每隔一
定时间间隔就将所有原子的坐标写入一次dcd文件。而这个参数就是设定写入的时间间隔。当然，dcd文件会随着模拟的进行而越来越大，如果写入很频繁或者模拟进行的时间很长，就会得到一个很大的dcd文件。另外，如果不需要得到模拟后的轨迹也可以不设定这一参数，这样NAMD将不会生成dcd文件。
除了以上叙述的这些输出文件，namd还会产生一个日志文件，后缀名是“.log”。这一文件的内容将在以后的内容讲到。
- outputEnergies: 设定每隔多少步在日志文件中输出系统的各种能量（每种立场如范德华力，静电力分别对应一种能量）。这里我们的设定是每隔100步输出一次。 - outputPressure: 同样地，这个值是为了设定每多少步在日志文件中输出一次系统压力。
5、附加参数（Extra Parameters） ? Spherical Boundary Conditions
- sphericalBC: 设定是否要设置球形边界条件。
- sphericalBCcenter: 设定球形体系的中心。输入你记下的球状水体中心的坐标。在这里我们已经给出了所需要的坐标值。
为了使球形边界条件可以维持，需要设定一个边界势能，使得球状水体得以保持形状而不会扩散到真空中去。以下三行参数就是设定了边界势能。
- sphericalBCr1: 设定第一个边界势能起作用的起始半径。以A为单位。 - sphericalBCk1: 设定边界势能的force constant。单位是 kcal/mol·A。 - sphericalBCexp1: 设定边界势能函数方程的指数值。必须是正偶数。 6、执行脚本（Execution Script）:最后一个部分，包含三个参数设定：
? Minimization: 在本次模拟时，一开始NAMD将不断改变各个原子的位置，搜索
整个体系势能的最低点（此时个原子的动能均为0），以作为动力学开始的初态。这就是能量最小化（minimization）
- minimize:。这一参数设定的是能量最小化时反复改变原子位置的次数。 - reinitvels: 能量最小化时，各个原子的速度还是0。这个参数设定的是能量最小化完成之后体系升温所至的温度。在这个例子中是$temperature,即为310K。 ? run：设置分子动力学模拟进行的时间。以步长为单位，这里设定为2500步，即
2500 fs x 2 = 5000fs，或5 ns（nanoseconds，纳秒）。
7、现在关闭写字板。 2.3.2 进行动力学模拟
前面说过，NAMD没有图形界面，需要打开命令行之后找到namd所在目录，然后输入namd2运行程序。运行时输入命令的格式是：
namd2 配置文件位置 > 输出日志文件位置
单击 开始→运行，输入cmd后回车，打开命令行窗口。使用cd命令改换目录到namd所在的目录，（本教程中默认目录是namd-tutorial/namd），然后输入： namd2 ../1-2-sphere/ubq_ws_eq.conf > ../1-2-sphere/ubq_ws_eq.log （注意斜线“/”的方向，不可以输入“\\”）
回车之后，系统不返回任何信息，光标位于下一行行首并不断闪动，说明NAMD已经开始运行。按Ctrl+Alt+Del 打开Windows 任务管理器，在“进程”选项卡中可以看到有namd2这一进程运行，并且CPU占用率一直是100% 。为了使动力学模拟能够顺利完成，在此期间不要进行任何操作。
提示：Windows命令行中改变当前目录的方法 在Windows命令行模式中和VMD TK Console 中都是用cd命令改换当前目录的。但是注意二者的使用方法不同。尤其注意的是：在Windows命令行中以“\\”而不是 “/”作为目录分隔符。 在Windows命令行模式中，假设当前目录是C:\\Windows： 1、改换到当前盘C盘的任何目录(以xxx表示)只需要输入： cd c:\\xxx\\xxx\\ 2、改换到其它盘中的目录（以D:\\xxx\\xxx表示），可以先在当前目录输入： cd d:\\xxx\\xxx\\ 注意此时目录不会改变。只需要再输入： d: 就可以发现目录改变到了目的目录 另外，“..”表示的是“当前目录的上一级目录”。
提示：输入和输出文件 这里, ../1-2-sphere/ubq_ws_eq.conf是我们的配置文件的所在位置，../1-2-sphere/ubq_ws_eq.log是我们想要得到的日志文件的所在位置。注意：除了日志文件“.log”，NAMD还会输出其他文件，但不需要指定NAMD输出的其他文件(如ubq_ws_eq.coor,ubq_ws_eq.vex等)的位置，因为它们默认输出到配置文件ubq_ws_eq.conf所在文件夹中(namd-tutorial/1-2-sphere)，而日志文件ubq_ws_eq.log必须指定位置，否则将会默认输出到NAMD主程序所在的文件夹中（namd-tutorial/namd）。 大约10至20min后程序运行完毕，命令行窗口中，光标跳到下一行，并显示当前目录，但是不会在窗口中返回任何信息。
2.4 立方水体中泛素的分子动力学模拟
在这一节中，我们将对立方水体中的泛素分子进行分子动力学模拟。
对立方水体进行动力学模拟的时候，很大的一点不同之处就是需要设定周期性边界条件。如果不设定，那么NAMD将认为立方水体四周是真空。进行模拟的时候由于水分子的表面作用，真空中的立方体会逐渐变成球形（因为球体表面积最小）。这样我们设定立方水体也就没有什么实际意义了。
之所以要把水体设定成立方体，是因为我们要将该水体视作结构单元（cell），然后在模拟的时候假定立方体四周像搭积木一样包围着和该水体同样的一个个立方体。换句话说，我们所研究的立方水体只是一个巨大、连续水体的一个小单元。很显然，球体不可能作为结构单元。只要想像一下为什么砖块都是立方体而不是球体就明白了。 2.4.1 得到配置文件
和上一节球形水体的分子动力学模拟类似，我们提供了动力学模拟所需的配置文件。 由于整个体系的边界（立方体单元的边界）是周期性重复的，所以我们在配置文件中要设定相应的周期性边界条件（Periodic Boundary Conditions）。使用写字板打开namd-tutorial/1-3-box 目录下的配置文件ubq_wb_eq.conf，我们将要讲解这一配置文件和上一节进行球状水体动力学模拟时的配置文件的不同之处。
1、比较一下可以知道，唯一的不同之处在Simulation Parameters一部分。新增加了下面几项：
? Periodic Boundary Conditions
- 首先设定了三个单元基向量（cell basis vectors）。这三个向量是：
cellBasisVector1, cellBasisVector2, cellBasisVector3。后面的数值是它们的坐标。从坐标中可以知道，这三个向量是两两垂直的。比如cellBasisVector1的坐标是{42，0，0}。平行于x轴，而cellBasisVector2则为{0，44，0}，平行于y轴。三个向量设定之后，就可以建立起一个三维立方体空间。
- cellOrigin: 设定初始立方单元的中心。其他的立方单元就以此为中心，开始周期性的重复。因此在2.2一节进行蛋白质的溶质化的时候，我们需要记录下生成的立方水体的中心（这里在配置文件中我们已经给出了中心坐标）
- wrapWater: 设定开启周期性边界时需要设定这一参数。这个参数的意思是：如果出于边界区域的水分子运动超出了立方体边界，则对它的坐标做关于边界的镜像变换，从而使水分子回到了立方体单元之中。这样可以防止立方体边界区域的水分子不断向外扩散。
- wrapAll:顾名思义，就是对所有超出边界的原子进行镜像。当这一参数设定为on时，wrapWater可以不用设定。 ? PME
PME的全称是Particle Mesh Ewald。当设定了周期性边界的时候，使用PME可以方便地处理静电力相互作用。particle mesh是一个假想的三维网格，系统中各个原子所带的电荷被分布到这一网格上，然后可以计算出网格上面的各点的静电势，从而计算出各个原子的受力情况。因此，网格的尺寸应当足够精细，可以精确反映电荷在系统中的分布。
- PME：设定PME为on“开”或off“关”。
- PMEGridSizeX: 设定网格尺寸。注意：虽然参数名中有一个“X”，但并不是指X轴方向的尺寸，而是沿单元基向量1（cellBasisVector1）方向的尺寸。尺寸的定义方法是：沿cellBasisVector1所定义的方向作直线，直线跨越的网格的节点的数目就是网格的尺寸。所以，这一个参数设定越大，网格就越密，越能够精确地反映电荷分布的真实情况，但计算耗时也越多。一般来说，尺寸略小于10nm的网格能够比较好的反映生物系统中的电荷分布，因为生物系统中原子的最近距离大约是10nm左右。 - PMEGridSizeY: 设定网格沿单元基向量2（cellBasisVector2）方向的尺寸。定义方式如上。
- PMEGridSizeZ: 设定网格沿单元基向量3（cellBasisVector3）方向的尺寸。定
义方式如上。
需要注意的是，单元基向量1，2，3的长度略有不同，分别为42，44，47。因此虽然
PMEGridSizeX,Y,Z设定都为32，实际上每个方向上的网格尺寸也是不同的。
? Constant Pressure Control (variable volume)
- useGroupPressure：设定是否在计算系统压力时考虑每个氢原子之间的相互作用。NAMD根据原子之间的作用力和各个原子的动能计算系统的压力。如果rigidBonds 设定为on，这一参数必须设定为yes。
- useFlexibleCell：设定是否允许周期性立方体单元的边长独立变化。 - useConstantArea：设定是否保持立方单元的x-y横截面积不变，而沿z轴的边长发生变化。
- langevinPiston：设定是否使用Langevin piston 控制系统的压力。 如果使用Langevin piston, 那么需要设定以下选项：
- langevinPistonTarget： 设定使用Langevin piston维持系统压力为多少巴（1大气压 = 1.013巴）
- langevinPistonPeriod： 设定使用Langevin piston时的振荡时间常数（oscillation time constant）。单位是飞秒。
- langevinPistonDecay： 设定使用Langevin piston时的衰减时间常数（damping time constant）。单位是飞秒。
- langevinPistonTemp： 设定使用Langevin piston时的噪音温度（noise temperature）。单位是K。这个温度设定应和环境温度相等。 ? Output
- xstFreq：xst代表“extended system trajectory”文件。这一文件储存了周期性单元参数，特别是记录了周期性单元边界在动力学模拟时的变化。这一参数设定每隔多少步记录一次周期性单元参数。设定这一参数后，结果将会输出3个文件:2个恢复文件和1个结果文件。 2.4.2 进行动力学模拟
打开命令行，用cd命令改换目录到NAMD所在的目录下(namd-tutorial/namd)，输入： namd2 ../1-3-box/ubq_wb_eq.conf > ../1-3-box/ubq_wb_eq.log 动力学模拟应当在20min之内结束。
至此，我们已经完成了两次动力学模拟。
2.5 简单的结果分析
在这一节中，我们将对球状水体中泛素的分子动力学模拟进行简单的结果分析。立方水体中的结果分析留待下一章完成。 2.5.1 获得结果
球状水体的动力学模拟的结果在1-2-sphere目录中。在windows资源管理器中找到该目录，该目录中应当包括以下输出文件，共11个（注意：如果选择隐藏已知文件的扩展名，那么ubq_ws_eq.log文件将被显示成“文本文件”，没有扩展名“.log”）。
? ubq_ws_eq.log ? ubq_ws_eq.coor ? ubq_ws_eq.vel ? ubq_ws_eq.xsc ? ubq_ws_eq.dcd
? ubq_ws_eq.restart.coor ? ubq_ws_eq.restart.vel ? ubq_ws_eq.restart.xsc ? ubq_ws_eq.restart.coor.old ? ubq_ws_eq.restart.vel.old ? ubq_ws_eq.restart.xsc.old
此外该目录下还包括配置文件ubq_ws_eq.conf
在这11个输出文件中，有7个文件是以二进制编码的形式存储的，以便精确记录各个原子的坐标和速度。非二进制编码的文件包括：
? ubq_ws_eq.xsc, ubq_ws_eq.restart.xsc, ubq_ws_eq.restart.xsc.old。它们是扩展
系统配置文件。这些文件存储的是整个系统的周期性单元维度（periodic cell dimension）以及存储时间。这里我们使用球形水体，没有设定周期性边界条件，因此这些文件事实上没有用。
? ubq_ws_eq.log：日志文件。从不同的立场计算出的能量就储存在这个文件中。由于日
志文件非常重要，我们将仔细讲述其中的内容。下面使用写字板打开日志文件。方法是右键单击日志文件，在弹出的菜单中选择 打开方式→写字板 。
（如果直接双击打开，默认是使用记事本打开的。）
1、日志文件一共分成三部分：系统信息，能量最小化，能量平衡。下面我们将分别讲述： 2、第一部分系统信息以“Info：”为每行的开头（图）。它包括了运行动力学模拟所需的各种参数，以及其他关于整个系统的信息，比如原子数，各种键的类型和数目，总质量和总电荷量。注意：一直到最后一个“Info：”都是第一部分的内容，记录了分子动力学模拟过程的基本信息。
2、找到最后一个“Info：”，之后就是第二部分，给出的是关于系统能量最小化的信息（图）。
? 第二部分开头是“TCL: Minimizing for 100 steps”，说明能量最小化每循环进行100
步。每次循环输出一次图中矩形标注出的内容，包括：
PRESSURE
GPRESSURE ETITLE ENERGY INITIAL STEP GRADIENT TOLERANCE
? PRESSURE和GPRESSURE：给出的是系统的压力。本次动力学模拟的日志文件中，压力值
除第一个为步数（即系统所走过的的timestep）外，都显示为一排“0”，这是因为只有在定义周期边界条件之后，NAMD才会计算压力。在立方水体动力学模拟的日志文件中，我们将能看到计算得的压力值。
? ENTITLE和ENERGY: 给出的是力场相互作用能量的列表。ENTITLE后面列出的是相互作用
能的名称。TS是TimeStep，记录了动力学模拟进行了多少步；后面的能量包括BOND（键能），ANGLE（角度扭曲能）等等。下面一行以ENERGY开头，列出的是对应的各项能量的精确数值。
? INITIAL STEP 和GRADIENT TOLERANCE：这两项是记录能量最小化过程的两个参数。
GRADIENT TOLERANCE在能量最小化过程中会一直下降，可以用于粗略地判断是否能量已经达到最小。
3、向下翻页，一直找到“TCL：Running for 2500 steps”，从这一行向后就是第三部分：能量平衡（图）。
图中a中的数值为能量最小化终态时的数值，b中的数值为能量平衡开始100步时的输出数据，注意到二者在各项能量上已经有较大差异了。能量平衡时数据每100步输出一次，这是我们在配置文件中设定的。参见2.3.1配置文件中 –outputEnergies 和 –outputPressure两项，我们设定的数值是100。
b中有三行数据，分别记录了压力（PRESSURE, GPRESSURE） 和能量（ENERGY）。压力仍然为0。拖动左下方滚动条向右，可以看到动能（KINETIC）和温度（TEMP）值（图）。可以看到两项值都已经不为0了。而之前的所有输出的动能和温度都为0。
日志文件的结尾显示：
WRITING EXTENDED SYSTEM TO OUTPUT FILE AT STEP 2600 WRITING COORDINATES TO OUTPUT FILE AT STEP 2600 CLOSING COORDINATE DCD FILE
WRITING VELOCITIES TO OUTPUT FILE AT STEP 2600 ==========================================
WallClock: 569.015991 CPUTime: 569.015991 Memory: 0 kB Program finished.
说明程序正常结束了，最终进行了2600步运算（不是我们在配置文件中设定的2500是因为能量最小化100步也包括在内）。并给出了CPU时间以供参考。
2.5.2 计算整个蛋白的RMSD值
在动力学模拟分析中，很基本的一项分析就是RMSD值的求算。RMSD值即均方根偏差（Root Mean Square Deviation）。在统计学上，这个量就相当于标准差，反映的是数据偏离平均值的程度。在蛋白质结构解析，模建，结构联配（structure alignment）以及分子动力学模拟中，RMSD值是非常常用的一项参数，用于衡量原子偏离比对位置的程度。
在本例中，RMSD值反映的是泛素各个部分原子偏离平均位置的程度，也就是各原子运动
幅度的大小。RMSD值越大，说明该原子的运动的空间范围越大，原子的空间位阻也就越小。原子α在某一段时间Nt内的RMSD值求算公式为：
上式中，Nα 是参与比较的原子数目，r?(tj)是原子α在tj时刻的位置，
Nt是时间段，以步长为单位。
VMD提供了很多分析工具。其中一个工具NAMD Plot 插件可以用于处理输出文件，然后由进行作图。我们下面就用这一工具分析动力学模拟过程中整个蛋白的RMSD值的变化。这个值是各原子RMSD值的总和。计算时，我们需要使用动力学模拟过程输出的轨迹文件：ubq_ws_eq.dcd
首先，让我们简单看一看NAMD plot是怎样作图的。 1、运行VMD主程序。
2、在VMD主窗口选择菜单项Extensions→Analysis→NAMD Plot。
3、在 NAMD Plot窗口选择File→Select NAMD Log File，在弹出的文件浏览窗口中找到namd-tutorial/1-2-sphere/bq_ws_eq.log，打开这个文件。
4、在NAMD plot 窗口中显示出很多热力学参量（图）。其中TS代表Time Step，BOND代表键能等等。这些变量都可以用于作图。例如：点击选择参数TEMP，然后选择 File → Plot Selected Data。这样我们就得到了随着时间变化温度的变化关系。
5、你可以试着用其他热力学数值对Time Step作图，熟悉使用NAMD plot。然后关闭NAMD plot。
6、首先，打开VMD，选择File→New Molecule菜单项，在弹出的文件浏览中找到 common目录下的文件usq_ws.psf，点击按钮Load。这样我们首先载入了泛素蛋白的结构信息。（注意图形窗口中是空白的，并没有显示出蛋白结构，因为psf文件不存储各个原子的坐标）
7、此时Molecule File Browser窗口上端一栏应该显示Load files for: 0: ubq_ws.psf （图）说明如果再次载入新文件，所载入的文件将与usq_wb.psf结合起来一同显示。我们将载入原子运动轨迹文件，和psf结构文件配合即可显示出蛋白质原子的运动轨迹。因此点击Browse，找到1-2-sphere文件夹下的ubq_wb_eq.dcd，点击Load，这时关闭Molecule File Browser窗口。可以在图形窗口中看到已经载入的泛素和球状水体。
8、我们将要使用一个脚本（script）计算RMSD值。用写字板打开1-2-sphere目录下的脚本rmsd.tcl，可以看到脚本的内容为：
set outfile [open rmsd.dat w]; set nf [molinfo top get numframes]
set frame0 [atomselect top \set sel [atomselect top \# rmsd calculation loop
for {set i 1 } {$i < $nf } { incr i } { $sel frame $i
$sel move [measure fit $sel $frame0] puts $outfile \}
close $outfile 脚本的每一句含义解释如下：
? ?
set outfile [open rmsd.dat w]： 新建文件rmsd.dat作为输出文件
set nf [molinfo top get numframes]：声明变量nf，并把轨迹文件中的帧数(frame)
赋值给变量nf ?
set frame0 [atomselect top \声明变量frame0，frame0储存第一帧中蛋白质骨架非氢原子的位置。 \and backbone and noh\表示我们选择的是蛋白质骨架上的所有原子，并且不包含氢。定义这个变量是为了以第一帧作为基准,其他各帧都和第一帧比较计算RMSD值(注意我们在定义的时候,
RMSD值是以原子运动的平均位置为基准的。这里为了计算方便我们以第一帧为基准) ?
set sel [atomselect top \：声明变量sel，变量暂时没有确定的值。以后每个循环依次将各帧的蛋白质骨架中的非氢原子放入这一变量中。
? ?
# rmsd calculation loop：注释：RMSD值循环计算开始
for {set i 1 } {$i < $nf } { incr i } {：和C语言中的for语句用法一样，这一语句表示：定义变量i=1，每次循环i的增量为1，i > nf 时跳出循环。相当于C语言中的语句：for（int i=1;i < nf;i++） 此后由{ }括起来的内容就是循环体。每循环执行一次
? ? ?
$sel frame $i : 将变量sel设定为所要比较的那一帧（第i帧）。
$sel move [measure fit $sel $frame0] ：将第i帧和第1帧的蛋白骨架联配（align）。 puts $outfile \$sel $frame0] :计算第i帧和第1帧之间比较的RMSD值，并输出到我们在第一条语句中新建的文件：rmsd.dat。
9、如果刚才关闭了VMD，请再次打开，重新载入common/ubq_ws.psf 和 1-2-sphere/ubq_ws_eq.dcd。下面，在VMD主窗口中选择Extensions →TK Console打开控制台，在控制台中以cd命令转换目录到namd-tutorial/1-2-sphere， 然后输入source rmsd.tcl 运行脚本。脚本会在1-2-sphere目录下生成文件rmsd.dat。这一文件记录了蛋白骨架的RMSD值随时间的变化。
得到rmsd.dat数据文件后，我们可以用多种绘图软件进一步作图。这里以最常见的Excel为例，做出蛋白质骨架RMSD随时间的变化。
10、关闭VMD，打开Excel，选择 “文件”→“打开” ，在弹出的文件浏览中找到namd-tutorial/1-2-sphere目录，文件类型选择所有文件，找到rmsd.dat打开它（如果选择不显示已知文件类型的扩展名，可能看不到后缀“.dat”）。
11、这时会弹出“文本导入向导”对话框（图）。直接单击“完成”导入文本到Excel。
导入后，数据位于A1至A9，就是蛋白质骨架的RMSD值。我们还需要输入时间以进行作图。前面在配置文件中，我们的动力学模拟时间设定是2500 timesteps，1 timestep = 2fs（飞秒），所以共计5000fs = 5ps(皮秒)，又因为我们设定dcd文件的输出频率是250 timesteps输出一次（参见2.3.1配置文件），所以每次输出对应0.5ps。因此在B1-B9中依次输入1, 1.5,…,4.5，5。注意：0.5ps是第一帧，被用作计算的基准帧。
12、为了作图方便，我们将A列的RMSD值数据剪切→粘贴到C列。A列空出即可（图）
13、选择“插入”→“图表”开始作图。在弹出的对话框中选择“XY 散点图”中的右下角子图表“折线散点图”。点击“完成”即可看到做出的初步图形（图）
0.90.80.70.60.50.40.30.20.100123456系列1 14、进一步进行修饰并加上图题等内容可以得到最终结果（图）。
RMSD0.90.80.70.60.50.40.30.20.1000.511.5球形水体中泛素RMSD(5ps)T/ps22.533.544.555.5 其实上面的这张RMSD图中的信息量很少，因为本次分子动力学模拟时，能量平衡进行的时间太短，写入dcd文件的频率又低，2500 timesteps 中，一共输出了10帧，因此作图呈离散点状。下面是一张典型的RMSD值图：
可以看到，RMSD值在最后逐渐趋于恒定，在一定值上下波动。这说明体系已经达到能量平衡。虽然从我们自己做的RMSD图中得不出明显的结论，但是基本分析方法我们已经掌握了。