专栏首页量子化学用Gaissian16中的GIC功能实现翻转过程的势能面扫描

用Gaissian16中的GIC功能实现翻转过程的势能面扫描

一、简介

势能面扫描在计算化学领域中有很多应用场景,本公众号之前也推送过相关介绍:

势能面扫描前需要用户对扫描坐标有一个明确的定义。在Gaussian 16的广义内坐标(GIC)功能出现之前,我们只能对一些简单的结构参数,如笛卡尔坐标、键长、键角、二面角做势能面扫描。GIC的出现让我们可以定义更加复杂的扫描坐标。

翻转过程在自由基和三取代的N、P体系中很常见。以氨气分子的翻转为例,三个H原子构成一个平面,翻转过程可以看成是N原子到此平面的距离不断变化的过程:

因此一个很自然的想法是:将N原子到三个H原子形成的平面的距离设置为扫描坐标,以此坐标做势能面扫描。实际操作中我并没有将距离作为扫描坐标,详细的扫描坐标的设置会在下文中解释。不使用GIC的话,这是很难实现的,使用GIC后问题就能迎刃而解。

二、相关数学公式的推导

为了得到扫描坐标在数学上的定义,我们需要对相关公式进行推导,先求出不共线三点的平面方程,再求出第四个点到这个平面的距离。

2.1 求过(不共线)三点的平面方程

已知三个点的坐标

得到两个向量

则垂直于此平面的法向量为

其中

则点法式平面方程可写为

将平面方程化为一般式

可求得

综上,不共线三点的平面方程为

其中

2.2 点到平面的距离

设空间中第四个点的坐标为

则此点到平面

的距离为

三、氨气分子翻转过程的势能面扫描

以氨气分子为例,介绍GIC在扫描翻转过程的势能面上的应用。

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

%nprocshared=24
%mem=24GB
#p opt=(modredundant,GIC) 6-31g(d,p) nosymm m062x

Title Card Required

0 1
 H     1.80953700   -0.74056700    0.00001400
 H     1.80955100    0.66571400    0.81195400
 H     1.80951800    0.66574100   -0.81197200
 N     1.41988000    0.19697800    0.00000500

X1=X(1)
Y1=Y(1)
Z1=Z(1)
X2=X(2)
Y2=Y(2)
Z2=Z(2)
X3=X(3)
Y3=Y(3)
Z3=Z(3)
X4=X(4)
Y4=Y(4)
Z4=Z(4)
NA=(Y2-Y1)*(Z3-Z1)-(Y3-Y1)*(Z2-Z1)
NB=(Z2-Z1)*(X3-X1)-(Z3-Z1)*(X2-X1)
NC=(X2-X1)*(Y3-Y1)-(X3-X1)*(Y2-Y1)
ND=-1.0*(NA*X1+NB*Y1+NC*Z1)
Dist(NSteps=20,StepSize=-0.07363)=(NA*X4+NB*Y4+NC*Z4+ND)/SQRT(NA**2+NB**2+NC**2)

其中

X1=X(1)

...

Z4=Z(4)

的作用是用自定义的变量X1,Y1,Z1等表示每个原子的笛卡尔坐标。

NA=(Y2-Y1)*(Z3-Z1)-(Y3-Y1)*(Z2-Z1)

NB=(Z2-Z1)*(X3-X1)-(Z3-Z1)*(X2-X1)

NC=(X2-X1)*(Y3-Y1)-(X3-X1)*(Y2-Y1)

ND=-1.0*(NA*X1+NB*Y1+NC*Z1)

的作用是计算三个H原子构成的平面方程。

Dist(NSteps=20,StepSize=-.07363)=(NA*X4+NB*Y4+NC*Z4+ND)/SQRT(NA**2+NB**2+NC**2)

的作用有两个,一是定义扫描坐标Dist,二是给定扫描的步数和步长。这是整个流程中最关键的部分。我们在2.2中计算的是点到平面的距离,是一个非负数。而这里Dist的表达式与2.2中有一个区别是去掉了绝对值符号,这使得Dist的值可以是负数。

易得按照新的Dist的定义,N原子在3个H原子构成的平面的两侧时,Dist有相反的符号。反之亦然,Dist改变大小和符号时,N原子就可以在3个H原子构成平面的两侧移动。

那么如何得到初始结构下Dist的大小和符号,让我们对扫描的步数和步长有一个定义呢?根据右手定则可以确定Dist的符号,不借助其他方法Dist的大小只能目测。一个比较简单的方法是用半经验方法算一个单点,Gaussian就会输出各GIC的值。

用半经验方法算单点的输入文件如下:

%nprocshared=24
%mem=24GB
#p nosymm pm6 geom=addGIC

Title Card Required

0 1
 H    1.80953700   -0.74056700    0.00001400
 H    1.80955100    0.66571400    0.81195400
 H    1.80951800    0.66574100   -0.81197200
 N    1.41988000    0.19697800    0.00000500

X1=X(1)
Y1=Y(1)
Z1=Z(1)
X2=X(2)
Y2=Y(2)
Z2=Z(2)
X3=X(3)
Y3=Y(3)
Z3=Z(3)
X4=X(4)
Y4=Y(4)
Z4=Z(4)
NA=(Y2-Y1)*(Z3-Z1)-(Y3-Y1)*(Z2-Z1)
NB=(Z2-Z1)*(X3-X1)-(Z3-Z1)*(X2-X1)
NC=(X2-X1)*(Y3-Y1)-(X3-X1)*(Y2-Y1)
ND=-1.0*(NA*X1+NB*Y1+NC*Z1)
Dist=(NA*X4+NB*Y4+NC*Z4+ND)/SQRT(NA**2+NB**2+NC**2)

注意,由于只需要做单点计算,RouteSection部分有改动。并且在Dist的定义部分去掉了势能面扫描的内容。

单点计算结束后,我们可以在输出文件中找到以下内容:

NA     GIC-1         -8.1553
NB     GIC-2          0.0
NC     GIC-3          0.0002
ND     GIC-4         27.8873
Dist   GIC-5          0.7363

可见Gaussian默认输出所有自定义的GIC。借助这个特性,我们也可以先随意给定NSteps和StepSize的值,提交Gaussian作业后马上kill掉进程,这样也可以在输出文件中找到初始结构下GIC的值。

在得到Dist数值之后,我们就可以设置扫描步数和步长为:

NSteps=20,StepSize=-0.07363

计算完毕后,用GaussView打开输出文件,可以看到氨气分子结构随扫描过程的变化及势能曲线。

取出势能曲线最高点的结构做过渡态搜索,优化两步就能收敛到翻转过渡态:

四、三取代膦翻转过程的势能面扫描

与三取代的胺的性质相似,三取代膦也可以有翻转过程。有文献指出,膦被氧化生成自由基正离子之后,翻转过程的能垒会降低。这里我们首先完成膦的翻转过程的势能面扫描,再做翻转过程的过渡态搜索。

中性膦的势能面扫描的输入文件如下:

%nprocshared=24
%mem=48GB
#p opt=(modredundant,GIC) m062x/6-31g(d,p) scrf=(smd,solvent=acetonitrile) nosymm

title

