| 查看: 5365 | 回复: 14 | ||||||||||||
| 【奖励】 本帖被评价13次，作者sobereva增加金币 10.4 个 | ||||||||||||
| 本帖产生 1 个 QC强帖 ，点击这里进行查看 | ||||||||||||
[资源]
使用Multiwfn研究分子动力学中的弱相互作用
|
||||||||||||
|
使用Multiwfn研究分子动力学中的弱相互作用 1 前言 曾经笔者在《使用Multiwfn图形化研究弱相互作用》（http://hi.baidu.com/sobereva/item/bb1162d881b8d5ee3cc2cbdc）中详细介绍了杨伟涛课题组在JACS,132,6498 (2010)提出的研究弱相互作用的约化密度梯度(RDG)方法。原先的RDG方法只考虑单个结构下的弱相互作用，然而由于热运动，实际体系的结构是不断波动的，其中的弱相互作用也因此不可能只通过单一的结构全面表现出来。在JCTC,9,2226 (2013)中，杨伟涛课题组提出了averaged RDG (aRDG)方法，它是原先RDG方法的扩展，使得RDG分析方法可以用于分析多帧结构，尤为适合结合分子动力学模拟技术来研究平衡的动态环境中的弱相互作用。虽然aRDG原则上说也可以结合诸如蒙特卡罗模拟，但下文讨论的只限于aRDG方法对全原子分子动力学轨迹的分析。 2 原理 如果不熟悉RDG方法，请务必先阅读《使用Multiwfn图形化研究弱相互作用》。aRDG方法和原先的RDG方法在原理上十分类似，同样是通过RDG函数的等值面来展现弱相互作用区域，并且将sign(lambda2)rho函数用不同的颜色投影到RDG等值面上来展现出弱相互作用类型。但是在aRDG方法中，计算RDG和sign(lambda2)rho函数所用的电子密度、电子密度的梯度以及电子密度的Hessian矩阵都是来自对多帧结构取平均得到的。也就是说，比如要研究的分子动力学轨迹有1000帧，那么aRDG所用的密度就是对这1000帧都依次计算密度然后取平均得到。因此aRDG中的RDG函数和sign(lambda2)rho函数都表现的是整条轨迹中平均的RDG和平均的sign(lambda2)rho。 除了这点区别外，aRDG方法还定义了一个热波动指数(thermal fluctuation index, TFI)，TFI=std(rho)/avg(rho)。其中avg(rho)就是上面说的整个轨迹中的平均密度，std(rho)就是整个轨迹中密度的方差，一个函数的方差越大就表明函数的波动越强。热波动指数的用处是展现不同区域弱相互作用的稳定程度，分析的时候将热波动指数通过不同颜色映射到平均的RDG等值面上，通过考察颜色，就能判断不同区域弱相互作用形成的稳定与否。 3 在Multiwfn中的使用 从Multiwfn 3.1版开始，aRDG方法可以在通过Multiwfn的主功能200里的功能1来实现。此方法的原文是通过在NCIplot程序上进行修改得到的，然而在本文撰文时，目前公开的NCIplot程序还未能实现aRDG分析。而即便是原先的RDG分析方法，相对于NCIplot，Multiwfn在功能、易用性以及速度上也都有非常显著的优势。Multiwfn 3.1版可在http://multiwfn.codeplex.com上免费下载。 由于aRDG分析用的是平均的密度及其梯度和Hessian，所以分析结果会随着考虑的帧数增加而逐渐趋于收敛。考虑的帧数越多，结果越精确，图像也越细腻。如果只考虑一帧的结构，那么和原先的RDG分析就没区别了。建议至少考虑500帧结构，如果想作出漂亮的图，应该用1000帧或更多。 使用Multiwfn的aRDG分析功能时，需要在一开始载入.xyz轨迹文件。.xyz格式很简单，.xyz轨迹文件只不过是简单地将一般的只记录单个结构的.xyz文件合并在一起而已，所以算是个比较通用的记录轨迹的格式。用户可以用任何自己最常用的分子动力学程序来跑出轨迹，然后载入VMD程序，之后直接转换为.xyz格式轨迹文件。 在轨迹中，感兴趣的分子的坐标应当保持固定不变，这可以在动力学程序中将其坐标冻结来实现。并且感兴趣的分子最好是在整个体系的靠中间的位置。计算格点数据时设定的计算区域应当能够容纳这个感兴趣的分子，且四周有一定延展区域，以避免重要区域的RDG等值面被截断。 由于需要考虑几百、上千帧，所以计算那么多帧的电子密度、梯度及Hessian的格点数据是不可能直接基于波函数的，因为无论是计算那么多次波函数还是计算那么多格点数据都是极其耗时的。aRDG原文用的是QM/MM方法，感兴趣的区域靠QM来描述并基于波函数获得密度，而外环境部分靠分子力学描述。但QM/MM动力学在实际操作上并不方便，对很多用户来说都有不小门槛，而且依然较耗时，且输出的格式与动力学程序的依赖性也太强因此不够普适。所以，Multiwfn中的aRDG分析一律使用promolecular近似，也就是基于自由原子密度的叠加来近似产生电子密度及其梯度和Hessian。这个近似显著简化了计算过程，大大节约了计算量，而且分析结果依然比较可靠。 4 实例：标况下纯水体系中水分子之间的平均弱相互作用 限于时间和精力，此文只给出一个应用aRDG的例子，也就是分析常温常压下的纯水体系中水分子之间的弱相互作用。对于其它体系的研究请根据本例举一反三。本节涉及到的所有文件包括轨迹都可以在此处下载：http://pan.baidu.com/share/link? ... p;amp;uk=1074012119。 4.1 生成轨迹 用任何程序来获得动力学轨迹都是可以的，本文用gromacs 4.5版来获得。如果读者习惯使用其它动力学程序，下文的gromacs的具体指令就不用管了，只要弄明白每一步的目的就行了。 首先建立一个空盒子文件emptybox.gro，盒子各边长都为2.5nm，内容如下： WAT 0 2.5 2.5 2.5 盒子边长当然可以更长，只不过计算要更费时间。2.5nm算是底限了，再小的话就难以保证盒子长度是非键相互作用cut-off距离的两倍以上了。 然后运行以下命令往盒子里面填充水得到纯水盒子，其中将含有511个水分子 genbox -cp emptybox.gro -cs spc216.gro -o water.gro 运行此命令，并且选GROMOS96 53a6 force field，得到水盒子体系的top文件。水分子用的是SPC/E模型。 pdb2gmx -f water.gro -o water.gro -p water.top -water spce 接下来跑100ps NPT动力学，令水盒子在298.15K 1atm下平衡。范德华作用的cut-off为1.2nm。静电相互作用用PME方法计算，实空间cut-off距离为1.2nm。步长0.001ps。 grompp -f pr.mdp -c water.gro -p water.top -o water-pr.tpr mdrun -v -deffnm water-pr 现在用VMD打开water-pr.gro，从中选一个比较靠近盒子中央的水分子，我们选取101号水分子，下图中将它用vdW方式突出显示  注意101号水分子的三个原子编号是301、302和303，氧原子是301号原子。 这个水分子的坐标在接下来的动力学过程中将保持冻结。为此，先生成index文件，输入以下命令 make_ndx -f water-pr.gro ri 101 // 之所以是ri而不是r是因为ri是从1开始编号的，r是从0开始编号的 q 得到的index.ndx里的group r_101便是101号水分子。 运行以下命令来做298.15K下的1ns的平衡动力学，每隔1ps保存一帧结构，共得到1000帧。注意用NVT系综，因为如果允许控压的话，会对盒子进行scale从而影响101号水的坐标。由于之前已经在1atm下平衡了，水的可压缩性又很低，所以改用NVT模拟依然合理。101号水的坐标靠freezegrps = r_101关键词来冻结住。 grompp -f md.mdp -c water-pr.gro -p water.top -o water-md.tpr -n index.ndx mdrun -v -deffnm water-md 将water-pr.gro和刚生成的water-md.xtc载入VMD，核对一下是否101号水确实固定住了。然后在VMD main窗口中点击与这个体系对应的项，并在其上点击右键，选Save coordinate，将File type设为xyz，Selected atoms选all，Frames框当中的First和Last分别写1和1000，然后点Save按钮，将轨迹保存为wat.xyz。 4.2 用Multiwfn生成格点数据 启动Multiwfn，输入以下命令 wat.xyz 200 // 主功能200 1 // aRDG分析 1,1000 // 分析的帧号范围是1~1000 7 // 通过原子编号自定义盒子中心，并自定义各方向盒子尺寸和格点数 301,301 // 输入两个原子编号，它们的中点即为盒子中心。这里输入的两个原子编号都是301（被固定的水的氧原子编号），也就是说将301号原子作为盒子中心 80,80,80 // 三个方向格点数都设为80 4.5,4.5,4.5 // 盒子在xyz方向的前后延展距离都为4.5埃，因此盒子尺寸是9*9*9埃。盒子设太大的话会使格点间距较大，导致图像粗糙；如果太小，则会截断RDG等值面。所以设置的大小应该根据实际体系凭借直觉来定，如果不确定的话可以先试几个值，且每次都只考虑很少的帧数（如前20帧）来节约时间 之后，Multiwfn开始计算每一帧的电子密度及其梯度和Hessian，并取平均，之后用平均量计算出平均的RDG和sign(lambda2)rho。整个过程需要花一阵子时间，在目前一般的四核机子上耗时约半个小时到一个小时。如果是更大的动力学体系则更费时间，所以可以先暂时干别的去。 计算完毕后可以选选项1来查看平均RDG vs 平均sign(lambda2)rho散点图（分析方法和《使用Multiwfn图形化研究弱相互作用》介绍的一致，这里就不谈了），也可以用相应选项将其保存到图形文件中或者把数据点输出到文本文件里。此体系的散点图如下  选择6将平均RDG和平均的sign(lambda2)rho分别导出到当前目录下的avgRDG.cub和avgsl2r.cub。 由于我们要考察弱相互作用的稳定性，所以选7来计算热波动指数并将之输出到当前目录下的thermflu.cub里。计算热波动指数需要重新计算一遍每一帧的密度，所以也是比较耗时的，和计算平均RDG及sign(lambda2)rho过程耗时基本一样。 4.3 分析结果 本文使用的是VMD 1.9来观看结果。在本节开头提到的压缩包中提供了avgRDG.vmd和avgRDG_TFI.vmd，这是VMD作图脚本。将它拷到VMD程序目录下，并且将avgRDG.cub、avgsl2r.cub和thermflu.cub也都拷到VMD程序目录下。 启动VMD，在控制台窗口输入source avgRDG.vmd，就可以显示出平均RDG等值面，isovalue为0.25。同时，平均sign(lambda2)rho函数也通过不同颜色投影到了等值面上。默认情况下所有分子都以CPK模式显示了出来，这给分析101号水分子周围的弱相互作用带来了不便，所以进入Graphics-Representation，选Style为CPK的那一栏，将Selected Atoms里的内容改为101号水分子对应的三个原子，即改为serial 301 302 303然后点回车，此时就只有101号水分子显示出来了。经过视角的调整，结果如下所示（在VMD里点t之后可以平移视角，点r可以恢复旋转视角模式。点c再点击一个原子可以将视角旋转中心设在相应原子上）  期望的结果在图上都展现了，等值面也没有被截断，说明盒子设得合适。不过不幸的是，在格点数据计算区域的边缘，有一大堆零零碎碎乱七八糟的等值面，这和我们感兴趣的弱相互作用无关却十分碍眼。利用Multiwfn强大的格点数据处理功能，可以将这些离特定的一些原子（此例即101号水分子的三个原子）比较远的区域的等值面给消掉。做法是先写一个文本文件，假设叫index.txt，内容如下 3 //原子数 301 302 303 //原子编号 然后启动Multiwfn并输入： avgRDG.cub // 即刚才生成的平均RDG的格点数据 13 // 主功能13，用于处理格点数据 13 // 设定距离特定原子比较远的格点的数值 1.5 // 如果格点距离特定原子的距离超过相应原子范德华半径的1.5倍，则格点的数值将被设为下面用户输入的值。范德华半径倍数的设定是依靠经验的，可以多试几次并检验结果以达到满意效果 100 // 新设定的值为100，这个值是随意的，只要是个明显大于作RDG等值面图时的isovalue的值就可以将这些点屏蔽掉 index.txt // 从这个文件中读入“特定原子”的编号 0 // 将当前已经被改过的格点数据保存到新的cube文件中 avgRDG.cub // 新输出的格点数据文件名 然后将新得到的avgRDG.cub移动到VMD目录下覆盖原先的avgRDG.cub（可以将原先的avgRDG.cub先备份一下）。然后重新用avgRDG.vmd脚本作图并做前述的调整，结果如下所示（给出了两个视角）  这次的图非常清楚。默认的色彩刻度是-0.25~0.25，对应“蓝色-绿色-红色”的变化。越蓝的地方说明静电、氢键效应越强，越红的地方说明位阻效应越明显，而绿色表明相应位置平均密度值较低，对应于很弱的作用，即色散作用。在两个O-H键对着的地方都有两个蓝色小圆片，表现出在动态环境下，O-H能与其它水形成很显著的氢键。图中有个显著的，胖次形状的绿色等值面，这清楚地展现了水分子倾向于以这些方向和其它水形成范德华作用。在氧的上方有一大坨等值面，其中两侧是蓝色，中间偏红，其蓝色区域体现出氧原子倾向于以这样的方位作为氢键受体，而红色则在一定程度体现出实际环境中氧容易以这样的方向和其它的水接触而表现出互斥效应。氧上方的那一大坨等值面可以看出是由两层值面相连构成的，如果将isovalue设低一些，如0.18，那么两层等值面就分开了，如下所示  之所以会有两层等值面原因尚不明，但可能分别与第一配位层、第二配位层分子的作用有联系。最外层的等值面的形状偏离对称性相对显著一些，而且等值面上投影的颜色有些“斑驳”，这说明对于这样靠外区域的弱相互作用，只考虑1000帧还不足以使结果充分收敛。 关掉之前的等值面图，然后运行source avgRDG_TFI.vmd，并作适当调整，结果如下  默认色彩刻度范围是0~1.5，还是对应蓝-绿-红的变化。越偏蓝的地方，说明热波动指数越小，弱相互作用在动力学过程中越稳定；越红的地方热波动指数越大，弱相互作用在动力学过程中就越不稳定。可见O-H对着的小圆片颜色稍微发蓝，这表明水分子靠O-H作为氢键给体与其它水分子构成的氢键在动力学过程中是很稳固的，很不容易被破坏。而胖次形状的等值面整体呈红色，说明水分子间的范德华作用相对于氢键是明显不易维持稳定的，这容易理解，毕竟范德华作用是很弱的。氧上方的等值面的内层部分主体都是绿色，表现出水分子以氧作为受体与其它水分子构成的氢键还算比较稳定。图中最上方的一大片红色等值面表现出这部分弱相互作用很难保持得住，波动性极大，因此将之忽略了也无妨（aRDG原文中干脆直接把这个区域给截掉了，以避免说不清楚）。 最后，我们来看看如果我们只考虑动力学轨迹中的一帧得到的结果会是如何。下面的图显示的是随机选取的一帧得到的RDG等值面填色图  可见，如果只考虑一帧，表现出的弱相互作用和平均的弱相互作用有很大差距。图中只显示了两个氢键区域，和零星的范德华作用区域，从这样的单帧的图上我们显然是没法直接了解到整个动力学过程中这个水分子是怎么和其它水分子作用的。 5 总结与展望 aRDG方法给研究动态体系中的分子间相互作用增加了一个很重要的新途径。以往研究这类问题往往就是计算径向分布函数、计算概率分布、计算原子间距离/角度数值分布之类，是间接表征手段，远没有aRDG方法直观。aRDG不仅可以研究动力学过程中小分子与其它小分子的平均弱相互作用，还可以研究诸如受体-蛋白间、小肽-小肽之间、分子与固体表面间的平均弱相互作用。aRDG原文中还研究了过渡态结构下体系与溶剂的平均弱相互作用。 相对于单一结构，动力学过程中包含了很多可以影响弱相互作用的因素，特别是温度。比如研究不同温度下平均RDG、平均sign(lambda2)rho以及热稳定指数的差异，对于研究温度对弱相互作用的影响十分有益。其它因素，诸如电场、浓度、压力等等因素对平均弱相互作用的影响也可以尝试通过aRDG分析来研究讨论。 |
回复此楼» 收录本帖的淘贴专辑推荐
 第一性原理 | 资源 | 量化软件学习 | SOB 集锦 |
| 高斯/MOPAC | 计算 | 量化软件-Multiwfn | 文献 |
| 我的淘帖 | Roubeijia |
» 猜你喜欢
深圳湾实验室杨奕课题组招聘博士后、助理研究员及软件工程师
已经有18人回复
-
山东省基金：作为项目负责人承担面上项目数量累计不得超过2项
已经有29人回复
-
物理化学论文润色/翻译怎么收费?
已经有293人回复
-
美国夏威夷大学生物胶体与界面实验室诚招全奖博士生和博士后
已经有1人回复
-
再次祈福国自然基金高中！
已经有20人回复
-
fluent带有化学反应的模拟
已经有0人回复
-
拉曼光谱。。。。。。
已经有1人回复
-
求《物理化学学习进阶与考研攻略》pdf文件
已经有0人回复
-
红外光谱出现了一个很宽的峰是怎么回事来？
已经有0人回复
» 本主题相关商家推荐: (我也要在这里推广)
» 本主题相关价值贴推荐，对您同样有帮助:
- Multiwfn 3.0波函数分析程序的意义、功能与用途
已经有38人回复
- 使用Multiwfn做电子密度、ELF、静电势、密度差等函数的盆分析
已经有30人回复
- 衡量芳香性的方法以及在Multiwfn中的计算
已经有12人回复
- 建议multiwfn开发GUI界面 文件打开什么的效率太低了
已经有10人回复
- 使用Multiwfn做电荷分解分析(CDA)、绘制轨道相互作用图
已经有24人回复
- 如何用Multiwfn软件做出P DOS?
已经有4人回复
- 使用Multiwfn绘制原子轨道图形、研究原子壳层结构及相对论效应的影响
已经有17人回复
- Multiwfn的初步学习，越详细越好。
已经有6人回复
- 使用Multiwfn绘制NBO及相关轨道
已经有30人回复
- 使用Multiwfn作电子密度差图
已经有37人回复
- 使用Multiwfn做拓扑分析以及计算孤对电子角度
已经有18人回复
- 回答一些关于Multiwfn的疑问以及未来Multiwfn的发展打算
已经有37人回复
- multiwfn画图
已经有5人回复
- multiwfn使用中的一些问答
已经有5人回复
- 【原创】Multiwfn 2.0 正式发布
已经有35人回复
- 【sobereva个人文集】电子定域性的图形分析
已经有56人回复
2楼2013-05-18 20:28:39
4楼2013-05-18 21:48:18
5楼2013-05-18 22:01:34
8楼2013-05-18 22:34:03
10楼2013-05-19 19:39:23
11楼2013-05-20 17:43:25
12楼2013-05-20 22:24:14
13楼2013-05-21 14:23:30
15楼2018-06-07 18:21:58
简单回复
yjr3楼
2013-05-18 20:44
回复
五星好评 牛！！！！！
czyzsu6楼
2013-05-18 22:02
回复
五星好评 顶！
安德7楼
2013-05-18 22:25
回复
五星好评 
贺仪9楼
2013-05-19 10:29
回复
五星好评 支持！
Quan.14楼
2013-06-20 05:11
回复
五星好评 谢谢分享！！！！！！！！！！！！！
5