摘抄:同时扫描两个坐标的方法

前言

本文综合摘抄自不同文章和帖子,标注在一级标题后。

第一篇摘抄

《在Gaussian16中同时扫描两个反应坐标》

前言

本博客之前摘抄过同样来自《量子化学》公众号的,高斯中的两种常见势能面扫描到文章:

用高斯做势能面扫描(一):刚性扫描

用高斯做势能面扫描(二):柔性扫描

可能大家都熟知,在柔性扫描中如果写了两个扫描坐标,如

B 15 S 7 0.1

B 16 S 7 0.1

是依次扫描两个坐标,无法做到同时,因此得到的是一张二维势能面,总扫描点数是两个坐标扫描点数的乘积,计算量较大。然而有时候我们只想同时扫描两个反应坐标,即两个坐标同时改变,得到一条曲线。

例如,找 [2+2] 环加成反应的过渡态经常会碰到这种问题,对于复杂的分子结构,手动调整过渡态初猜很难合适,此时使用 opt=ts 找到过渡态的成功率自然也不高,这时候我们可能就想,取柔性扫描势能曲线(面)上的突跃点作为过渡态初猜,然而去扫两个坐标得到一张二维势能面未免过于耗时。

对于这种问题,笔者以往采用了两种做法:(1)写了一个小程序来产生调整键长后新的结构(不仅仅是拉进 / 远两个原子),然而产生的结构我并不满意,算法还需改进或者仍有 bug。(2)若仅算一两步反应,那么就手动在 GaussView 里调整好两个键长,每次算完下载下来再调键长,这样扫描 5 个点就要下载、调整 5 次,甚是麻烦。

由于 G16 推出了广义内坐标(GIC)功能,于是笔者便研究了一下官网的说明

http://gaussian.com/gic/

琢磨出了一个同时扫描两个键长的文件模板。可能(大概率)这不是最优的方案(第二部分介绍另一种方案)。

为简洁起见,本文以甲醛和水的加成反应(真空中)为例,示范如何写输入文件。

示例

在这个反应中有两个主反应坐标(C−O 键和 O−H 键)同时在动,单独去扫描 C−O 键或者 O−H 键能量都会一直升高,并不会有突跃点。为防止混淆,短横线−符号左边始终表示甲醛分子中的原子,符号右边则表示水分子中的原子。

当然,化学直觉较强的同学可以直接构造出这个简单反应的合理过渡态初始结构,但这招对复杂结构很难见效,因此有必要掌握同时扫描两个反应坐标的技巧。假设我们希望 C−O 键和 O−H 键按如下距离逐渐减小

C−O: 1.8, 1.7, 1.6, 1.5

O−H: 1.7, 1.5, 1.3, 1.1

则利用 GIC 功能的输入文件如下:

输入文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p M062X/6-31G(d) nosymm geom=addgic

Title Card Required

0 1
 C    -1.13256115   -0.04504504    0.00000000
 H    -0.59939740   -0.97274996    0.00000000
 H    -2.20256115   -0.04504504    0.00000000
 O    -0.50552072    1.04600530    0.00000000
 O    -0.71081810   -0.28742611    1.83667749
 H    -0.42043150   -1.07235298    2.30695978
 H    -0.23246558    0.54490727    1.83667749

R15=R(1,5)
R47=R(4,7)
R15(value=1.8)
R47(value=1.7)

--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p opt M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read

R15(freeze)
R47(freeze)

--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read

R15(value=1.7)
R47(value=1.5)

--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p opt M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read

R15(freeze)
R47(freeze)

--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read

R15(value=1.6)
R47(value=1.3)

--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p opt M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read

R15(freeze)
R47(freeze)

--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read

R15(value=1.5)
R47(value=1.1)

--Link1--
%chk=methanal_water.chk
%mem=24GB
%nprocshared=24
#p opt M062X/6-31G(d) nosymm geom(allcheck,addGIC) guess=read

R15(freeze)
R47(freeze)

其中 R15 和 R47 是新起的变量名,可以改成其他名称。而 R(1,5) 则表示 1 号和 5 号原子之间的距离。注意 value 并不能限制键长严格等于给定的值,只能是十分接近,这是高斯官网提到的(不过这点对普通用户并不重要)。

第一步是先用 GIC 调整两根化学键至目标键长,调整完成后会做单点计算,下一个任务 –Link1– 读取之,将调整好的键长冻结、做限制性优化;完成后再用 GIC 调整两根键至新的键长,算单点,再 –Link1– 读取之,将调整好的键长冻结、做限制性优化。。。如此几步下来就能同时让两根键缩短了。从始至终只用了一个 gjf 文件,中途无需人为干预,可以认为是达到了目的。完成后用 GV 6.0 打开 log 文件