0 1
 P     -0.08670900    1.26287500   -0.95222100
 C     -1.47074400    0.32462200   -0.16515700
 C     -2.64121000    0.11500200   -0.89937700
 C     -1.39311100   -0.16786900    1.14315500
 C     -3.71993900   -0.56673800   -0.33663100
 H     -2.70680000    0.48251800   -1.92025200
 C     -2.46631700   -0.85288100    1.70502000
 H     -0.48207200   -0.02125600    1.71892100
 C     -3.63247400   -1.05223900    0.96516400
 H     -4.62347800   -0.72323100   -0.91717100
 H     -2.39458900   -1.23297200    2.71923500
 H     -4.46832100   -1.58825600    1.40318700
 C      1.36920700    0.36444500   -0.26127900
 C      1.61072200   -0.96464400   -0.67382300
 C      2.27789400    0.98437200    0.60369100
 C      2.74461700   -1.62383700   -0.19552200
 C      3.40744700    0.31243600    1.06745000
 H      2.11239800    2.00712700    0.92566900
 C      3.64063200   -0.99858600    0.66843900
 H      2.92760800   -2.64729400   -0.51196700
 H      4.09752100    0.81484900    1.73762200
 H      4.51618400   -1.53337000    1.02229700
 C     -0.13426100    2.78566600    0.10023900
 H     -1.07659600    3.29649900   -0.11053200
 H     -0.08645200    2.55787100    1.16860700
 H      0.68341700    3.45840600   -0.16864500
 C      0.66647000   -1.68997800   -1.59928300
 H     -0.25773600   -1.96742600   -1.08128400
 H      0.37874200   -1.06348500   -2.45041700
 H      1.12722100   -2.60172900   -1.98447600

X1=X(2)
Y1=Y(2)
Z1=Z(2)
X2=X(23)
Y2=Y(23)
Z2=Z(23)
X3=X(13)
Y3=Y(13)
Z3=Z(13)
X4=X(1)
Y4=Y(1)
Z4=Z(1)
NA=(Y2-Y1)*(Z3-Z1)-(Y3-Y1)*(Z2-Z1)
NB=(Z2-Z1)*(X3-X1)-(Z3-Z1)*(X2-X1)
NC=(X2-X1)*(Y3-Y1)-(X3-X1)*(Y2-Y1)
ND=-1.0*(NA*X1+NB*Y1+NC*Z1)
Dist(NSteps=20,StepSize=-0.16)=(NA*X4+NB*Y4+NC*Z4+ND)/SQRT(NA**2+NB**2+NC**2)

计算完毕后,用GaussView打开势能面扫描的输出文件,可以看到结构随扫描过程的变化及势能曲线。

取出势能曲线最高点的结构做过渡态搜索,优化三步就能收敛到翻转过渡态:

自由基正离子膦的势能面扫描的输入文件如下:

%nprocshared=24
%mem=48GB
#p opt=(modredundant,GIC) um062x/6-31g(d,p) scrf=(smd,solvent=acetonitrile) nosymm

title

1 2
 P      0.11694400    1.13577500    0.56859500
 C      1.56222700    0.22089900    0.07077600
 C      2.82938400    0.69566800    0.44800900
 C      1.43525800   -0.95614100   -0.68416900
 C      3.96295700   -0.00178100    0.05363800
 H      2.92984100    1.59576800    1.04687900
 C      2.57899900   -1.64232500   -1.07070700
 H      0.45518900   -1.31362900   -0.98503900
 C      3.83790800   -1.16810400   -0.70181100
 H      4.94362000    0.36112400    0.33926900
 H      2.48792500   -2.54443000   -1.66494700
 H      4.72693800   -1.71044700   -1.00496300
 C     -1.42137300    0.35948300    0.11779600
 C     -1.79142200   -0.86909900    0.71024200
 C     -2.26044300    1.00942000   -0.80274800
 C     -3.00255300   -1.43274500    0.31044200
 C     -3.45971400    0.42077300   -1.17677800
 H     -1.97302900    1.96156900   -1.23585100
 C     -3.82502300   -0.80361900   -0.62103400
 H     -3.31110200   -2.37555300    0.75066700
 H     -4.10351500    0.91537000   -1.89475800
 H     -4.76363900   -1.26757500   -0.90466800
 C      0.21246100    2.87741700    0.09391000
 H      1.10895300    3.31314400    0.53803100
 H      0.27677500    2.95886300   -0.99456000
 H     -0.66737000    3.40215200    0.46823700
 C     -0.94138000   -1.55987300    1.74591200
 H     -0.08734000   -2.07014000    1.28988800
 H     -0.55120600   -0.85359000    2.48523900
 H     -1.52990900   -2.30908900    2.27677400

