Hsiun YuBin's Blog

熊煜彬的个人博客

0%

Gaussian + orca + xtb + molclus 对阿司匹林做构象搜索

发表于 更新于 分类于
本文字数： 6.8k 阅读时长 ≈ 6 分钟
本文将使用计算化学程序 Gaussian、orca、xtb 以及 molclus 对阿司匹林做构象搜索。

前言

在计算化学中，构象搜索是研究一个体系经常涉及到的问题，特别是有机分子的构象。我们都知道，很多有机体系中都存在着可以任意旋转的键，一个简单的乙烷都有两种构象，更别说那种更加复杂的有机物了。为了寻找一个体系的能量最低构象，通常需要依靠计算化学的手段。本文将介绍笔者最愿意使用的一种方法，对常见有机物阿司匹林进行了构象搜索（阿司匹林的结构如下图所示）。

对一个有机物做构象搜索有非常多的方法，在计算化学公社上有很多例子。最原始的方法莫过于使用 Gaussian 对体系进行刚性或柔性扫描，这种方法在 Sob 老师《详谈使用 Gaussian 做势能面扫描》一文中有详细的介绍，感兴趣的可以浏览，本文不介绍这种方法。本文主要使用的是由 Sob 老师开发的免费构象搜索程序 molclus 对体系进行构象搜索。根据 Sob 老师的帖子，molclus 主要有以下几种使用方式：

  • genmer 结合 molclus：用于团簇分子或原子团簇构型搜索。使用简便
  • gentor 结合 molclus：用于分子构象搜索。使用简便，但不支持环状区域构象搜索（环状区域是指诸如环己烷这样有多种构象的环状区域）
  • Confab 或 Frog2 构象生成程序结合 molclus：用于分子构象搜索。使用最为傻瓜化，但不支持环状区域构象搜索、不支持普通有机分子以外的情况
  • 经典力场分子动力学程序结合 molclus：用于分子团簇搜索和分子构象搜索。使用稍繁琐，被动力学模拟的体系必须有恰当的力场
  • xtb 程序跑动力学结合 molclus：普适性极强，用于构象搜索、原子/分子团簇构型搜索皆可，使用容易（不受力场约束，不需要准备拓扑文件）

本文使用的方式为最后一种，xtb 程序跑动力学结合 molclus 进行构象搜索。本文中所有使用的程序以及版本如下所示，所使用的计算环境为 Rocky Linux 8.8 64bit。

  • Gaussian 16 A.03 Linux
  • orca 5.0.4 Linux
  • Molclus 1.9.9.9 Linux
  • xtb 6.6.0 Linux
  • Shermo 2.3.6
  • Multiwfn 3.8(dev) Linux

计算流程说明

根据 Sob 老师的帖子，使用 xtb 程序结合 molclus 进行构象搜索的流程应该为以下流程：

  • 用 xtb 在 GFN0-xTB 下跑动力学
  • 通过 Molclus 调用 xtb，在 GFN0-xTB 下对动力学的每一帧进行批量优化
  • 通过 Molclus 调用 xtb，对第二步优化产生的每一个结构用 GFN2-xTB 在溶剂模型下进行批量优化
  • 通过 Molclus 调用 Gaussian，对第三步筛选出来的结构用 B3LYP-D3(BJ)/6-31G* 在 PCM 的水模型下做优化和振动分析得到较可靠的结构和自由能热校正量，再调用 orca 对每个结构在 PWPB95-D3(BJ)/def2-TZVPP 结合 SMD 模型表现的水环境下算高精度单点能，二者相加得到水环境下的高精度自由能。

用 xtb 在 GFN0-xTB 下做动力学模拟

首先创建一个含有 xtb 的 MD 任务设置的文件 md.inp，内容如下，双斜杠后面的是注释。

1
2
3
4
5
6
7
8
$md
   temp= 400  //温度(K)
   time= 100.0  //模拟的总时间(ps)
   dump= 50.0  //每多少fs往轨迹文件里写入一次
   step=  1.0  //步长(fs)
   hmass=1  //氢原子的质量是实际的多少倍
   shake=1  //将与氢有关的化学键距离都用SHAKE算法约束住
$end

确保 xtb 已正常安装了，然后运行命令：

1
xtb Aspirin.xyz --input md.inp --omd --gfn 0