会发现有很多任务,这是用了 –Link1– 造成的,只能打开一个个 Optimization 查看结构和能量,看在哪个任务处发生了突跃。这点不是很方便,后续如果没有发现更优的方案,笔者会考虑写个脚本从 log 里提取信息,避免这一步操作。

回到这个问题,经过查看后发现第三个任务是能量最高点,下一步能量就会下降了。将第三个任务中最后一帧结构另存为、提取出来作为过渡态初猜,经过 10 步可以收敛到正确的过渡态上。

总结:本文用一个简单的反应展示了如何同时扫描两个反应坐标,该反应若仅扫描其中任何一个反应坐标都是得不到突跃点的。对于更密的扫描步长、及扫描键角等等,读者可根据文中提供的示例文件自己举一反三。

PS1:作为练习,感兴趣的读者可以自己试试找溴化氢 HBr 与乙烯或乙炔加成的过渡态,考验一下自己的化学直觉。若找不到的话,可试试本文的技巧。

PS2: GIC 功能仅在 >= G16 A 版本才有,G09 无法使用该功能。

第二篇摘抄

《一种用Gaussian 16中的GIC功能实现同时扫描多个坐标的方法》

简介

势能面扫描是我们用 Gaussian 常做的计算,一般可以分为刚性扫描和柔性扫描。

如果在柔性扫描中给定两个坐标,那么我们将会得到二维势能面。

但是有时候我们只希望两个坐标同时变化得到一条势能曲线,这可以通过使用 Gaussian 中的 GIC(广义内坐标)实现,上文给出了一个可行的解决方案。

但是之前方案的缺点是使用了 Link1,在用 GaussView 打开输出文件时不能很方便地显示能量的变化趋势,这在找能量极大,极小点时会带来困难。因此这里给出了一个新的方案,不使用 Link1,让势能曲线可以直观地显示出来。

这个新的方案同样需要使用 GIC,没有接触过 GIC 的同学可以在 Gaussian 官网上学习相关的资料。

GIC使 Gaussian 用户可以自定义一些结构参数,如键长、键角、二面角等。在基础的结构参数上还可以用数学运算,如加减乘除、平方开方、三角函数等,定义更复杂的结构参数。根据定义出的结构参数,我们可以做限制性优化,势能面扫描等等。

甲醛与水的加成反应

还是以甲醛与水的加成反应为例,介绍 GIC 输入文件的写法。

我们需要扫描的是 O4−H7 和 C1−O5 间的键长,因此需要对相关结构参数有一个定义:

RCO=R(1,5)

ROH=R(4,7)

RCO=R(1,5) 定义了 RCO 为 1 号原子和 5 号原子之间的距离。ROH=R(4,7) 定义了 ROH 为 4 号原子和 7 号原子之间的距离。注意新建变量名时,不建议用 R12、R34 这种变量名,因为可能会被高斯程序内部占用,引起冲突报错。与上文相同,我们希望 RCO 和 ROH 按照如下距离逐渐缩小:

RCO=1.8, 1.7, 1.6, 1.5, 1.4

ROH=1.7, 1.5, 1.3, 1.1, 0.9

正如前面提到的,我们不希望做二维势能面扫描,因此我们只能选择一个扫描坐标。这里选择 RCO,我们只需要对 RCO=R(1,5) 做一些改动即可:

RCO(NSteps=4,StepSize=-0.1)=R(1,5)

NSteps=4 指扫描步数为 4,StepSize=−0.1指扫描步长为−0.1。我们先将初始结构调整为 RCO=1.8,ROH=1.7,上面的命令可以使 RCO 按每一步 0.1 的步长减小,共减小 5 步,这正是我们所希望的。

在完成了 RCO 扫描的设置后,我们设置 ROH,让 ROH 能随着 RCO 变化而变化。这可以用 GIC 中的 Frozen 功能完成。使用 Frozen 功能的前提是找出我们需要固定的量。我们将 RCO 和 ROH 的值输入到 Excel 中,并作出趋势线:

得到 RCO 和 ROH 间满足的关系为

ROH=2*RCO-1.9

这个等式说明 2.0*RCO-ROH 在整个扫描过程中可以作为一个不变量。

我们定义一个量 F=2.0*RCO-ROH,并将其固定:

F(Frozen)=2.0*RCO-ROH

