- Dec 24, 2021 更新关于耦合常数计算的介绍
- Jul 16, 2019 初版
本文简单介绍一下怎么用高斯的DFT计算NMR化学位移和耦合常数，以及如何在GaussView里观看计算生成的图谱。计算和可视化所用的软件分别是G16 A.03, GV 6.0.16。
用量子化学方法计算NMR化学位移现在已经有了常用套路，DFT泛函一般选取B97-2、revTPSS（注：高斯关键词分别写作B972、revTPSSrevTPSS），基组常用pcSseg-1或pcSseg-2，更换其他泛函（如B3LYP、 M062X）或基组（如6-311G**、def2-TZVP），结果精度在统计意义上是不会比上述常用搭配要好的。当然，如果体系很小或土豪机器配置很高，用得起cc-pVQZ基组，那自然是用QZ更好。
标准物质TMS的计算
与实验上的1H-NMR化学位移需要用到标准物质四甲基硅烷（TMS）一样，理论计算上也需要分别计算TMS与目标分子中氢的核磁屏蔽数值，然后相减就是目标分子中氢的核磁化学位移。下面展示的是TMS分子的高斯输入文件：
代码语言：javascript复制%chk=TMS_def2TZVP_opt.chk
%mem=16GB
%nprocshared=28
#p M062X/def2TZVP opt=calcfc freq scrf=(solvent=Chloroform)
[空行]
geometry opt of TMS at Td point group [标题行，可瞎写]
[空行]
0 1
Si 0.00000000 0.00000000 0.00000000
C 0.00000000 0.00000000 1.94000006
H -0.87364980 -0.50440195 2.29666647
H 0.87364980 -0.50440195 2.29666647
H 0.00000000 1.00880390 2.29666647
C -0.00000000 -1.82904960 -0.64666669
H -0.87364980 -2.33345189 -0.29000077
H -0.00000000 -1.82904994 -1.71666493
H 0.87364980 -2.33345189 -0.29000077
C 1.58400341 0.91452480 -0.64666669
H 1.58400372 0.91452497 -1.71666493
H 1.58400372 1.92332887 -0.29000077
H 2.45765352 0.41012302 -0.29000077
C -1.58400341 0.91452480 -0.64666669
H -1.58400372 1.92332887 -0.29000077
H -1.58400372 0.91452497 -1.71666493
H -2.45765352 0.41012302 -0.29000077
[空行]
--Link1--
%chk=TMS_def2TZVP_opt.chk
%mem=16GB
%nprocshared=28
#p B972/gen NMR guess=read geom=allcheck scrf=(solvent=Chloroform)
[空行]
[自行去EMSL网站上拷贝pcSseg-2基组数据，贴于此]在GV里构建TMS分子结构后，可以利用右键Tools -> Point Group工具，然后勾选Enable Point Group Symmetry，点击右中区域Symmetrize，使结构符合Td点群对称性，这样所有氢是核磁化学等价的，算出来的核磁屏蔽数值相同。
从上面的输入文件可以看出，先用M062X/def2TZVP进行结构优化，后用B972/pcSseg-2算NMR化学位移。计算核磁的关键词为NMR，相当于NMR=GIAO，表示用默认的GIAO方法计算。两步任务所使用的泛函和基组不要求一致。这有两点原因：（1）几何结构的描述上B97-2泛函未必是最好的；（2）pcSseg-2基组挺大的（比cc-pVTZ还大），拿它算单点还行，做结构优化时间就比较长了。如果将来研究的体系比较大，无论是TMS还是目标分子，基组建议统一降为pcSseg-1。隐式溶剂化模型的溶剂用的是氯仿，相当于实验上的氘代氯仿。若不用隐式溶剂化模型，气相中的结果一般也还行。
算完以后用文本编辑器打开log（或out）文件，搜索Isotropic，会发现每个原子有一处，每处等号后面的小数就是对应原子的核磁屏蔽数值，氢、碳、硅全都有（见下图），这里氢的核磁屏蔽各向同性数值为31.6753。Linux下可以直接用命令grep 'H Isotropic =' xx.log搜索、打印到屏幕上。
目标分子的计算
接下来我们以丙酮分子为例，计算其1H-NMR化学位移。丙酮分子的输入文件与TMS类似，这里就不贴出了。同样，在构造丙酮分子之后可以选择点群工具使结构对称化为C2v点群；不做对称化、直接优化也行。要注意的是，对称化这一步并不能保证最后的优化结果被高斯识别为C2v点群，可能会有微小的偏差；严格的C2v点群需要使用内坐标定义分子结构。上面的TMS分子优化后仍保持Td点群是凑巧的。下图展示的是从丙酮输出文件grep出的结果：
代码语言：javascript复制3 H Isotropic = 29.3619 Anisotropy = 4.8253
4 H Isotropic = 29.2212 Anisotropy = 4.2041
5 H Isotropic = 29.8594 Anisotropy = 6.6485
7 H Isotropic = 29.3618 Anisotropy = 4.8249
8 H Isotropic = 29.2211 Anisotropy = 4.2044
9 H Isotropic = 29.8594 Anisotropy = 6.6487可以看出左侧3个氢（3,4,5）与右侧3个氢（7,8,9）有十分微小的差异，这是略微偏离C2v点群造成的，影响可以忽略。实验测得甲基氢的峰只有一个，但是计算得到的3,4,5三个氢的数值却是不一样的。这是由于我们计算的只是一个静态结构，室温（至0 ℃范围均可）时丙酮分子中的甲基是可以自由旋转的，因此实验测得的是大量不同构象（符合玻尔兹曼分布）下的统计平均结果。如果是凝聚相下大分子中化学位移的计算，是要通过跑MD从轨迹中选取帧数进行统计平均的；不过对于常见有机小分子，当精度要求不是很高的时候，我们只要简单地把优化得到的结构中的三个氢平均一下再与TMS进行比较即可：
这相当于氢原子之间交换序号得到的local minima的构象平均。实验上丙酮在氘代氯仿溶剂中的化学位移是多少呢？笔者不会做实验，只好去找数据库了，比如AIST的SDBS有机分子谱学数据库
https://sdbs.db.aist.go.jp/sdbs/cgi-bin/cre_index.cgi
常见小分子各种谱的数据挺详细的。同意网站搜索协议之后，搜索acetone
接下来选择对应的细分类、溶剂、氢谱，可以看到实验值为2.162 ppm。说明上述计算结果还是挺接近的。按照相应的流程，也可以计算13C-NMR碳谱、19F-NMR氟谱等等，只要注意对应的标准物分别是什么就行了。不过非氢元素的化学位移值本身往往很大，因此计算与实验相差2 ppm左右也是正常的，DFT结合隐式溶剂化模型不能苛求更高的精度了。
用GV看NMR图谱
接下来介绍一下怎么在GV里看NMR图谱（这纯粹是一种直观的图示功能，有的同学直接把化学位移值标在相应的键线式原子旁也很直观）。进入GV安装目录的data文件夹，找到nmr.data文件用文本编辑器打开（在修改前可以先备份之），可以看到里面已经预存了一些标准物的NMR屏蔽值，但都是过时的方法或基组，不符合现今需要。我们可以在底下添加一些常用的标准物质
双引号里面是注释，记录是在什么水平下算的，以防日后忘记。修改之后保存，重启GV，打开丙酮的log文件，面板上选择Results -> NMR，会弹出NMR谱图对话框，依次选择Element: H, Reference: TMS B972/pcSseg-2 ultrafineGIAO （对应刚添加的nmr.data文件内容）。
横坐标表示NMR化学位移，纵坐标是简并度，即表示有几个相同的氢。上图可以看出丙酮左右两边甲基的氢是对称的；同时由于我们算的是静态结构，同一个甲基上的三个氢不等价。如果想显示等价的图示效果，可以用平均值将log文件内6个氢的不同数值替换掉，再打开就会显示等价了。GV的NMR谱图对话框Plots菜单中还有导出图片、导出数据及在对应原子上显示化学位移等各种选项，这里不就一一介绍了。
耦合常数的计算
至于耦合常数的计算，只需要加上spinspin子关键词，即NMR=spinspin。不过要注意的是，耦合常数计算时间是频率计算时间的大概两倍[1]，计算量较大。在高斯的输出文件中，最后会以矩阵的形式，输出体系中任意两个原子之间的耦合常数，以Total nuclear spin-spin coupling J (Hz)开头，如下所示：
截图中D是一种Fortran双精度浮点数的格式，相当于Python和C语言中的E，表示乘以10的多少次方的意思。红框处的数据即为2号原子与3号原子之间的耦合常数。在总耦合常数的上方，还会输出耦合常数的四个组成部分，依次为
- Fermi Contact (FC) contribution to J
- Spin-dipolar (SD) contribution to J
- Paramagnetic spin-orbit (PSO) contribution to J
- Diamagnetic spin-orbit (DSO) contribution to J
自旋-自旋耦合常数的计算还需要注意基组的选择。 由于在一般化学问题中，我们往往比较关心价电子的性质，因此量子化学中大部分基组都是对价电子进行充分优化的。而核磁性质与核附近的电子密切相关，尤其是上述计算中的Fermi Contact项，因此需要使用专门为此优化过的基组。推荐使用Frank Jensen开发的pcJ-n系列基组（n可取0~4，所含极化函数依次增多），一般可用pcJ-1基组。
以下我们计算1,1-二氯乙烷的甲基氢（2、3、4号）与α-氢之间的耦合常数，分子结构如下：
输入文件如下：
代码语言：javascript复制%oldchk=C2H4Cl2.chk
%chk=C2H4Cl2-pcJ.chk
#p b972/gen nmr=spinspin geom=allcheck guess=read scrf(solvent=chloroform)
[此处略去pcJ-1基组的数据]其中，C2H4Cl2.chk为几何结构优化的chk文件。
计算得到的结果如下：
代码语言：javascript复制Total nuclear spin-spin coupling J (Hz):
1 2 3 4 5
1 0.000000D 00
2 0.123224D 03 0.000000D 00
3 0.123226D 03 -0.127554D 02 0.000000D 00
4 0.126772D 03 -0.123457D 02 -0.123470D 02 0.000000D 00
5 0.384660D 02 -0.735589D 00 -0.740798D 00 -0.871558D 01 0.000000D 00
6 0.353153D 00 0.307238D 01 0.306465D 01 0.122549D 02 0.172186D 03
7 0.238009D 01 0.913651D 00 0.726357D 01 0.171017D 01 -0.245092D 02
8 0.238029D 01 0.726520D 01 0.917600D 00 0.170527D 01 -0.245084D 02
6 7 8
6 0.000000D 00
7 0.532547D 01 0.000000D 00
8 0.532486D 01 0.469992D 01 0.000000D 00将(6,2)、(6,3)、 (6,4)三个矩阵元取平均值，得6.1。从SDBS数据库查得的实验值为6.0，非常吻合。
在高斯中，若使用普通的基组（如def2-TZVP、cc-pVTZ等），可以使用nmr=mixed关键词，来改善FC项的计算。此时高斯内部会对用户所给基组进去收缩，并加上紧缩的基函数（体现在基函数的轨道指数很大）来描述内核电子。而对其余三项，依然使用用户设定的原始基函数来进行计算。
总结
在结构优化之后，使用常见套路B972/pcSseg-2和B972/pcJ-1计算NMR化学位移和耦合常数，结果与实验值吻合不错。
PS1：标准物质的NMR屏蔽数值显然只需计算一次即可，以后对同一计算级别的目标分子可以重复使用。
PS2：如果是量化新手，不妨练习一下如何找到丙酮中甲基转动的过渡态（在优化好的结构上旋转一个甲基60°即可作为过渡态初猜，笔者算的Gibbs自由能垒是0.07 kcal/mol）。
PS3：计算核磁时需要分别引用泛函B97-2（见高斯官网http://gaussian.com/dft）、基组pcSseg-2和pcJ-1（在EMSL上打开基组数据时就有引用文献），及此搭配计算NMR化学位移的文章（任一篇高水平文章，用于支持这么做的理由）。
PS4：ChemDraw中也有预测氢谱和碳谱化学位移的功能，画好键线式之后选中，菜单栏上点击Structure -> predict 1H-NMR shifts，即点即出，对于定性、半定量地与实验值对比，还是很方便的。
PS5：推荐一篇计算NMR化学位移的实例应用paper，DOI: 10.1021/acs.jctc.7b00380。
PS6：几个常见谱图数据库https://mason.gmu.edu/~sslayden/Lab/spec-db.htm
参考资料
1. http://gaussian.com/nmr/
2. Sobereva. 《谈谈如何又好又快地计算NMR化学位移》
http://bbs.keinsci.com/thread-4434-1-1.html
