计算化学公社
标题: 利用Multiwfn计算倾斜、扭曲环的NICS_ZZ [打印本页]
作者Author: sobereva 时间: 2014-11-20 07:26
标题: 利用Multiwfn计算倾斜、扭曲环的NICS_ZZ
利用Multiwfn计算倾斜、扭曲环的NICS_ZZ
Using Multiwfn to calculate NICS_ZZ of tilted andtwisted rings
First release: 2014-Nov-19 Last update: 2023-Aug-11
1 前言
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这个六元环，它既是斜着的，又不是纯平面。
(, 下载次数 Times of downloads: 196)
虽然也可以旋转这个团簇，让这个环正好基本平行于XY平面上，然后按照常规步骤算NICS(1)_ZZ，但是如果要研究的环很多，每研究一个环都这么转一次结构，实在太麻烦。而利用Multiwfn (
http://sobereva.com/multiwfn)，则可以很方便地解决计算倾斜非平面环的NICS(1)_ZZ的问题。下面就以上面那个环作为例子，来演示一遍操作。
2 实例
我们先计算出环中心。有很多方法计算环中心，两种比较常用，其一是取环的质心，其二是取环的AIM理论中的环临界点(RCP)位置。由于前者不需要波函数文件，有结构就行了，为了省事我们此例就用质心。后者也可以十分方便快速地通过Multiwfn的拓扑分析功能得到，可参见《使用Multiwfn做拓扑分析以及计算孤对电子角度》（
http://sobereva.com/108 ）。
启动Multiwfn，输入以下命令
Fullerene_No.1-C2.pdb // C36结构文件
100 //其它功能（Part 1）
21 // 计算各种基于几何结构的属性
5,6,25,27-29 // 组成环的原子的序号，这些原子将会被用来计算各种基于几何结构的属性
从输出文件中我们看到了此环的质心（对此例等价于几何中心）
Center of mass (X/Y/Z): 0.96300000 -1.89100000 0.43750000 Angstrom
我们把质心坐标复制下来，退回到Multiwfn主菜单，然后输入
25 //离域性与芳香性分析
4 //辅助计算非平面体系的NICS_ZZ
0.96300000 -1.89100000 0.43750000 //粘贴刚才的质心坐标
因为这个环本身不是完全平面的，因此也没有办法唯一地定义一个平面来代表它。但可以通过环上的原子根据最小二乘法拟合出一个最具有代表性的平面。因此我们输入5,6,25,27-29，使Multiwfn根据环上的这些原子的坐标进行拟合。
屏幕上显示
RMS error of the plane fitting: 0.194520 Angstrom
The unit normal vector is 0.31391747 -0.90899632 0.27419245
The X,Y,Z coordinate of the points below and above 1 Angstrom of the plane from
the point you defined, respectively:
0.6490825274 -0.9820036752 0.1633075459 Angstrom
1.2769174726 -2.7999963248 0.7116924541 Angstrom
告诉了平面的拟合误差（如果平面本身就是纯平的则误差为0）、此平面的单位法向量，以及从刚才输入的环中心沿着法向量向上和向下分别移动1埃后的坐标。这两个坐标就可以设为Bq原子来计算NICS(1)_ZZ了。假设整个体系大体是个平面，那么用这两个点的NICS(1)_ZZ值取平均是比较合理的。不过当前体系是笼状体系，所以我们只应当取处在环外的那个点来算NICS(1)_ZZ。
建议先别关Multiwfn，因为待会儿还得把算好的Bq点的屏蔽张量输进去。如果关了的话，之后还得重新做上述操作。
我们基于C36的结构创建一个Gaussian输入文件，然后把刚才得到的1.2769174726 -2.7999963248 0.7116924541这个点作为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.2769174726 -2.7999963248 0.7116924541
用gview查看一下，可见Bq原子位置很合适，确实是在环平面上方1埃处
(, 下载次数 Times of downloads: 157)
注意输入文件里的nosymm很重要，如果不写的话，Gaussian可能会自动旋转体系到标准朝向下，导致垂直于环平面的方向在计算前后不同，Multiwfn也就没法正确提取垂直于环平面方向的分量了。
用Gaussian计算此文件，看到NMR部分Bq原子的输出
37 Bq Isotropic = -6.1406 Anisotropy = 25.7238
XX= 0.0370 YX= -0.4791 ZX= -1.2514
XY= 3.5065 YY= -23.1631 ZY= 7.3440
XZ= 7.0901 YZ= -35.5524 ZZ= 4.7043
Eigenvalues: -29.2734 -0.1570 11.0086
切换回Multiwfn窗口，将屏蔽张量依次输入进去：
0.0370,-0.4791,-1.2514
3.5065,-23.1631,7.3440
7.0901,-35.5524,4.7043
（秘籍：想图省事的话，把Gaussian输出文件的屏蔽张量那三行一次性直接粘贴到Multiwfn窗口里也可以）
最后会显示
The shielding value normal to the plane is -12.1124014284
The NICS_ZZ value is thus 12.1124014284
即NICS(1)_ZZ值为12.11（单位是ppm）。
其实本例在进入主功能25的子功能4之前可以不去特意用主功能100的子功能21。因为当前用的环中心就是几何中心，而且定义拟合平面用的原子就是环上的所有原子。对于这种情况，直接进入子功能4，让你输入环中心的时候直接敲回车，之后输入环上的原子序号后，程序自动就会将这些原子的几何中心作为环中心。
3 其它事项
对于平面上方和下方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算一次即可，而不需要每研究一个环都单独算一次。
作者Author: 赵云跳槽 时间: 2015-4-22 09:25
我想问下，平时说的x, y, z方向跟xx,yy,zz方向 有什么可以联系的地方？
作者Author: sobereva 时间: 2015-4-22 14:22
X,Y,Z是方向，而XX,YY,ZZ都是二阶张量的分量。
下图中，B'是感生磁场，B0是外加磁场，σ是磁屏蔽张量。
(, 下载次数 Times of downloads: 190)
根据矩阵运算关系可见，比如σXY就体现了在Y方向外加B0y磁场时，对X方向的感生磁场会有多少贡献
作者Author: haipwei 时间: 2016-7-9 20:17
老师，您好，我想问下除了用Multiwfn还没有没别的确定要计算的Bq原子的位置的方法，因为我下了个Multiwfn不太会用，谢谢老师。
作者Author: sobereva 时间: 2016-7-10 07:09
没有
不会就学啊！谁天生就会。你用别的软件天生就会？其它哪个程序像multiwfn一样手册里又有一大堆例子，还特意写博文去详细示例介绍？
Multiwfn这么简单好用的程序都不会用，其它连个像样英文手册都没有的程序能会用才怪。
作者Author: weizheng007 时间: 2017-10-20 08:49
请问老师几个问题，
(, 下载次数 Times of downloads: 170)
(, 下载次数 Times of downloads: 178)
一、（图一）本征值与屏蔽张量的关系以及与XX、YY、ZZ的关系，二、（图二）G09计算出Bq的Iso值，是不是代表NICS（1）的值， 如果用G09输出的NMR结果如何得到NICS（1）_ZZ的值以及与Multiwfn软件结果的区别
作者Author: sobereva 时间: 2017-10-20 08:52
1 看这里的描述
将核独立化学位移(NICS)分解为sigma和pi轨道的贡献
http://sobereva.com/145
2 如果Bq那个位置是在平面上方1埃处，是NICS(1)的值
其余问题看
衡量芳香性的方法以及在Multiwfn中的计算
http://sobereva.com/176
里关于NICS的描述
作者Author: weizheng007 时间: 2017-10-20 08:58
看了才迷糊的， 1、博文有的XX、YY、ZZ以外的数值为零， 有的不为零，但Iso均为对角元的平均值，Aniso也有讲述，但没有看到本征值讲法，2、Gaussian输出结果和Multiwfn结果不一样更迷糊
作者Author: sobereva 时间: 2017-10-20 17:52
矩阵本征值是线性代数里的基本内容，对磁屏蔽张量对角化一下就完了
正因为倾斜的环你直接从高斯输出文件里读不到垂直于平面的分量，所以才需要用Multiwfn基于高斯的输出数据进行变换来得到
作者Author: 克里斯保 时间: 2019-3-12 14:24
本帖最后由 克里斯保 于 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这个值？是因为初始输入的时候并没有将环调到平面上的原因吧
作者Author: sobereva 时间: 2019-3-12 14:49
纯平面体系，而且平行于XY的时候，直接读ZZ取负值和用Multiwfn给出来的完全一样
本帖重点是非纯平面、倾斜环的NICS_ZZ的计算
作者Author: trevormx 时间: 2019-8-28 09:24
sobereva老师您好，我想知道NICS_ZZ计算的数学原理是什么？我计算一个苯环六个原子不在一个平面上，用间隔的三个原子去定义一个平面也感觉不太科学，我是通过六个原子的坐标最小二乘法来拟合这个平面的。能否将这个平面方程或者它的法向量代入到Multiwfn来计算？或者有没有数学公式，我用matlab之类的软件来计算？
作者Author: sobereva 时间: 2019-8-28 23:25
公式在手册3.100.24节提到了，你拟合了平面方程就有了法矢量，按照公式对磁屏蔽张量做个简单的矩阵运算就能得到
作者Author: sobereva 时间: 2019-8-30 04:29
对Multiwfn程序刚做了更新，现已可以在官网下载。之前版本定义环平面是通过三个原子定义的，但这样不准确，有时候很难找合适的三个原子。目前Multiwfn官网上最新的Multiwfn改进了这点，用户可以输入一批原子（比如环上的所有原子），程序会用这些点通过最小二乘法拟合出最有代表性的平面，明显合理多了。
相应地，本文的内容和例子都进行了更新。
作者Author: trevormx 时间: 2019-9-1 19:37
辛苦SOBEREVA老师了！谢谢啦
作者Author: trevormx 时间: 2019-9-1 19:39
非常感谢sobereva老师，公式解释的很明白啊，我自己也用mathematica编了一个程，验证了一下multiwfn，原谅我比较多疑，喜欢明白自己在算什么。
作者Author: QuantumicGuy 时间: 2022-1-12 18:09
本帖最后由 QuantumicGuy 于 2022-1-12 18:50 编辑
请问sob老师，螺烯体系算了多个环平面上下1 A的NICS(1)_ZZ，有差距，请问是取苯环上下的平均值吗
作者Author: sobereva 时间: 2022-1-13 08:19
把立体结构图贴出来，标上环的序号
作者Author: QuantumicGuy 时间: 2022-1-13 09:50
图片左侧贴出的是4号分子的主视图、X、Y方向的分子结构，谢谢老师！
作者Author: sobereva 时间: 2022-1-13 11:08
取平均
作者Author: QuantumicGuy 时间: 2022-1-13 13:14
谢谢sob老师！
作者Author: 乐平 时间: 2022-3-31 11:22
本帖最后由 乐平 于 2022-3-31 05:30 编辑
看了你的结果受到一点启发。
我发现你只标记了一部分的环（1,2,3,4,5，或者 1,2,3,4,5,6），似乎是在立体图左侧的那些。
能不能所有的环同时计算呢（比如立体图左侧，和立体图右侧，它们分别在两个倾斜方向不同的平面）？
作者Author: QuantumicGuy 时间: 2022-4-24 12:49
这个double helicene是有对称性的，算一半就可以了。我这里算的是环的NICS(1)zz、NICS(-1)zz（对应表格里的up、down）
作者Author: 乐平 时间: 2022-4-24 20:15
本帖最后由 乐平 于 2022-4-24 17:15 编辑
感谢回复和讨论！
我有个类似的非平面结构，一直在纠结怎么取平面，看了你的帖子才明白可以取各个小苯环作为小平面来算（之前一直以为必须取结构几何中心的平面）。
另外，我还想问，这些编号 1, 2, 3, 4, 5 的小苯环平面的 Bq 上方以及下方 1 A 距离的 Bq 都在同一个输入文件里一起计算？
再次另外，Sob 老师的原贴里在准备 Gaussian 的 .gjf 输入文件时，写着
建议先别关Multiwfn，因为待会儿还得把算好的Bq点的屏蔽张量输进去。如果关了的话，之后还得重新做上述操作。
对于非平面体系，具有多个苯环的稠环大体系，不可能一直等到计算结束一直不关闭 Multiwfn。既然输入文件里已经有了 Bq 的坐标，Multiwfn 在处理计算结果的时候，能否直接从输出文件（.out）的开头读取 Bq 的坐标信息呢？
以及，
用Gaussian计算此文件，看到NMR部分Bq原子的输出
37 Bq Isotropic = -6.1406 Anisotropy = 25.7238
XX= 0.0370 YX= -0.4791 ZX= -1.2514
XY= 3.5065 YY= -23.1631 ZY= 7.3440
XZ= 7.0901 YZ= -35.5524 ZZ= 4.7043
Eigenvalues: -29.2734 -0.1570 11.0086
切换回Multiwfn窗口，将屏蔽张量依次输入进去：
0.0370,-0.4791,-1.2514
3.5065,-23.1631,7.3440
7.0901,-35.5524,4.7043
（秘籍：想图省事的话，把Gaussian输出文件的屏蔽张量那三行一次性直接粘贴到Multiwfn窗口里也可以）
Multiwfn 能否直接从输出文件（.out）里读取与之对应的 Bq 的屏蔽张量信息呢？ 对于稠环的大体系，一个一个手动输入很容易看错或遗漏……
作者Author: sobereva 时间: 2022-4-25 02:25
计算任务可能同时涉及多个Bq，原理上程序难以总是自动正确判断哪个Bq对应哪个环，所以还是得用户自己输入
需要反复输入进Multiwfn里的一批命令可以记在一个文本文件里，每次直接粘贴进窗口里
作者Author: 乐平 时间: 2022-4-25 10:14
本帖最后由 乐平 于 2022-4-25 04:19 编辑
谢谢 Sob 老师的回复。
我还是有疑惑
计算任务可能同时涉及多个Bq，原理上程序难以总是自动正确判断哪个Bq对应哪个环，所以还是得用户自己输入
既然可以同时计算多个 Bq，那么 Gaussian 的输出文件（.out）里的结果按理说应该是依照输入文件（.gjf）里 Bq 的顺序依次输出的。
如果不是依次输出的话，我在同一个输入文件 (.gjf) 里设置多个 Bq 的坐标并提交了计算任务，得到结果之后，再手动粘贴屏蔽张量也不可能一一对应上 Bq 呀。换句话说，如果输出文件（.out）里的屏蔽张量不是按照输入文件 (.gjf) 里设置多个 Bq 的顺序输出的话，我怎么知道 .out 里的第一组屏蔽张量对应的是我设置的第一个环的屏蔽张量呢？
作者Author: QuantumicGuy 时间: 2022-4-25 20:13
本帖最后由 QuantumicGuy 于 2022-4-25 20:14 编辑
可以同时算多个，ICSS就是计算一大堆Bq原子(通过Multiwfn绘制等化学屏蔽表面(ICSS)研究芳香性 http://sobereva.com/216)。写好输入文件等算完拿回来用Multiwfn处理。可以用可视化软件（如GaussView）显示原子序号lable，输出文件里Bq原子序号与之对应，这样就清楚的知道算的是哪个位置的NICS了。
作者Author: sobereva 时间: 2022-4-26 01:22
是依次输出的，但是很难确保用户不会做了什么额外的添油加醋的事情（解答大量用户的提问时，发现老有用户做一些匪夷所思的事情）
我觉得手动输入也没什么费劲的，这样也杜绝了额外的可能想象不到的隐患
作者Author: wocaishiniu 时间: 2025-3-31 18:08
sob老师好，我新下载的Multiwfn，但是输入100后，只显示到23，输入25没有反应（如图一），然后我用24算的平面上下一A的鬼原子的坐标，算完以后用gaussview打开，发现鬼原子之间有连接，七元环上下1A的鬼原子和七元环上的碳原子没有连接（如图二画圈位置），您知道这是怎么回事吗？
作者Author: sobereva 时间: 2025-4-1 00:07
25是指主菜单里的主功能25
GaussView根据原子间距离判断成键，不用管
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |