请选择 进入手机版 | 继续访问电脑版

计算化学公社

 找回密码
 现在注册!
查看: 5159|回复: 15

[Multiwfn资源与经验] 利用Multiwfn计算倾斜、扭曲环的NICS_ZZ

[复制链接]

2万

帖子

25

威望

2万

eV
积分
46931

管理员

公社社长

发表于 2014-11-20 07:26:06 | 显示全部楼层 |阅读模式
利用Multiwfn计算倾斜、扭曲环的NICS_ZZ

文/Sobereva @北京科音
First release: 2014-Nov-19  Last update: 2019-Aug-30

NICS、NICS(1)_ZZ的定义在《衡量芳香性的方法以及在Multiwfn中的计算》(http://sobereva.com/176)的3.1节有详细介绍。NICS(1)_ZZ具体指的是垂直于环平面上方1埃处在垂直于环平面方向上的磁屏蔽值的负值。对于平面体系,这很好算,比如体系在XY平面上,那么先计算出环的中心,然后把坐标的Z值加1,放置个Bq原子(对于Gaussian来说),然后算NMR,读取Bq原子的磁屏蔽张量的ZZ分量取负号即可。但是,被研究的环如果不在XY/YZ/XZ平面上,或者环本身就是扭曲的而非完美平面,那么计算NICS(1)_ZZ就很麻烦了,光是确定要计算的Bq原子的位置就不好搞。比如下面的例子,C36的一个异构体,假设要研究红点标注的5,6,25,28,27,29这个六元环,它既是斜着的,又不是纯平面。

1.png

虽然也可以旋转这个团簇,让这个环正好基本平行于XY平面上,然后按照常规步骤算NICS(1)_ZZ,但是如果要研究的环很多,每研究一个环都这么转一次结构,实在太麻烦。而利用Multiwfn (http://sobereva.com/multiwfn),则可以很方便地解决计算倾斜非平面环的NICS(1)_ZZ的问题。下面就以上面那个环作为例子,来演示一遍操作。


我们先计算出环中心。有很多方法计算环中心,两种比较常用,其一是取环的质心,其二是取环的AIM理论中的环临界点(RCP)位置。由于前者不需要波函数文件,有结构就行了,为了省事我们此例就用质心。后者也可以十分方便快速地通过Multiwfn的拓扑分析功能得到,可参见《使用Multiwfn做拓扑分析以及计算孤对电子角度》(http://sobereva.com/108 )。

启动Multiwfn,输入以下命令
Fullerene_No.1-C2.pdb  // C36结构文件
100
21  // 计算各种基于几何结构的属性
5,6,25,27-29  // 组成环的原子的序号,这些原子将会被用来计算各种基于几何结构的属性
从输出文件中我们看到了此环的质心(对此例等价于几何中心)
Center of mass (X/Y/Z): 0.96300000 -1.89100000 0.43750000 Angstrom

我们把质心坐标复制下来,然后输入
q  //返回
24  //辅助计算非平面体系的NICS_ZZ
0.96300000 -1.89100000 0.43750000 //粘贴刚才的质心坐标
因为这个环本身不是完全平面的,因此也没有办法唯一地定义一个平面来代表它。但可以通过环上的原子根据最小二乘法拟合出一个最具有代表性的平面。因此我们输入5,6,25,28-29,使Multiwfn根据环上的这些原子的坐标进行拟合。

屏幕上显示
RMS error of the plane fitting:    0.296080 Angstrom
The unit normal vector is   -0.29531754    0.91625276   -0.27068142

The X,Y,Z coordinate of the points below and above 1 Angstrom of the plane from
the point you defined, respectively:
    1.2583175359   -2.8072527609    0.7081814199 Angstrom
    0.6676824641   -0.9747472391    0.1668185801 Angstrom

告诉了平面的拟合误差(如果平面本身就是纯平的则误差为0)、此平面的单位法向量,以及从刚才输入的环中心沿着法向量向上和向下分别移动1埃后的坐标。这两个坐标就可以设为Bq原子来计算NICS(1)_ZZ了。假设整个体系大体是个平面,那么用这两个点的NICS(1)_ZZ值取平均是比较合理的。不过当前体系是笼状体系,所以我们只应当取处在环外的那个点来算NICS(1)_ZZ。

建议先别关Multiwfn,因为待会儿还得把算好的Bq点的屏蔽张量输进去。如果关了的话,之后还得重新做上述操作。

我们基于C36的结构创建一个Gaussian输入文件,然后把刚才得到的1.2583175359   -2.8072527609    0.7081814199这个点作为Bq原子的位置,此时输入文件如下(此文仅作为示例,纯粹为节省计算量用STO-3G,实际研究中决不能低于6-31G*,否则结果没法用)
# B3LYP/STO-3G NMR nosymm
[空行]
Generated by Multiwfn
[空行]
0 1
C     -3.06300000   -0.54400000    0.11800000
C     -2.20500000   -1.61200000   -0.25800000
...[略]
C      2.55900000   -0.14400000   -1.38400000
C      2.29700000    1.21800000   -1.08100000
Bq 1.2583175359   -2.8072527609    0.7081814199


用gview查看一下,可见Bq原子位置很合适,确实是在环平面上方1埃处
2.png

注意输入文件里的nosymm很重要,如果不写的话,Gaussian可能会自动旋转体系到标准朝向下,导致垂直于环平面的方向在计算前后不同,Multiwfn也就没法正确提取垂直于环平面方向的分量了。

用Gaussian计算此文件,看到NMR部分Bq原子的输出
     37  Bq   Isotropic =    -6.1166   Anisotropy =    25.9219
   XX=     0.1309   YX=    -0.5201   ZX=    -1.2347
   XY=     3.2617   YY=   -23.2564   ZY=     7.2520
   XZ=     6.9045   YZ=   -35.8046   ZZ=     4.7756
   Eigenvalues:   -29.4361    -0.0784    11.1646


切换回Multiwfn窗口,将屏蔽张量依次输入进去:
0.1309,-0.5201,-1.2347
3.2617,-23.2564,7.2520
6.9045,-35.8046,4.7756

(秘籍:想图省事的话,把Gaussian输出文件的屏蔽张量那三行一次性直接粘贴到Multiwfn窗口里也可以)

最后会显示
The shielding value normal to the plane is      -12.3700817601
The NICS_ZZ value is thus       12.3700817601

即NICS(1)_ZZ值为12.37。

其实本例在进入主功能100的子功能24之前可以不去特意用子功能21。因为当前用的环中心就是几何中心,而且定义拟合平面用的原子就是环上的所有原子。对于这种情况,直接进入子功能24,让你输入环中心的时候直接敲回车,之后输入环上的原子序号后,程序自动就会将这些原子的几何中心作为环中心。

对于平面上方和下方1埃处的NICS(1)_ZZ都需要计算的情况,你在输入屏蔽张量的界面里,把平面上方1埃处的屏蔽张量粘贴进去,得到的就是上方1埃处的NICS_ZZ(1),如果把下方1埃处的屏蔽张量粘贴进去,得到的就是下方1埃处的NICS_ZZ(1)。

顺带一提,虽然NICS(0)_ZZ我很不推荐,但如果非要算的话,步骤和上述一样,只不过把Bq放在环中心即可。要计算比如NICS(5)_ZZ也可以,也就是把环中心坐标加上5*单位法向量,就得到了要算的垂直于平面上方5埃处的点的坐标。之后的步骤和文中一样,输入那个点的屏蔽张量即可。

如果要研究的环有多个,那么按照上述步骤,确定出各个环要计算的Bq原子位置,都写进Gaussian输入文件里,只需要Gaussian算一次即可,而不需要每研究一个环都单独算一次。

评分

参与人数 3eV +8 收起 理由
zsu007 + 5 赞!
aqhuangry + 2
captain + 1 赞!

查看全部评分

北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。约6000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

185

帖子

0

威望

1990

eV
积分
2175

Level 5 (御坂)

发表于 2015-4-22 09:25:33 | 显示全部楼层
我想问下,平时说的x, y, z方向跟xx,yy,zz方向 有什么可以联系的地方?

2万

帖子

25

威望

2万

eV
积分
46931

管理员

公社社长

 楼主| 发表于 2015-4-22 14:22:59 | 显示全部楼层
赵云跳槽 发表于 2015-4-22 09:25
我想问下,平时说的x, y, z方向跟xx,yy,zz方向 有什么可以联系的地方?

X,Y,Z是方向,而XX,YY,ZZ都是二阶张量的分量。
下图中,B'是感生磁场,B0是外加磁场,σ是磁屏蔽张量。
1.png

根据矩阵运算关系可见,比如σXY就体现了在Y方向外加B0y磁场时,对X方向的感生磁场会有多少贡献
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。约6000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

10

帖子

0

威望

87

eV
积分
97

Level 2 能力者

发表于 2016-7-9 20:17:15 | 显示全部楼层
老师,您好,我想问下除了用Multiwfn还没有没别的确定要计算的Bq原子的位置的方法,因为我下了个Multiwfn不太会用,谢谢老师。

2万

帖子

25

威望

2万

eV
积分
46931

管理员

公社社长

 楼主| 发表于 2016-7-10 07:09:57 | 显示全部楼层
haipwei 发表于 2016-7-9 20:17
老师,您好,我想问下除了用Multiwfn还没有没别的确定要计算的Bq原子的位置的方法,因为我下了个Multiwfn不 ...


没有
不会就学啊!谁天生就会。你用别的软件天生就会?其它哪个程序像multiwfn一样手册里又有一大堆例子,还特意写博文去详细示例介绍?
Multiwfn这么简单好用的程序都不会用,其它连个像样英文手册都没有的程序能会用才怪。
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。约6000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

37

帖子

0

威望

254

eV
积分
291

Level 3 能力者

发表于 2017-10-20 08:49:40 | 显示全部楼层
请问老师几个问题, 2 (2).png 2 (1).png 一、(图一)本征值与屏蔽张量的关系以及与XX、YY、ZZ的关系,二、(图二)G09计算出Bq的Iso值,是不是代表NICS(1)的值,  如果用G09输出的NMR结果如何得到NICS(1)_ZZ的值以及与Multiwfn软件结果的区别

2万

帖子

25

威望

2万

eV
积分
46931

管理员

公社社长

 楼主| 发表于 2017-10-20 08:52:25 | 显示全部楼层
weizheng007 发表于 2017-10-20 08:49
请问老师几个问题,一、(图一)本征值与屏蔽张量的关系以及与XX、YY、ZZ的关系,二、(图二) ...

1 看这里的描述
将核独立化学位移(NICS)分解为sigma和pi轨道的贡献
http://sobereva.com/145
2 如果Bq那个位置是在平面上方1埃处,是NICS(1)的值
其余问题看
衡量芳香性的方法以及在Multiwfn中的计算
http://sobereva.com/176
里关于NICS的描述
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。约6000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

37

帖子

0

威望

254

eV
积分
291

Level 3 能力者

发表于 2017-10-20 08:58:47 | 显示全部楼层
sobereva 发表于 2017-10-20 08:52
1 看这里的描述
将核独立化学位移(NICS)分解为sigma和pi轨道的贡献
http://sobereva.com/145

看了才迷糊的, 1、博文有的XX、YY、ZZ以外的数值为零,   有的不为零,但Iso均为对角元的平均值,Aniso也有讲述,但没有看到本征值讲法,2、Gaussian输出结果和Multiwfn结果不一样更迷糊

2万

帖子

25

威望

2万

eV
积分
46931

管理员

公社社长

 楼主| 发表于 2017-10-20 17:52:11 | 显示全部楼层
weizheng007 发表于 2017-10-20 08:58
看了才迷糊的, 1、博文有的XX、YY、ZZ以外的数值为零,   有的不为零,但Iso均为对角元的平均值,Aniso ...

矩阵本征值是线性代数里的基本内容,对磁屏蔽张量对角化一下就完了
正因为倾斜的环你直接从高斯输出文件里读不到垂直于平面的分量,所以才需要用Multiwfn基于高斯的输出数据进行变换来得到
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。约6000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

3

帖子

0

威望

128

eV
积分
131

Level 2 能力者

发表于 2019-3-12 14:24:34 | 显示全部楼层
本帖最后由 克里斯保 于 2019-3-12 14:31 编辑

sob老师,您好,我想问一下,为什么
37 Bq Isotropic = -5.9828 Anisotropy = 26.8438
XX= 0.4617 YX= -0.6483 ZX= -1.1534
XY= 2.4525 YY= -23.5551 ZY= 6.8922
XZ= 6.4117 YZ= -36.8502 ZZ= 5.1448
Eigenvalues: -30.0625 0.2009 11.9131

这个计算出来为什么不能将zz那个值取负呢?而是再导入Multiwfn中导出-13.38这个值?是因为初始输入的时候并没有将环调到平面上的原因吧

2万

帖子

25

威望

2万

eV
积分
46931

管理员

公社社长

 楼主| 发表于 2019-3-12 14:49:46 | 显示全部楼层
克里斯保 发表于 2019-3-12 14:24
sob老师,您好,我想问一下,为什么
37 Bq Isotropic = -5.9828 Anisotropy = 26.8438
XX= 0.4617 YX= - ...

纯平面体系,而且平行于XY的时候,直接读ZZ取负值和用Multiwfn给出来的完全一样
本帖重点是非纯平面倾斜环的NICS_ZZ的计算
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。约6000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

5

帖子

0

威望

121

eV
积分
126

Level 2 能力者

发表于 2019-8-28 09:24:03 | 显示全部楼层
sobereva老师您好,我想知道NICS_ZZ计算的数学原理是什么?我计算一个苯环六个原子不在一个平面上,用间隔的三个原子去定义一个平面也感觉不太科学,我是通过六个原子的坐标最小二乘法来拟合这个平面的。能否将这个平面方程或者它的法向量代入到Multiwfn来计算?或者有没有数学公式,我用matlab之类的软件来计算?

2万

帖子

25

威望

2万

eV
积分
46931

管理员

公社社长

 楼主| 发表于 2019-8-28 23:25:18 | 显示全部楼层
trevormx 发表于 2019-8-28 09:24
sobereva老师您好,我想知道NICS_ZZ计算的数学原理是什么?我计算一个苯环六个原子不在一个平面上,用间隔 ...

公式在手册3.100.24节提到了,你拟合了平面方程就有了法矢量,按照公式对磁屏蔽张量做个简单的矩阵运算就能得到
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。约6000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

2万

帖子

25

威望

2万

eV
积分
46931

管理员

公社社长

 楼主| 发表于 2019-8-30 04:29:10 | 显示全部楼层
trevormx 发表于 2019-8-28 09:24
sobereva老师您好,我想知道NICS_ZZ计算的数学原理是什么?我计算一个苯环六个原子不在一个平面上,用间隔 ...

对Multiwfn程序刚做了更新,现已可以在官网下载。之前版本定义环平面是通过三个原子定义的,但这样不准确,有时候很难找合适的三个原子。目前Multiwfn官网上最新的Multiwfn改进了这点,用户可以输入一批原子(比如环上的所有原子),程序会用这些点通过最小二乘法拟合出最有代表性的平面,明显合理多了。

相应地,本文的内容和例子都进行了更新。
北京科音自然科学研究中心:http://www.keinsci.com  致力于计算化学的发展和传播,不定期开办各层次量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训。欢迎加入“北京科音”微信公众号获取培训最新消息和计算化学资讯
思想家公社QQ群,1号:18616395,2号:466017436。约6000人,专门交流理论、计算化学。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一概不批。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最流行的量子化学波函数分析程序)

此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!

5

帖子

0

威望

121

eV
积分
126

Level 2 能力者

发表于 2019-9-1 19:37:07 | 显示全部楼层
sobereva 发表于 2019-8-30 04:29
对Multiwfn程序刚做了更新,现已可以在官网下载。之前版本定义环平面是通过三个原子定义的,但这样不准确 ...

辛苦SOBEREVA老师了!谢谢啦
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949-1号 )

GMT+8, 2019-10-21 09:10 , Processed in 0.174773 second(s), 28 queries .

快速回复 返回顶部 返回列表