将核独立化学位移(NICS)分解为sigma和pi轨道的贡献
2023-Jun-3注：笔者后来又写了《基于Gaussian的NMR=CSGT任务得到各个轨道对NICS贡献的方法》（http://sobereva.com/670）介绍另一种分解NICS的方法，可以得到各个分子轨道的贡献，只需要用Gaussian而不依赖于NBO
将核独立化学位移(NICS)分解为sigma和pi轨道的贡献
Decomposition of nuclear independent chemical shift (NICS) into contributions from sigma and pi orbitals
文/Sobereva @北京科音
First release: 2012-May-31 Last update: 2014-Aug-1
1 前言
衡量芳香性的指标很多，也没有统一的定义。现今存在几十种指标，从磁性质、电子离域性、能量、结构等诸多方面对芳香性进行描述。NICS是被使用得最广泛的衡量芳香性的指标。它的含义是在某个人为设定的不在原子核位置上的磁屏蔽值的负值，越负（对磁场屏蔽越强）则芳香性越强。这个位置设定得有一定任意性。原文(JACS,118,6317)是取在共轭环的几何中心，为了明确起见，后来被特称为NICS(0)。有人认为取在环中心会将sigma和pi轨道的贡献叠加在一起而说不清楚，于是提出取在平面上方或下方1埃的位置，称为NICS(1)，这个指标体现的主要是pi电子的贡献。后来有人更进一步认为，NICS值不该取各项同性值（磁屏蔽张量XX、YY、ZZ之和的平均），而应该只考察ZZ值（Z垂直于环平面），这样最适合体现pi芳香性，这被称为NICS(1)_ZZ。虽然NICS(1)_ZZ很适合讨论pi芳香性，但是还是稍微掺入了一点（尽管可忽略）的sigma轨道的影响，而且也不可能通过选取某个空间位置的来单独讨论sigma芳香性。
磁屏蔽值可以是分解为轨道的贡献的，因此，直接研究各个sigma轨道和各个pi轨道对NICS值的贡献就可以直接研究sigma和pi芳香性。然而一般量子化学程序都不支持分解为轨道的贡献，虽然一些文献中使用DeMon-Master程序去做，但是此程序貌似是私下传播，难以获取（注意这并不是公开的DeMon2k，DeMon2k是做不到的）。相对明朗一些的方法是使用NBO程序的NCS（自然化学屏蔽）功能来做，它既可以将磁屏蔽张量分解为定域化轨道的贡献，也可以分解为正则分子轨道(CMO)的贡献。
实际上，NBO程序3.1版及后续版本的NCS功能有很大不同。NBO 3.1能做的是将磁屏蔽张量分解为NLMO的贡献，而不能分解为NBO或MO的贡献。而从后续版本，比如常用的5.0版，能做的是将磁屏蔽张量分解为NBO或MO的贡献，据我所知，反倒不能再分解为NLMO的贡献了。
能做NCS的NBO程序必须是内嵌到Gaussian作为l607模块的版本，独立运行的NBO不行。在第二节，将介绍如何用NBO5.0版l607做NCS，这也被称为NBO5.G。NBO5.G是收费的，但网上能找到源码，可以自行编译，方法见《编译NBO5.0独立运行版和嵌入Gaussian03 C02版的方法》（http://sobereva.com/144）。在第三节将介绍如何直接用Gaussian自带的NBO3.1版的L607做NCS。
2 使用NBO 5.0分解磁屏蔽张量
NBO5.0和G03的C及以后版本之间有点兼容性问题。也就是，直接通过Bq虚原子定义NICS计算的位置的话，NBO5.0发现Bq原子没有S基函数，就会报错终止。为解决这个问题，要么用G03的A或B版本，要么用更新版本的NBO5.x程序。虽说可以把Bq写成比如H-Bq，这样的话这个虚原子上就有了当前基组中H的基函数，可以避免报错，但这额外的基函数会明显影响NICS值。更好的办法是给这个虚原子只赋予一个弥散范围很大的S基函数，这个基函数基本不会对体系能量、NICS值带来什么影响，却能够由此绕过这个不兼容问题。
本文都用B3LYP/6-31+G*优化的苯作为例子，NICS、NCS计算也是在这个级别下进行。输入文件如下，Z方向垂直于苯环平面
#P b3lyp/gen NMR pop=nboread //虽然相关资料说还需要写IOp(10/46=1)，但笔者发现没必要
b3lyp/6-31+g(d) opted
0 1
C 0.00000000 1.39864200 0.00000000
C 1.21125900 0.69932100 0.00000000
C 1.21125900 -0.69932100 0.00000000
C 0.00000000 -1.39864200 0.00000000
C -1.21125900 -0.69932100 0.00000000
C -1.21125900 0.69932100 0.00000000
H 0.00000000 2.48606800 0.00000000
H 2.15299800 1.24303400 0.00000000
H 2.15299800 -1.24303400 0.00000000
H 0.00000000 -2.48606800 0.00000000
H -2.15299800 -1.24303400 0.00000000
H -2.15299800 1.24303400 0.00000000
bq 0. 0. 0. //环中心
bq 0. 0. 1. //环平面中心上方1埃
C H 0
6-31+g(d)
****
13 14 0 //设定13、14号原子，也就是虚原子的基函数
S 1 1.0 //都只含有一个收缩度为1的S型基函数
0.005 1. //第一个数值越小，表明这个基函数越弥散，它的存在对结果影响也会越小，0.005已经是十分弥散了。如果设得太小会导致报错。第二个数是收缩系数，对于当前基函数收缩度为1的情况这个值是随意的，一般就设1.0。
****
$NBO NCS <XYZ MO> $END // $NBO和$END夹着的是传递给NBO模块的关键词。被<>括住的关键词是NCS分析的具体选项
Gaussian输出的两个虚原子的磁屏蔽张量信息为
13 Bq Isotropic = 7.8997 Anisotropy = 7.6878
XX= 5.3444 YX= 0.0000 ZX= 0.0000
XY= 0.0000 YY= 5.3298 ZY= 0.0000
XZ= 0.0000 YZ= 0.0000 ZZ= 13.0249
Eigenvalues: 5.3298 5.3444 13.0249
14 Bq Isotropic = 10.0608 Anisotropy = 27.9171
XX= 0.7644 YX= 0.0000 ZX= 0.0000
XY= 0.0000 YY= 0.7457 ZY= 0.0000
XZ= 0.0000 YZ= 0.0000 ZZ= 28.6722
Eigenvalues: 0.7457 0.7644 28.6722
这表明NICS(0)为-7.8997，NICS(0)_ZZ=-13.0249。NICS(1)为-10.0608，NICS(1)_ZZ=-28.6722。可见，无论是在环中心还是其上方1埃处，环上离域电子对垂直于环平面的磁场屏蔽强度比起平行于环平面的要明显大得多。
Isotropic（各向同性屏蔽值）就是一般所说的绝对化学位移，单位是ppm，和溶液中的普通核磁共振实验数据是对应的（相对于混乱运动的分子，外磁场相当于没有特定方向，故视为各向同性），这个值是磁屏蔽张量矩阵的迹的平均，比如7.8997=(5.3444+5.3298+13.0249)/3。Anisotropy用来衡量各项异性的程度，计算方法是Eig(3) - (Eig(2) + Eig(1))*0.5，Eig(i)代表第i个本征值，从小到大排序。如果屏蔽张量在各方向相同，此值即为0。
为了具体了解各个NBO及其离域效应对屏蔽值的影响，需要在NCS关键词后的<>内写上XYZ关键词。这样就会依次输出每个原子的磁屏蔽张量的来自各个NBO的贡献。比如对于环中心的虚原子：
NBO
L NL XX YX ZX XY YY ZY XZ YZ ZZ
===============================================================================
1. -4.16 2.65 0.00 2.65 -1.10 0.00 0.00 0.00 -3.73
50. -0.95 0.30 0.00 -0.43 0.14 0.00 0.00 0.00 -1.52
51. 0.66 -0.22 0.00 0.12 -0.17 0.00 0.00 0.00 0.88
54. 0.33 -0.14 0.00 0.40 -0.21 0.00 0.00 0.00 0.23
55. 0.15 -0.05 0.00 0.03 -0.03 0.00 0.00 0.00 0.13
56. 0.39 -0.21 0.00 0.57 -0.32 0.00 0.00 0.00 0.15
......略
21. -0.08 -0.51 0.00 -0.51 -0.67 0.00 0.00 0.00 -1.76
24. 0.00 0.00 0.00 0.02 0.04 0.00 0.00 0.00 0.15
28. -0.01 -0.03 0.00 0.03 0.05 0.00 0.00 0.00 0.11
76. 0.02 0.04 0.00 0.01 0.02 0.00 0.00 0.00 0.15
80. 0.03 0.05 0.00 0.00 0.00 0.00 0.00 0.00 0.11
-------------------------------------------------------------------------------
L -3.84 -0.02 0.00 -0.01 -3.85 0.00 0.00 0.00 -22.66
NL 9.25 -0.02 0.00 -0.03 9.24 0.00 0.00 0.00 35.69
-------------------------------------------------------------------------------
Total 5.41 -0.04 0.00 -0.04 5.40 0.00 0.00 0.00 13.03
-------------------------------------------------------------------------------
第一列是NBO的编号，如果靠左边写就说明它是占据数较高的NBO轨道，也被称为Lewis型(L)。靠右写的是非Lewis型(NL)，绝大多数都是占据数很低的NBO，但也可以是有一定占据但不那么接近2的NBO。比如XX标签下面的第三个数0.66，可以解释为1号NBO向51号NBO离域的微量电子对屏蔽张量XX分量的影响为0.66ppm（更确切来说，是51号NBO对应的NLMO由1号NBO引入的“离域性尾巴”的贡献）。最后Total显示的总XX/YY/ZZ值5.41/5.40/13.03和Gaussian前面给出的XX/YY/ZZ方向屏蔽值5.3444/5.3298/13.0249是十分吻合的。Total当中L就是所有L项的加和，NL就是所有NL项的加和。注意只有数值大于特定阈值的项才列在表里，阈值可以通过比如NCS=0.02关键词来设定为0.02。仔细观看ZZ列，并结合轨道图形观察，会发现三个C-C pi型NBO对磁场有很强屏蔽作用，每个贡献都高达9.88ppm。而C-C和C-H的sigma键都会产生一定的去屏蔽作用，每个键都贡献负三点几。
紧接着输出的是这个原子屏蔽张量主轴方向的值，33对应整体屏蔽磁场最强的方向，11对应最弱的方向。（如果表中有任何一项的12、13或23方向的值不很小，那么非对角元也会被输出出来，不过此例所有项的非对角元都很小所以没显示）。主轴与笛卡尔坐标的变换矩阵也会输出。ISO就是各向同性屏蔽值，CSA就是各向异性值，和前面给出的计算式子一样，比如第一行，-1.09=-3.73-0.5*(0.24-5.50)。由于列数比较少，所以每行开头都留出了空间直接给出了NBO编号对应的含义。如NL 54. C 3(ry*)就代表这是54号NBO，是3号碳上的里德堡轨道。
NBO 11 22 33 CSA ISO
============================================================
1. C 1- C 2 0.24 -5.50 -3.73 -1.09 -3.00
NL 50. C 3(ry*) -0.37 -0.43 -1.52 -1.11 -0.77
NL 51. C 3(ry*) 0.12 0.36 0.88 0.64 0.46
NL 54. C 3(ry*) 0.14 -0.02 0.23 0.17 0.11
NL 55. C 3(ry*) 0.03 0.09 0.13 0.07 0.08
NL 56. C 3(ry*) 0.15 -0.08 0.15 0.11 0.07
......略
20. C 5(cr) 0.08 -0.83 -1.76 -1.38 -0.84
NL 63. C 4(ry*) 0.01 0.03 0.15 0.12 0.06
NL 67. C 4(ry*) 0.02 0.01 0.11 0.10 0.05
NL 89. C 6(ry*) 0.00 0.05 0.15 0.12 0.06
NL 93. C 6(ry*) -0.01 0.04 0.11 0.10 0.05
21. C 6(cr) -0.94 0.18 -1.76 -1.38 -0.84
NL 24. C 1(ry*) 0.04 0.01 0.15 0.12 0.06
NL 28. C 1(ry*) 0.02 0.01 0.11 0.10 0.05
NL 76. C 5(ry*) 0.05 0.00 0.15 0.12 0.06
NL 80. C 5(ry*) 0.04 -0.01 0.11 0.10 0.05
------------------------------------------------------------
Lewis -3.86 -3.83 -22.66 -18.82 -10.12
non-Lewis 9.23 9.27 35.69 26.44 18.06
------------------------------------------------------------
Total 5.37 5.44 13.03 7.62 7.95
------------------------------------------------------------
最后NCS分析会输出以下汇总内容，如果不写<XYZ>的话NCS分析也只会输出这些而不会输出上述每个原子的成分细节
Summary of isotropic NMR chemical shielding
Total Lewis (L) and non-Lewis (NL) contributions: (ppm)
NBO C( 1) C( 2) C( 3) C( 4) C( 5) C( 6) H( 7)
------------ ------- ------- ------- ------- ------- ------- -------
1. C 1- C 2 L -49.21 -49.22 -4.75 -0.06 -0.06 -4.75 -1.67
NL -2.91 -2.90 3.56 -0.07 -0.08 3.56 -0.26
2. C 1- C 2 L 0.27 0.28 3.24 1.61 1.61 3.25 0.58
NL 1.47 1.47 -3.04 2.87 2.88 -3.07 0.18
3. C 1- C 6 L -49.21 -4.75 -0.06 -0.06 -4.75 -49.22 -1.67
NL -2.91 3.56 -0.08 -0.07 3.56 -2.90 -0.26
4. C 1- H 7 L -35.94 -4.61 -0.82 -0.45 -0.82 -4.61 28.45
NL 7.33 2.22 0.49 0.19 0.49 2.22 -1.23
......略
20. C 5(cr) L -0.26 -0.21 -0.26 -0.41 203.56 -0.41 -0.10
NL 0.08 0.05 0.08 0.05 0.11 0.05 0.02
21. C 6(cr) L -0.41 -0.26 -0.21 -0.26 -0.41 203.56 -0.02
NL 0.05 0.08 0.05 0.08 0.05 0.11 0.00
------------ ------- ------- ------- ------- ------- ------- -------
Lewis 51.84 51.86 51.91 51.87 51.91 51.93 26.16
non-Lewis 15.86 15.89 15.87 15.84 15.84 15.84 -1.38
------------ ------- ------- ------- ------- ------- ------- -------
Total 67.69 67.75 67.78 67.71 67.76 67.77 24.77
这里C1-C2出现两次是因为它们之间被认为是双键。
如果在NCS后面的<>里写上MO，则会输出每个原子屏蔽张量来自于每个占据正则分子轨道的贡献，如下所示。还会紧跟着输出在主轴方向下的值，最后也会输出汇总信息。
Full Cartesian NMR shielding tensor (ppm) for atom gh(13):
Canonical MO contributions
MO XX YZ ZX XY YY ZY XZ YZ ZZ
===============================================================================
1. 0.12 0.00 0.00 0.00 0.12 0.00 0.00 0.00 0.22
2. 0.12 0.00 0.00 0.00 0.12 0.00 0.00 0.00 -1.20
3. 0.12 0.00 0.00 0.00 0.12 0.00 0.00 0.00 -1.20
4. 0.07 0.00 0.00 0.00 0.22 0.00 0.00 0.00 -2.99
5. 0.22 0.00 0.00 0.00 0.07 0.00 0.00 0.00 -2.99
......略
19. 1.27 0.00 0.00 0.00 8.03 0.00 0.00 0.00 21.90
20. -0.44 0.00 0.00 0.00 4.87 0.00 0.00 0.00 -5.07
21. 4.85 0.00 0.00 0.00 -0.45 0.00 0.00 0.00 -5.08
-------------------------------------------------------------------------------
Total 41.91 0.00 0.00 0.00 41.93 0.00 0.00 0.00 70.75
-------------------------------------------------------------------------------
遗憾的是NBO5.0结合G03 C.02版给出的结果我发现是明显错误的，根本没法用，总的XX、YY、ZZ值明显和Gaussian算出来的以及分解为NBO贡献时的总值不符。结合G98得到的结果肯定是正确的，如果是G03/09的用户，可能得用更新的NBO版本才行，但官方主页上宣称只支持G9X。
3 使用NBO 3.1分解磁屏蔽张量
在G03/09中只要用了NMR关键词，并且同时用了比如pop=nbo关键词调用了NBO 3.1分析模块，则NBO 3.1模块就会自动进行NCS分析，例如：
#P b3lyp/6-31+G* NMR pop=nbo
b3lyp/6-31+g(d) opted
0 1
C 0.00000000 1.39864200 0.00000000
C 1.21125900 0.69932100 0.00000000
C 1.21125900 -0.69932100 0.00000000
C 0.00000000 -1.39864200 0.00000000
C -1.21125900 -0.69932100 0.00000000
C -1.21125900 0.69932100 0.00000000
H 0.00000000 2.48606800 0.00000000
H 2.15299800 1.24303400 0.00000000
H 2.15299800 -1.24303400 0.00000000
H 0.00000000 -2.48606800 0.00000000
H -2.15299800 -1.24303400 0.00000000
H -2.15299800 1.24303400 0.00000000
bq 0. 0. 0.
bq 0. 0. 1.
这里不必像NBO5.0那样用混合基组的技巧解决与Gaussian兼容性问题，NBO3.1对虚原子的处理和Gaussian是兼容的。
例如，在NBO模块输出信息末尾可以找到各个NLMO对苯环中心的虚原子的磁屏蔽张量：
Contributions to Isotropic shielding for atom 13 ( 7.9768 ppm):
Standard Origin Atomic Origin
NLMO Total Paramagnetic Diamagnetic Paramagnetic Diamagnetic
1 -1.3247 -3.4073 2.0825 -1.3247 0.0000
2 -1.3248 -3.4072 2.0825 -1.3248 0.0000
3 8.1896 3.0323 5.1573 8.1896 0.0000
4 -0.8489 -1.4446 0.5957 -0.8489 0.0000
5 -1.3247 -3.4080 2.0833 -1.3247 0.0000
6 8.1936 3.0366 5.1570 8.1936 0.0000
7 -0.8492 -1.4451 0.5959 -0.8492 0.0000
8 -1.3247 -3.4073 2.0825 -1.3247 0.0000
9 -0.8492 -1.4451 0.5959 -0.8492 0.0000
10 -1.3248 -3.4072 2.0825 -1.3248 0.0000
11 8.1896 3.0323 5.1573 8.1896 0.0000
12 -0.8489 -1.4446 0.5957 -0.8489 0.0000
13 -1.3246 -3.4080 2.0835 -1.3246 0.0000
14 -0.8492 -1.4451 0.5959 -0.8492 0.0000
15 -0.8492 -1.4451 0.5959 -0.8492 0.0000
Sum 11.5298 -20.0134 31.5433 11.5298 0.0000
顺磁和抗磁成份这里也单独给了出来，但是一般只需要关注它们的总效应，即Total那列就行了。
可见这里Sum的值11.5298ppm和Gaussian算出来的各向同性磁屏蔽值7.9768对应不上，这是因为默认情况下，贡献较小的NLMO不会被输出出来，而Sum值仅仅是所有全部显示出来的贡献值的加和，这里只显示出了全部21个NLMO中的15个的贡献。由于没有算进去那些没显示出来的，所以总值自然对应不上。这个默认设定非常傻。为了让所有NLMO的贡献都显示出来，需要加上IOp(6/65=-1)。
默认情况下输出的只是各向同性磁屏蔽值来自各NLMO的贡献，如果想要将NLMO对磁屏蔽张量的各个元素的贡献也输出出来，就要用pop=NCSall关键词，而pop=nbo就不必写了。例如上面例子的route section写成
#P b3lyp/6-31+G* NMR pop=NCSall IOp(6/65=-1)
在输出文件中就能找到诸如下面这样的信息，全部NLMO对环中心的磁屏蔽张量ZZ分量的贡献都给出了。各个NLMO对NICS(0)_ZZ的贡献也就是Total那列的负值。
Contributions to Atom 13 Sigma-ZZ ( 13.1863 ppm):
Standard Origin Atomic Origin
NLMO Total Paramagnetic Diamagnetic Paramagnetic Diamagnetic
1 -0.8999 -5.9092 5.0093 -0.8999 0.0000
2 12.0045 5.7848 6.2197 12.0045 0.0000
3 -0.8999 -5.9092 5.0093 -0.8999 0.0000
4 -1.6607 -1.3298 -0.3309 -1.6607 0.0000
5 -0.9002 -5.9111 5.0108 -0.9002 0.0000
6 -1.6618 -1.3315 -0.3303 -1.6618 0.0000
7 -0.8999 -5.9092 5.0093 -0.8999 0.0000
8 12.0045 5.7848 6.2197 12.0045 0.0000
9 -1.6618 -1.3315 -0.3303 -1.6618 0.0000
10 -0.8999 -5.9092 5.0093 -0.8999 0.0000
11 -1.6607 -1.3298 -0.3309 -1.6607 0.0000
12 -0.9003 -5.9111 5.0108 -0.9003 0.0000
13 12.0057 5.7863 6.2194 12.0057 0.0000
14 -1.6618 -1.3315 -0.3303 -1.6618 0.0000
15 -1.6618 -1.3315 -0.3303 -1.6618 0.0000
16 -1.2438 -0.7266 -0.5172 -1.2438 0.0000
17 -1.2431 -0.7258 -0.5172 -1.2431 0.0000
18 -1.2431 -0.7258 -0.5172 -1.2431 0.0000
19 -1.2438 -0.7266 -0.5172 -1.2438 0.0000
20 -1.2431 -0.7258 -0.5172 -1.2431 0.0000
21 -1.2431 -0.7258 -0.5172 -1.2431 0.0000
Sum 13.1863 -30.4450 43.6313 13.1863 0.0000
可以按照《使用Multiwfn绘制NBO及相关轨道》（http://sobereva.com/134）来观看这些NLMO的图形，找出sigma和pi轨道。最好不要用pop=saveNLMO将NLMO保存到check文件中去观看，而建议用plot文件方式来观看，以免自动对轨道重新排序而不便与NCS输出信息对应。通过图形观察可知2、8、13号NLMO都是C-C pi轨道，磁屏蔽作用很强，每个都达到12ppm，因此NICS(0)_ZZ_pi值即是3*-12=-36ppm。有兴趣的读者可以参看Org. Lett., 8, 863文章的补充材料的末尾，作者是利用NBO5.G结合G98将磁屏蔽张量分解为CMO的贡献并计算的NICS(0)_ZZ_pi，他的结果和本文的值相符很好。
总的来说，用Gaussian自带的NBO3.1来做NCS比起用NBO5.0来做明显更方便，而且也不用买NBO5.0，也不涉及Gaussian源代码版权问题。