首页/文章/ 详情

分子动力学原理与蛋白配体GROMACS建模仿真(1月10日)

精品
作者优秀平台推荐
详细信息
文章亮点
作者优秀
优秀教师/意见领袖/博士学历/特邀专家/独家讲师
平台推荐
内容稀缺
4月前浏览3344
导读:GROMACS是用于研究生物、化学、物理等学科中分子微观机理的分子动力学程序包。它可以用分子动力学、随机动力学或者路径积分方法模拟溶液或晶体中的任意分子,进行分子能量的最小化,分析构象等。
它的模拟程序包含GROMACS力场、AMBER力场等诸多力场,研究的范围可以包括玻璃和液晶、聚合物、晶体、界面体系和生物分子溶液(蛋白质、核苷酸、糖等)。GROMACS是一个功能强大的分子动力学的模拟软件,其在模拟大量分子系统的牛顿运动方面具有极大的优势。
自2023年起,笔者原创分子动力学模拟GROMACS系列课程都已经上架仿真秀平台,后续计划更新更多新的课程,还可以为学员提供咨询服务等。为了帮助读者朋友更好的学习GROMACS分子模拟1月10日20时,笔者受邀在仿真秀2024设计仿真新年报告会第二期给读者朋友带来《分子动力学原理与蛋白配体建模仿真》公开课,这也是笔者在仿真秀的首次直播,详情见后文,欢迎各位老师,同学和工程师朋友前来交流互动。
一、分子动力学:微观世界的经典力学方法

分子动力学(Molecular Dynamics,简称MD),是从经典物理的统计力学出发的计算方法,它通过对分子间相互作用势函数及运动方程的求解,分析其分子运动的行为规律,模拟体系的动力学演化过程,给出微观量(如:分子的坐标与速度等)与宏观可观测量(如:体系的温度、压强、热容等)之间的关系,从而研究复合体系的平衡态性质和力学性质,是研究材料内部流体行为、通道运输等现象有效的研究手段。

分子动力学的基本思想是:根据玻恩-奥本海默近似,可以将原子核与电子分离开,再将原子核假设为组成体系的经典粒子进行研究。按照经典力学牛顿定律,体系中的经典粒子受力为:

其中,F为粒子所受的力,a为粒子的加速度,而m是粒子的质量。我们知道,力可以表示为势能函数对坐标的一阶导数,加速度可以表示为速度对坐标的一阶导数,我们将这组方程改写为以下两组方程:

其中,v是粒子的速度,U为相应的势能函数。显然,解这两组方程是非常繁杂的,我们通常会采用数值解进行求解,因为严格的解析解只可能出现在简单的势函数作用的体系之下,这在实际的研究中是几乎不可能出现的。有理论力学基础的同学应该对这个解法不陌生,它与解相空间的哈密顿方程有点相似。这样一个随时间演化的体系, 我们可以根据坐标描述其运动轨迹,根据速度描述其运动轨迹上的变化,从而完全描述其演化过程,这是所谓的经典力学的“确定性”。通过给定势函数,赋予体系初始的坐标和速度,我们会得到一系列包含了整个分子动力学过程的坐标与速度,再通过对坐标与速度的统计,我们可以得到需要的体系相应的热力学与动力学性质,流程草图如下:

通过方程可以看出,势函数(也称之为“力场”)的选取对分子动力学的结果非常重要,加之分子动力学对计算资源的消耗非常大,一个常见的分子动力学程序可能往往要跑一周以上,因此合理选取势函数、正确使用算法也成了分子动力学实际操作中的重点。对于势函数的选取,通常选用一些经典势(如:伦纳德-琼斯势,学过固体物理的同学肯定不会陌生)。当然,势函数也可以通过第一性原理计算,通过电荷密度得到Hellmann-Feynman力得到,这种方法又被称之为“第一性原理分子动力学”,目前应用也非常广泛,此处不再展开。 

二、OPLS-AA力场的聚合物建模及模拟

第一步:OPLS-AA力场简介/高分子模拟常用力场