F(Frozen) 可以让 Gaussian 在优化过程中固定 F 这个值不变。在扫描过程中,每当我们有一个新的 RCO,由于 F 的限制,我们总会得到相应的 ROH。这就让 ROH 可以随 RCO 变化而变化。

使用 GIC 功能的完整的输入文件如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
%mem=24GB
%nprocshared=24
#p opt=(modredundant,GIC) M062X/6-31G(d) nosymm

Title Card Required

0 1
 C       -1.12689502   -0.06090992    0.04839780
 H       -0.59373127   -0.98861484    0.04839780
 H       -2.19689502   -0.06090992    0.04839780
 O       -0.49985459    1.03014042    0.04839780
 O       -0.71648423   -0.27156123    1.78827969
 H       -0.42609763   -1.05648810    2.25856198
 H       -0.25552070    0.59195719    1.67268115

RCO(NSteps=4,StepSize=-0.1)=R(1,5)
ROH=R(4,7)
F(Frozen)=2.0*RCO-ROH

由于没有使用 Link1,用 GaussView 打开输出文件不会显示多步任务。查看 Scan 的结果,可以很快地找到最高点,如下图所示:

将最高点作为初猜做过渡态搜索,即可找到甲醛和水加成的过渡态。另外,若把 opt=(modredundant,GIC) 替换成两个关键词 opt geom=addGIC,效果一样。

简单总结一下,写同时扫描多个坐标所需 Gaussian 输入文件的通用步骤为:

1. 首先指定第一个扫描坐标, 例如

RCO(NSteps=4,StepSize=-0.1)=R(1,5)

2. 定义下一个坐标,并根据新坐标和旧坐标的关系定义出一个不变量,

ROH=R(4,7)

F(Frozen)=2.0*RCO-ROH

这一步可以通过函数的拟合来完成,相当于找到一个函数 F(RCO,ROH)=0。

3. 重复 2,直到所有坐标都定义完毕。

三个水分子间的质子转移反应

下面以三个水分子间的质子转移反应为例,演示如何同时扫描多个坐标。

我们首先对三个水分子的团簇做结构优化,得到稳定结构。

得到 O−H 间键长为 0.97819,不成键的 O−H 原子间距离为 1.83887。

为了得到 3 个质子同时转移的过渡态,需要同时缩短 H2−O4,H6−O7,H9−O1 间距离。扫描步数预设为 10 步,可得到步长约为−0.08。将 H2−O4 间距离定为扫描坐标,为了让 H6−O7,H9−O1 间距离随着 H2−O4 间距离减小而减小,需要引入两个不变量 F1、F2:

F1(Frozen)=RH2O4-RH6O7

F2(Frozen)=RH2O4-RH9O1

根据前文的介绍,很容易写出 GIC 计算的输入文件:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
%nprocshared=24
%mem=24GB
#p opt=(modredundant,GIC) m062x/6-31g(d,p) nosymm

Title Card Required

0 1
 O      -2.12109900    0.15275700    0.36236700
 H      -1.26182000   -0.15720700    0.01249000
 H      -2.74456300    0.04217400   -0.36151100
 O       0.46851300    0.12482800   -0.54235200
 H       1.12745300   -0.25541100    0.04672800
 H       0.37017800    1.05352900   -0.25131900
 O      -0.56877500    2.41551400    0.54801000
 H      -0.86226500    3.12772800   -0.02774600
 H      -1.33253600    1.81105900    0.63231100

RH2O4(NSteps=10,StepSize=-0.08)=R(2,4)
RH6O7=R(6,7)
F1(Frozen)=RH2O4-RH6O7
RH9O1=R(9,1)
F2(Frozen)=RH2O4-RH9O1

用 GaussiView 打开输出文件,可以看到结构随扫描过程的变化及势能曲线。

取出能量最高点的结构做过渡态搜索,优化 3 步之后就能收敛到相应过渡态。可见柔性扫描可以给出相当好的过渡态初猜。

综上,GIC 是一个 Gaussian 中很有用的功能,GIC 结合柔性扫描对过渡态搜索有很大的帮助。

第三篇摘抄

http://bbs.keinsci.com/thread-13552-1-1.html

这个帖子里面写了同时缩短两个键长时使用的输入文件。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# opt=(modredundant,gdiis), ub3lyp/genecp scf=xqc

Title Card Required

(略)

B 60 56 S 40 -0.05000
B 61 54 S 40 -0.05000

-N O C H
6-31G(d)
****
Fe 0
LanL2DZ
****

Fe 0
LanL2DZ
_
_

我看sobereva老师没提出啥意见,看样子可以这么写。