其中 --omd 告诉程序先做几何优化，然后再跑 MD。做几何优化是为了避免由于初始结构太烂，导致一开始某些原子受到过大的斥力使其速度太大而令动力学崩溃。--gfn 0 代表用 GFN0-xTB。计算完成后一共能得到 2000 个帧，这 2000 帧的结构都在 xtb.trj 文件中，之后将程序生成的 xtb.trj 改名为 traj.xyz，用于下一阶段的计算。Aspirin.xyz 是一种记载阿司匹林各原子坐标的文件格式，可以使用 Multiwfn 生成，具体的生成方法可以浏览《谈谈记录化学体系结构的xyz文件》，本例的 Aspirin.xyz 如下所示。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
    21
 Aspirin
O      -0.86120000     -1.22620000     -0.26620000
O       3.06170000     -1.68520000      0.19610000
O       1.44470000     -2.56350000     -0.69000000
O      -2.41890000     -2.59030000     -0.07290000
C      -0.14120000     -0.17850000     -0.41320000
C       1.27840000     -0.29560000     -0.34330000
C      -0.73020000      1.09250000     -0.62320000
C       2.04150000      0.90200000     -0.34200000
C       0.04910000      2.25900000     -0.64080000
C       1.43780000      2.16480000     -0.47240000
C       1.93690000     -1.52710000     -0.27010000
C      -1.97000000     -1.49550000      0.23250000
C      -2.73710000     -0.68480000      1.15950000
H      -1.73960000      1.18500000     -0.75870000
H       3.06100000      0.88610000     -0.25320000
H      -0.39460000      3.17150000     -0.76840000
H       2.00640000      3.01490000     -0.45520000
H      -3.02830000     -1.26930000      2.03410000
H      -3.63710000     -0.31380000      0.66640000
H      -2.17310000      0.17070000      1.53410000
H       3.51360000     -1.01700000      0.54670000

用 Molclus 调用 xtb 在 GFN0-xTB 下做批量优化

从 Molclus 官网下载最新版本并解压。把 traj.xyz 放到 Molclus 目录下。把 Molclus 目录下 settings.ini 里的 iprog 设为 4 代表调用 xtb，确保 itask 为 0 代表是做优化任务。把 xtb_arg 参数设为 --gfn 0 --chrg 0 --uhf 0，代表这是中性单重态体系且用 GFN0-xTB 计算。之后输入 ./molclus 命令启动当前目录下的 Molclus 程序可执行文件，程序就开始调用 xtb 优化 traj.xyz 里的每一帧了。

任务结束后，优化后的每一帧的结构和能量就都存到当前目录下的 isomers.xyz 里了。使用 isostat 对 isomers.xyz 里的结构去重和做能量排序，即输入：

1
2
3
4
5
./isostat
[Enter]     // 处理当前目录下的 isomers.xyz
[Enter]     // 使用默认的能量去重阈值(kcal/mol)
[Enter]     // 使用默认的结构去重阈值(埃)
[Enter]     // 不计算Boltzmann分布比例

最后一共得到 20 多个非重复结构产生在了当前目录下的 cluster.xyz 里。删掉之前的 traj.xyz，然后把 cluster.xyz 改名为 traj.xyz，使得这些去重后的结构作为接下来 Molclus 任务的输入结构，在此之前可以将 cluster.xyz 拖入 VMD 里浏览产生的构象。

用 Molclus 调用 xtb 在 GFN2-xTB 结合隐式水模型下做批量优化

接着把 settings.ini 里的 xtb_arg 参数设为 --gfn 2 --gbsa h2o --chrg 0 --uhf 0，代表用 GFN2-xTB 级别结合 GBSA 模型表现的水环境进行计算。启动 molclus 来调用 xtb 对当前所有结构进行优化。之后再次运行 isostat 进行去重和排序，做法同前。

1
2
3
4
5
./isostat
[Enter]     // 处理当前目录下的 isomers.xyz
0.5         // 能量去重阈值(kcal/mol)
0.25        // 结构去重阈值(埃)
[Enter]     // 不计算Boltzmann分布比例

程序运行后，一共得到了 7 个非重复的构象，同时载入 VMD 中观看这些构象。

用 Molclus 调用 Gaussian 和 orca 得到水环境下的高精度自由能

把 settings.ini 里的 iprog 设为 1 代表调用 Gaussian。itask 设为 3 代表做优化+振动分析来得到自由能。把 gaussian_path 和 orca_path 分别设为调用 Gaussian 和 orca 的实际命令。同时把 Shermo_path 修改为 Shermo 的具体路径，Shermo 将被用在最后计算溶解自由能。另外 Sob 老师还建议把 ibkout 设为 1 代表把算每个体系的输出文件都进行备份，便于之后必要时进行检查。