X1=X(2)
Y1=Y(2)
Z1=Z(2)
X2=X(23)
Y2=Y(23)
Z2=Z(23)
X3=X(13)
Y3=Y(13)
Z3=Z(13)
X4=X(1)
Y4=Y(1)
Z4=Z(1)
NA=(Y2-Y1)*(Z3-Z1)-(Y3-Y1)*(Z2-Z1)
NB=(Z2-Z1)*(X3-X1)-(Z3-Z1)*(X2-X1)
NC=(X2-X1)*(Y3-Y1)-(X3-X1)*(Y2-Y1)
ND=-1.0*(NA*X1+NB*Y1+NC*Z1)
Dist(NSteps=20,StepSize=-0.08964)=(NA*X4+NB*Y4+NC*Z4+ND)/SQRT(NA**2+NB**2+NC**2)

计算完毕后,用GaussView打开势能面扫描的输出文件,可以看到结构随扫描过程的变化及势能曲线:

取出势能曲线最高点的结构做过渡态搜索,优化三步就能收敛到翻转过渡态:

易知膦被氧化为自由基正离子后翻转过程变得更加容易了:

五、总结

使用GIC可以让我们快速实现复杂反应坐标的势能面扫描。面对翻转过程这一类扫描坐标很复杂的情况,如果不使用GIC,势能面扫描是相当难实现的。势能面扫描给出的初猜也可以让我们很快地找到翻转过程的过渡态。

参考文献

Kyle D. R., Daniel H. E. , and Alexander T. R.* J. Am. Chem. Soc. 2013, 135, 9354−9357

本文分享自微信公众号 - 量子化学(quantumchemistry),作者:zxli

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-10-29

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 【双清论坛】国内计算化学方法与软件开发(Ⅲ)

    本系列的三篇推文已经推送完成,小编其实只是担任了一个搬运工。通过这三篇文章,我们能够对国内的理论化学方法的研究有个大概的了解。国内很多课题组的方法和软件在自己的...

    用户7592569
  • 二氢卡宾单重态和三重态稳定性的比较

    当然,教材上的图是重绘的,原图来自J. Am. Chem. Soc. 1996, 118, 9908-9914。本文就试着重复该图,并将关于自旋极化单重态的计算...

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

    势能面扫描可以用来研究势能随少数几个坐标变化的情况,分为刚性扫描和柔性扫描两种。刚性扫描是指被扫描的坐标变化后,不改变未被扫描的坐标,相当于做一系列单点能计算。...

    用户7592569
  • 快速学习Redis-Jedis客户端

    版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明...

    cwl_java
  • 【音视频连载-009】第二季 FFmpeg 打造简易播放器

    随着以抖音为代表的短视频应用逐渐深入到我们日常生活中来,直播带货、网络会议、在线教育等音视频场景也变得更加常见,可以说音视频技术已经和我们的生活息息相关了。

    glumes
  • Redis开发规范

    以业务名(或数据库名)为前缀(防止key冲突),用冒号分隔,比如业务名:表名:id

    业余草
  • Java连接Redis

    乐心湖
  • Java Redis系列3(Jedis的使用+jedis连接池技术)

    一只胡说八道的猴子
  • Tuxedo域连接错误引出的netstat不同平台参数的区别

    最近测试同事有一个需求,搭建一套Tuxedo域,连接网关,通过其和其他域交互,搭建过程中配置正确,但域连接报错,模拟过程如下,假设本地域名称为LOCALDOM,...

    bisal
  • java操作redis: 将string、list、map、自定义的对象保存到redis中

    似水的流年

扫码关注云+社区

领取腾讯云代金券