分类
新手都可掌握的蜡烛图分析法

在XTB进行交易

010-89755166转810

巨大体系的范德华表面静电势图的快速绘制方法

重要提示 : 下文是很早以前写的,如今已经没意义了!如今使用Multiwfn官网上的最新版本,结合下文提到的xtb产生的molden文件作为输入文件,直接用《使用Multiwfn+VMD快速地绘制静电势着色的分子范德华表面图和分子间穿透图》( http://sobereva.com/443 在XTB进行交易 )介绍的ESPiso脚本绘制才是最好的做法。因为后来得益于《Multiwfn使用的高效的静电势算法的介绍文章已于PCCP期刊发表!》( http://sobereva.com/614 )中介绍的Multiwfn计算静电势速度的巨大提升,以及《巨幅降低Multiwfn结合VMD绘制分子表面静电势图耗时的一个关键技巧》( http://sobereva.com/602 )所描述的巨幅降低耗时的技巧,用443号文中的做法比下文的做法不仅更省事、省时,还不需要调用收费的Gaussian里的cubegen。

提示:如果你要绘制1000原子以上的特大体系的表面静电势图,应当使用《基于原子电荷极快速绘制超大体系的分子表面静电势图》( http://sobereva.com/639 )里的做法,用Multiwfn基于原子坐标和原子电荷来计算,耗时超级低!

巨大体系的范德华表面静电势图的快速绘制方法

First release: 2019-May-7 Last update: 2020-Apr-26

本文笔者的条件是:CPU为Intel i7-2630QM(4核),操作系统是Win7-64bit,Multiwfn是2019-May-5版,xtb是2019-Apr-16版,Gaussian是G16 A.03版,VMD是1.9.3版,这几个程序里只有Gaussian是收费的。笔者跑xtb时是在VMware装的CentOS 7.4虚拟机下跑的,因为此程序目前只有Linux版,不会装Linux的话看《在VMware 15中安装CentOS 7.6的完整过程视频演示》(http://sobereva.com/454)。

先把那个336原子的大分子的结构文件载入Multiwfn(支持mol、mol2、pdb、xyz、gjf等众多格式),进入主功能100的子功能2,选择相应选项导出为xyz文件,比如叫336.xyz。然后把此文件随便放到一个空目录下,在此目录下运行此命令通过xtb做单点计算:xtb 336.xyz --molden(如果这个结构是未经优化过的,可再加上--opt让xtb优化此结构)。在笔者的4核机子上只花了20秒钟就算完了,在当前目录下出现了molden.input,这是Molden输入文件,可以被Multiwfn读取并做波函数分析。

下面产生此体系的电子密度格点数据。启动Multiwfn,载入molden.input,然后依次输入以下命令:
5 //计算格点数据
1 在XTB进行交易 //计算电子密度
3 //高质量格点
计算过程花费仅不到1分钟,然后选2导出格点数据,此时当前目录下出现了density.cub,将之改名为density1.cub备用。

然后我们要计算静电势格点数据,用Multiwfn调用Gaussian目录下的cubegen速度会很快,但是这要求必须以fch作为输入文件,所以我们要把当前的波函数导出为fch文件。接着在Multiwfn里输入
0 //返回主菜单
100 //其它功能part 1
2 //导出文件
7 //导出fch文件
336.fch //输出文件名
马上当前目录下就有了336.fch。现在退出Multiwfn。

启动Multiwfn,依次输入
336.fch
5 //计算格点数据
12 //计算静电势
3 //高质量格点 在XTB进行交易
现在会看到Multiwfn正在调用cubegen计算静电势。这一步是整个过程最耗时的一步,花了50分钟(而在笔者的2*E5-2696v3 36核服务器上三分多钟就能算完)。然后选2导出格点数据,此时当前目录下出现了totesp.cub,将之改名为ESP1.cub。

然后进入Graphics - Representation,点击对应显示分子结构的rep,把CPK显示方式切换为Licorice,并让键的粗细从默认的0.3改为0.2。然后点击对应Isosurface显示方式的那个rep,选择Trajectory标签页,把Color Scale Data Range范围设成比默认-0.03~0.03更宽的-0.04~0.04使得颜色变化更柔和些。最后图像效果如下,两个视角都给了,可见效果相当令人满意!(可能有人觉得这个图还没一开始看到的酷炫,但此时的图更容易分辨分子表面上不同区域静电势的差异)

如何判断化合物优势构象?

注意这里的G为吉布斯自由能,包括了零点能和有限温下的熵贡献。显然,G越低,比例越高。对于简单分子(且各个构型的熵差异不大),我们可以很容易地通过某种 NCI 的存在与否来得到定性正确的优势构象(比如乙烷的 staggered conformation 是优势构象)。这个亚子应付一下考试倒是没问题,但对于涉及多种 NCI 的、稍复杂一点的分子,朴素的空间想象就很容易抓瞎。比如题主给出的 (1R,2S)-1,2-dibromo-1,2diphenylpropane,苯/溴、苯/甲基、苯/苯、溴/甲基、溴/溴,究竟哪一种主导,苯应该处于什么朝向?

所以正道当然不是背诵特例,而是来动手算一算,即进行构象搜索。虽然猛一听似乎穷举+算能量即可,但由于复杂分子的自由度非常多,扫描整个PES不现实,所以这是一个 non-trivial task,本质是一个全局优化问题。现在领域内的解决方法有模拟退火、遗传算法、粒子群优化、basin hopping等等等等。本文中我们通过分子动力学(molecular dynamics,MD)模拟,以高温下的轨迹作为取样。

结果与讨论

首先搭建一个(1R,2S)-1,2-dibromo-1,2diphenylpropane 的初始结构~

  • 用 xtb 在 GFN0-xtB 下做一个简单的几何优化,拿优化后的结构去跑一个 MD:500 K,100 ps,步长 2 fs,约束氢原子的化学键,每 100 fs 保存一帧。16分钟跑完。
  • 用 xtb 在XTB进行交易 在 GFN0-xtB 下对轨迹(总共1000帧)进行半经验粗优化,每个点大约 3 s,50 分钟跑完。用 molclus 对粗优后的构象进行归类(阈值:能量 0.25 kcal/mol,结构 0.25 Å),筛除重复。
  • 在 GFN2-xtB 下对分类和去重后的构象(总共91帧)进行半经验精优化,10分钟以内跑完。用 molclus 对精优后的构象进行归类(阈值:能量 0.25 kcal/mol,结构 0.25 Å),筛除重复。
  • 用 Gaussian 在 B3LYP-D3(BJ)/6-31G* 下对距离全局最优点 4 kcal/mol 以内的构象进行优化与频率计算,接着在 M06-2X-D3/6-311+G** 下算稍高精度单点;每个构象耗时20分钟左右。用 molclus 读取优化后的结构、单点能和自由能热矫正,进行归类(阈值:能量 0.25 kcal/mol,结构 0.25 Å),筛除重复。最终只保留了 3 种构象:

其实也就是 C1-C2 在XTB进行交易 旋转的三种 staggered conformation,记作 A、B、C。这时候我们用得到的单点能和热矫正项计算出吉布斯自由能,通过玻耳兹曼分布计算出 A、B、C 在 298.15 K 这三种构象的比例分别为 34.86%、46.94%、18.21%,即 B > 在XTB进行交易 A > C。可是 B 与 C 的结构明明更接近呀,为甚么插进去一个第三者 在XTB进行交易 路人A?于是我们又算了个 B3LYP-D3(BJ)/6-31G* 下的二面角扫描(16核跑了1h,但其实没必要):

Y 轴单位 kcal/mol,X 轴单位°

可见,构象转换的能垒均在 10 kcal/mol 左右(小于20,室温可以跨越),而且 DFT 能量的顺序为 C > B > A。这样看来 A才应该是优势构象呀?为了更深刻理解结构中存在的 NCI 以及对能量的贡献,我们用读取波函数计算了 sign(λ_2)ρ 与 RDG 函数,做出了散点图,并画出了空间中0.在XTB进行交易 5的等值面:

  1. A 结构中两个苯环处于相离最远的位置,除了同碳原子上的溴/苯基氢位阻之外,还存在有较弱的C1、C2取代基之间的溴/苯、甲基/苯、氢/苯位阻互斥。
  2. B 的散点图上可以明显看出相比 A 更多的吸引和互斥作用。左起第一条虚线对应的是两个苯基之间的π-π堆积,但是 C1、C2取代基之间的溴/溴、甲基/溴、溴/苯位阻互斥却大大增强了,形成了一个三叶草形状的较大等值面。
  3. C 结构中,C2上甲基与C1上苯基的位阻互斥导致两个苯基之间的二面角有 65.在XTB进行交易 4°,而 B 的只有54°。因此,C的两个苯基之间的π-π堆积被削弱了,散点图上显影的spike也减弱到了左起第二条虚线处。在散点图的右侧,也出现了一个较大的甲基/苯位阻排斥作用(右起第一条虚线)。

通过RDG分析,我们理解了电子能的 C > B > A 趋势。加上 ZPE 之后,该趋势依然保持。那为什么考虑了室温下的振动熵之后,吉布斯自由能的顺序变为 C 在XTB进行交易 > A > B?三种结构的差异来自哪里?一定是来自于振动模式,于是我们读取了之前频率计算的结果:

B3LYP-D3(BJ)/6-31G*

A 在 612 左右的峰的 amplitude 远强于 B 和 C,而在振动熵(通过简谐近似)的计算时,他将会有较大的贡献,导致考虑热矫正之后 A 变得比 B 更不稳定了。该振动模式长这样:

在 A 中,两个苯基的氢和C1、C2可以不太受阻碍地垂直于苯平面摆动。而在 B 之中,较强的π-π堆积作用限制了两个苯基的活动空间,使得该振动的被大为削弱。C 中苯基之间的 π-π堆积也存在,但是相比A稍弱,所以热矫正项的值也稍大于B的,振动谱中对应的峰的高度也略低。

软件在哪里下载?

点赞 点赞

电 话:021-34750366

邮 箱:[email protected]

地 址:上海市青浦区徐泾镇诸光路1588号

Copyright 版权所有 ©北京华航唯实机器人科技股份有限公司 京ICP备 13042403号-5

010-89755166转810

Hare80/orca-oniom-accelerator

This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.

Use Git or checkout with SVN using the web URL.

Work fast with our official CLI. Learn more.

Launching GitHub Desktop

If nothing happens, download GitHub Desktop and try again.

Launching GitHub Desktop

If nothing happens, 在XTB进行交易 download GitHub Desktop and try again.

Launching Xcode

If nothing happens, download Xcode and try again.

Launching Visual Studio Code

Your codespace will open once ready.

There was a problem preparing your codespace, please try again.

Latest commit

Git stats

Files

Failed to load latest commit information.

README.md

  • 通过备份xtbrestart文件加速ORCA的ONIOM计算过程
  • 给xtb联用提供更多参数
  • 备份xtb输出文件,方便出错时查找问题

嘉然

[1]该程序可能会对原来文件进行修改,请注意备份
[2]本程序目前支持GFN2-xTB(即含有下列关键词),其他方法可能会有bug

[3]计算中减少时间大约为两次全体系采用GFN2-xTB方法计算单点时间,对于小体系加速效果并不明显(甚至可能降低效率)
[4]最好不要用ORCA算 QM / XTB