尽管Gromacs与Amber都包含适合蛋白质与核酸体系模拟的力场,但到目前为止罕有专为人工聚合物体系开发的力场。不同于蛋白质/核算具有有限种类的残基,聚合物原则上可以有无限种不同的单体,去为每一种单体都设定优化参数是不可能的。

因此在实际作业中,往往使用普通的全原子力场对高分子体系进行模拟,包括COMPASS(仅MS)、CVFF(仅MS),OPLS-AA, DREIDING,GAFF,MM3,AMBER等。其中OPLS-AA是作为高分子模拟中比较常用的一种,全称Optimized Potentials for Liquid Simulations,转为准确再现凝聚态物化性质开发。被Gromacs原生支持,并且也有成熟的拓扑文件生成工具(LigParGen,TPPmktop,GMXtop),使用起来较为便捷。

二步:聚合物建模方法

高分子模拟的一大难点在于如何构建具有较高聚合度的分子以及拓扑文件。

对于聚合度低的情况(10-40个重复单元),最简单的方法是直接在GaussView中构建单体判断,重复拼接,保存为pdb格式然后将得到的分子直接丢到TPPmktop得到拓扑文件。

当聚合度较高时,以上方法就不再适用,一方面手动拼接上百个单元过于耗时,另一方面过大的大分子也不被TPPmktop等拓扑文件生成程序支持。此时推荐的方法是添加非标准残基来生成拓扑文件。残基的概念是来自于蛋白质,gromacs处理蛋白质建模的方法是通过pdb2gmx插件读取蛋白质pdb文件记载的残基名,匹配数据库对于蛋白质每个残基成键参数的定义来生成拓扑文件。对于人工高分子也可以用这种方法,通过自定义片段的信息(RTP文件),也就是添加非标准残基,来用pdb2gmx自动生成高分子的拓扑文件。

第三步:Gromacs添加非标准残基建模聚合物的方法

本文所有的文本编辑(RTP,PDB文件等)都默认使用Notepad++,会使用到一些特殊功能,建议安装。

首先在GaussView中构建单体单元,本例中使用的是DMA单元(N,N-DiMethylAcrylamide)

全选原子,点击GaussView窗口菜单中的Custom Fragments(在DNA图标右边),选中一个Group右键点击New Fragment from active molecule,就得到了自定义的基团,然后构建一个三单元的聚合物用于获取头基、尾基及中间单元的RTP文件,保存为PDB格式。这里建议修改原子名称,将同元素原子用编号区别开,例如按照如下格式,记得要保持原有对齐,比原来多了几个字符就在其后删去多少空格(参考压缩文件1.pdb):

C1 H11 H21 H31 C2 H21 C3 O N C4 H41 H42 H43 C5 H51 H52 H53……

(第二个重复单元也可以C1开始编号)

将所得pdb文件提交到TPPmktop ,点击Load PDB file for making all-atom OPLS topology右边的按钮提交文件,得到生成的itp与rtp文件。我们只需要rtp文件。

rtp文件包含了指认的原子类型、成键类型、improper二面角,其它未指认的均使用opls-aa的默认参数。

在这里我们需要依据这些数据构建自己的rtp文件,由于头基包含17个原子,从得到的rtp文件顺序复 制前17个原子信息即可,注意核对电荷总和,对于中性的残基片段应为0;但bonds的顺序是乱的,这里我们就不使用得到的rtp文件,自己按照成键规则编辑,注意一项是,与下一个片段相连的键,需要在下个片段的原子前加上+,表示是下个片段的头原子,在下个片段时则在上个片段相连的原子前加上-号,表示与上个片段的末原子相连;improper项根据这个片段的原子从得到的rtp中复 制,如果出现与下个原子相连的情况,也需要加上+号。

最终效果(这里只列举了其中一个残基,完整见压缩文件DMA.rtp):

将编辑好的rtp文件放到 你的gromcas文件夹/share/gromacs/top/oplsaa.ff内即可使用

第四步:聚合物PDB文件生成