如果你认真了解过 molclus 的话，那么就知道调用 Gaussian 并创建自定义任务时，需要修改模板文件 template.gjf，以下是本例中的模板文件内容。

1
2
3
4
5
6
# opt freq b3lyp/6-31g* int=fine em=gd3bj scrf(solvent=water)

Template file

0 1
[GEOMETRY]

同时，在 molclus 目录下创建一个叫 template_SP.inp 的 orca 输入文件的模板文件，内容如下，对应在RI-PWPB95-D3(BJ)/def2-TZVPP 结合 SMD 模型表现的水环境下做高精度单点计算。

1
2
3
4
5
6
7
8
9
10
! PWPB95 D3 def2-TZVPP def2/J def2-TZVPP/C RIJCOSX tightSCF noautostart miniprint nopop
%maxcore 15000
%pal nprocs 8 end
%cpcm
smd true
SMDsolvent "water"
end
* xyz 0 1
[GEOMETRY]
*

在 Molclus 每次调用 Gaussian 做完优化和振动分析任务后，当程序发现当前目录下有 template_SP.inp 文件时，就会自动调用ORCA做单点计算，并将所得的高精度单点能自动加到振动分析得到的自由能热校正量上作为此结构的最终能量，如果将这个能量再加上 1.89 kcal/mol，就正是该溶剂下的自由能。详情请浏览《谈谈隐式溶剂模型下溶解自由能和体系自由能的计算》。

运行 Molclus，当每个结构振动分析完成后 Molclus 会在屏幕上提示有没有虚频，在 isomers.xyz 文件里每个结构的Nimag = 字样后面也可以看到虚频数目。当检查无虚频后，还是用 isostat 处理得到的 isomers.xyz，并将 298.15 K 下的 Boltzmann 分布一同输出。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
./isostat
[Enter]     // 处理当前目录下的 isomers.xyz
[Enter]     // 使用默认的能量去重阈值(kcal/mol)
[Enter]     // 使用默认的结构去重阈值(埃)

 #    1 Count:    2  E=   -648.424923 a.u.  DGmin=   0.20  DE=    0.00 kcal/mol
 #    2 Count:    1  E=   -648.424338 a.u.  DGmin=   0.18  DE=    0.37 kcal/mol
 #    3 Count:    1  E=   -648.423848 a.u.  DGmin=   0.18  DE=    0.67 kcal/mol
 #    4 Count:    1  E=   -648.421557 a.u.  DGmin=   0.19  DE=    2.11 kcal/mol
 #    5 Count:    1  E=   -648.420735 a.u.  DGmin=   1.15  DE=    2.63 kcal/mol

298.15      // 在 298.15 K 下计算 Boltzmann 分布
  
 #    1   Count:   2   Ratio:  52.67%
 #    2   Count:   1   Ratio:  28.35%
 #    3   Count:   1   Ratio:  16.87%
 #    4   Count:   1   Ratio:   1.49%
 #    5   Count:   1   Ratio:   0.62%

这里将 Boltzmann 分布中分布比例前三的三个构象单独拿出来，分别绘制 IRI 等值面图，来考察分子内部的相互作用，如果不清楚 IRI 是什么的可以阅读笔者之前写的博文《Multiwfn + VMD 绘制 IRI 以及 IRI-pi 等值面图》。从下图可以清晰的看出，分布最多的那种构象，存在着 O-H···O 氢键，这必定是这种构象能量最低的原因所在。如果对分布最少的构象感兴趣，可以自行绘制后两种构象的 IRI 图，在这里就不再多赘述了。

请注意！如果你要考察其它温度下的 Boltzmann 分布，不仅 isostat 里要输入相应的温度，在 template.gjf 里也应当加上 temperature=[温度] 关键词来让 Gaussian 输出的自由能热校正量对应相应的温度。

参考文章

本文所介绍的构象搜索流程，主要参考了计算化学公社社长 Sob 老师的帖子：

除了以上有关的帖子外，由于分子动力学不在笔者的研究范围之内，所以未提及使用分子动力学程序搭配 Molclus 和 Gaussian 进行团簇构象和分子构型的搜索。感兴趣的同学可以自行浏览以下帖子：

  • 本文作者： Hsiun YuBin
  • 本文链接： https://ikuns.icu/012/
  • 版权声明： 本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！