NAMD入门教程(一) (2)

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 使用 这些文件包括 蛋白质分子的蛋白质分子的 PDBPDB 文件 文件 该文件负责储存蛋白质中所有原子的坐标 在后续课程 中我们还会了解到 PDB 文件还可以储存原子运动的速度等信息 蛋白质分子的蛋白质分子的 PSFPSF 文件 文件 该文件负责储存蛋白质的结构信息 注意 PDB 文件只记 录原子的空间位置 并不储存蛋白质中原子之间的成键情况 成键情况由 PSF 文 件负责记录 力场参数文件力场参数文件 force force fieldfield file file 力场参数文件是分子动力学模拟的核心 文 件中的数学方程决定了原子在力场中的受力如何计算 常用的四种力场是 CHAEMM X PLOR AMBER 和 GROMACS NAMD 可以使用以上任何一种力场进行分子动力学模 拟 配置文件 配置文件 configurationconfiguration filefile 配置文件的目的是告知 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 每输入一行命令后按回车键 下同 另外 尤其要注意空格的有无和空格的位置 否则空格位置不对可能造成命令执行错误 提示 提示 VMDVMD TKTK ConsoleConsole VMDVMD 控制台 中改变当前目录的方法控制台 中改变当前目录的方法 在 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文件 该文件包含蛋白结构的全部信息 知识链接 组氨酸的三种离子模式知识链接 组氨酸的三种离子模式 知识链接 知识链接 PDBPDB文件中原子的命名方式文件中原子的命名方式 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 分子速 率分布给各原子赋予一定的初始速度 接下来就是要根据牛顿力学方程 求解个水分子 以及蛋白质中各原子的运动轨迹 我们得到的结果 就是模拟泛素这一小肽在溶液状态 下的运动状态 知识链接 能量最小化和能量平衡 知识链接 能量最小化和能量平衡 MinimizationMinimization andand EquilibrationEquilibration 本次动力学模拟实际包括两个过程 能量最小化和能量平衡 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 说明所有和氢原子相连的键都被认为是不振动的 知识链接 知识链接 RigidRigid BondsBonds 为什么要设定 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 为了使动力学模拟能够顺利完成 在 此期间不要进行任何操作 提示 提示 WindowsWindows 命令行中改变当前目录的方法命令行中改变当前目录的方法 在 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 横截面积不变