对于构建较高分子量的聚合物,通过脚本复 制重复单元来生成高分子链是可行的方法。编辑之前得到的1.pdb,删去头基尾基除了碳的所有原子(这个在notepad中删除更快点),然后修改头尾碳为氯原子(用于标识),并用GaussView修改C-Cl键长为原来的一半,然后另存为gjf文件(如果gjf文件带有残基名,那么先存为mol2再存为gjf除去残基名),用notepad打开gjf文件,只复 制原子坐标信息,粘贴到新文件中,在首行插入片段的原子数,再留一空行,保存为1.xyz,然后使用以上链接的脚本(或者压缩包的polymerize.sh)生成聚合物xyz文件。

直接生成的聚合物存在几何构型不合理的问题,往往会被VMD判定存在不合理的成键,但在能量最小化中可以解决。

要使用pdb2gmx生成拓扑文件,我们需要得到带有残基号的PDB文件。所以下一步是用得到的xyz文件得到带有残基名的pdb文件。

pdb文件有一点麻烦的是文件里面的每个记录都有着严格的格式。每个记录中的字段, 如标识, 原子名称, 原子序号, 残基名称, 残基序号等, 不仅要按照严格的顺序书写, 而且每个字段所占的字符串长度, 及其所处的位置都是严格规定好的。在这里我使用自己编辑的excel文档xyz2pdb.xlsx与notepad来转换xyz格式并确保得到的pdb文件格式合法

如图所示,在P列左端的为工作区,用于记载原始数据,黑粗闸线框住的是输出区

具体的使用方法如下:

首先另开表格导入高分子的xyz文件(包含原子名称,坐标),将原子名称一列复 制到O列,坐标复 制到KLM三列。C列是原子序号。E列的原子名参考RTP的命名方式。G列用于修改残基名,自己定义。I列是判断残基序号的,你需要根据自己体系的残基原子数来修改I列,这里我在MID残基处使用=IF(I18<10," ",(IF(I18<100," "," ")))函数来判断残基序号。

其它空列(BDFHJ)并非无用,包含了为了对齐pdb格式所用的空格及函数,根据自己需要修改(一般不需要动)。

修改完后,直接复 制整个黑粗闸线区的内容到notepad上并转换tab为空格(在编辑>空白字符操作一栏,这步很重要)即可得到标准pdb格式的文件。

第五步:开始模拟

做好以上工作后,就可以使用pdb2gmx指令得到top拓扑文件与gro坐标文件;不过由于指认二面角的函数类型有误,所得到的拓扑文件还需要经过修改。

打开拓扑文件,修改[ dihedrals ]项的函数类型为“3”(第五个个数字),improper[ dihedrals ]项的函数类型为“4”,这个操作在notepad上批量修改较为方便(例如替换“1 i”为“4 i”,含空格)。

gro文件可能也需要修改,当聚合度很高时,某个坐标数值也会变得很大,自动创建的模拟盒子大小也会不足够,这时修改gro文件最后一行为一个合适的数值(三个数字均为最长值的两倍),进行能量最小化与真空NVT模拟后,笔直的高分子就会蜷缩为一团,这时候再调整模拟盒子大小为你期望的值,然后用solvatebox 指令添加溶剂。

三、蛋白配体GROMACS建模仿真公开课    

为了帮助读者朋友更好的掌握GROMACS分子模拟能力, 1月10日20时, 仿真秀2024设计仿真新年报告会  第二期,我们邀请仿真秀平台优秀讲师,北京大学博士后刘十三老师做分子动力学原理与蛋白配体建模仿真公开课,欢迎各位老师,同学和工程师朋友前来交流互动。

以下是直播安排

《分子动力学原理与蛋白配体建模仿真》

扫码观看和回放

(完)

来源:仿真秀App
ACT化学电子理论材料分子动力学GROMACS
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-01-11
最近编辑:4月前
仿真圈
技术圈粉 知识付费 学习强国
获赞 9252粉丝 20533文章 3207课程 207
点赞
收藏
未登录
还没有评论

课程
培训
服务
行家

VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