查看: 8630|回复: 43
CASSCF计算双自由基以及双自由基特征的计算
在线时间2514 小时
帖子威望eV
CASSCF计算双自由基以及双自由基特征的计算
文/Sobereva(3)& &
双自由基就是指一个体系里有两个未成对儿电子，它们的自旋相反，因此是单重态的体系。
普通的单重态体系是闭壳层的，所有电子都是配对儿的。最理想的双自由基体系，是两个自旋彼此相反的单电子之间完全不配对儿，二者的空间分布彼此完全分离。而多数双自由基体系，情况介于二者之间，即两个单电子之间部分配对儿，分布没有完全分离。
本文主要讨论的是Yamanaka等人提出的双自由基特征y的计算。其数值从0到1，0对应闭壳层，越接近1说明体系双自由基特征越明显，1就是理想双自由基。本文还顺带讨论怎么用CASSCF计算双自由基体系。
文中Multiwfn用的是撰文时最新版3.3.7(dev)，计算用Gaussian09 B.01完成。
2 CI波函数下的双自由基特征
在笔者《小谈几种量子化学方法对氢分子解离过程的描述》（）一文中曾详细介绍了CI、RHF、UHF、VB方法对氢分子解离成两个氢原子过程的处理。文中提到，CI可以完美地描述这个过程。H2的CI波函数可以写为
ψCI=c1*ψG+c2*ψD
其中ψG是基态组态，即两个电子都占据H2的成键轨道（HOMO）的组态函数，ψD是两个电子都激发到H2的反键轨道（LUMO）的组态函数。ψCI也可以等价地用两个氢原子a、b的1s原子轨道来描述，记为Xa和Xb
ψCI=[(c1+c2)*(|Xaα Xaβ|+|Xbα Xbβ|)+(c1-c2)*(|Xaα Xbβ|-|Xaβ Xbα|)]/2
其中| |是slater行列式符号。诸如Xbα就代表Xb上占据α电子。
在H2平衡距离下，c1几乎为1，c2几乎为0。c1不完全为1是因为必须引入ψD才能体现电子动态相关。
在解离过程中，c1不断减小，c2不断增大。
当H2彻底解离为两个相距很远的氢原子后，成为完美的双自由基体系，两个氢原子各带一个电子且彼此自旋相反，空间上完全分离。此时c1=1/√2，c2=-1/√2，两个组态的权重（系数的平方）完全一致了。ψCI=(ψG-ψD)/√2这种波函数的形式不容易直接其物理含义，但是改写为用原子轨道表示，即(|Xaα Xbβ|-|Xaβ Xbα|)/√2，意义一下子就清楚了，这是正是两个氢原子各带的一个电子所构成的单重态自旋匹配波函数。
可见，双激发组态的权重c2^2越大，双自由基特征越强，和基态组态权重c1^2相同时就是完美双自由基。基于这个考虑，Yamanaka等人在Chem. Phys. Lett., 33, 330 (1975)中将y=2*c2^2定义为双自由基特征，完美双自由基时c^2=1/2，故y恰为1。
2.2 CASSCF计算双自由基C4H8实例
CI波函数描述双自由基在形式上是完美的，实际中用的更多的是CASSCF波函数，即不光变分CI系数，还变分轨道系数。描述双自由基至少用CAS(2,2)，把那两个单电子轨道以相位相同和相位相反方式组合构成的两个分子轨道纳入活性空间里，这样才能产生合理描述双自由基所需要的ψG和ψD。
解离后的H2可以说是最最最简单、最完美的双自由基的模型体系，实际的双自由基体系还有许许多多其它的轨道，但是我们通常还是能从中找出与构成双自由基直接相关的两个分子轨道并将之纳入活性空间的，很多情况下这两个轨道就是HOMO和LUMO。比如C4H8这个体系，分子两头各缺一个氢，故必然是单电子出现的位点，在HF/6-31G*下计算后，HOMO和LUMO分别是（等值面数值为0.1）
1.png (113.09 KB, 下载次数: 3)
16:35 上传
回想H2的成键轨道和反键轨道的形态，很明显，这里HOMO就是分子两端的单电子轨道以相位相同方式组合，而LUMO就是以相位相反方式组合产生的。因此将这两个轨道弄进CAS(2,2)活性空间就可以至少定性正确描述这个双自由基了。值得一提的是，LUMO轨道中在分子中间也出现了一坨，这和表现双自由基关系不直接，但是问题不大，反正之后CASSCF还会对轨道变分。
这里我们就完整地算一下C4H8，看看它的双自由基特征是多少。本文用的C4H8的结构：
C& &&&-0.& & 1.& & 0.
H& &&&-0.& & 2.& & 0.
H& &&&-0.& & 2.& &-0.
C& &&&-0.& & 0.& & 0.
H& &&&-1.& &-0.& & 0.
H& &&&-1.& &-0.& &-0.
C& && &0.& &-0.& & 0.
H& && &1.& & 0.& &-0.
H& && &1.& & 0.& & 0.
C& && &0.& &-1.& & 0.
H& && &0.& &-2.& &-0.
H& && &0.& &-2.& & 0.
整个过程是
第一步：#P HF/6-31G*
第二步：#P CAS(2,2)/6-31G* guess=read
我们做CASSCF时不需要调换轨道顺序，因为HOMO和LUMO就是要用的两个轨道。输出的信息中看到
& && && & Configuration& && && &1 Symmetry 1 10
& && && & Configuration& && && &2 Symmetry 4 ab
& && && & Configuration& && && &3 Symmetry 1 01
& &( 1)& &&&EIGENVALUE& & -156.
(& & 1) 0.7961444 (& & 3)-0.6051067 (& & 2) 0.0000000 (
其中组态1就是ψG，组态3就是ψD。c1=0.7961444，c2=-0.6051067，故此体系的双自由基特征=2*(-0.=0.732。
2.3 用UHF的自然轨道做CAS
也经常有人用非限制性计算(UHF/UDFT)得到的自然轨道作为CASSCF的活性空间的轨道算双自由基，这种轨道称为UNO (UHF natural orbital)。这里也用C4H8来演示一下。如果不熟悉怎么用非限制性计算来计算双自由基，务必先看看《谈谈片段组合波函数与自旋极化单重态》（）。
我们先看看UNO轨道形状是什么样。先用#P UHF/6-31G* guess=mix计算来得到对称破缺波函数，此时UHF波函数就存到了check文件里，然后再用#P HF/6-31G* guess=(naturalorbitals,read,save,only)算一次。这说明，从check里读取波函数(read)，转化成无自旋自然轨道(naturalorbitals)，保存到check文件里(save)，然后任务就停掉(only)。之后把chk文件转化成fch，用文本编辑器打开fch文件，把第一行内容改为saveNO，再用Multiwfn打开此fch文件并进入主功能0。在图形窗口上方选orbital info.-show all，我们看到UNO轨道信息
Orb:& & 13 Ene(a.u./eV):& &&&0.000000& && & 0.0000 Occ:&&1.996827 Type: A+B
Orb:& & 14 Ene(a.u./eV):& &&&0.000000& && & 0.0000 Occ:&&1.996534 Type: A+B
Orb:& & 15 Ene(a.u./eV):& &&&0.000000& && & 0.0000 Occ:&&1.996049 Type: A+B
Orb:& & 16 Ene(a.u./eV):& &&&0.000000& && & 0.0000 Occ:&&1.123073 Type: A+B
Orb:& & 17 Ene(a.u./eV):& &&&0.000000& && & 0.0000 Occ:&&0.876927 Type: A+B
Orb:& & 18 Ene(a.u./eV):& &&&0.000000& && & 0.0000 Occ:&&0.003951 Type: A+B
Orb:& & 19 Ene(a.u./eV):& &&&0.000000& && & 0.0000 Occ:&&0.003466 Type: A+B
（注：如果不把fch文件开头改成saveNO，那么占据数就会显示到能量ene那列去，如果只是为了看轨道倒是无妨，但如果还想用Multiwfn计算性质，如密度、静电势、ELF等，结果就错了）
从输出信息中看到，&=15号的UNO差不多都是双占据轨道，&=18号的几乎都是零占据的虚轨道。而RHF计算时原本是HOMO和LUMO的16、17号轨道现在占据数明显偏离整数，如下一节所示，这正是双自由基体系的明显特征。来看看它们的形状
2.png (119.44 KB, 下载次数: 2)
16:36 上传
对比HOMO和LUMO的图，会看到其实UNO 16和HOMO没明显差别，UNO 17和LUMO的主要差异在于分子中间的那一坨多余的区域没了（isovalue设小后还是会看到它的存在，只不过数值减小了很多），因此更纯粹地表现出分子两端的单电子轨道的组合。对于其它体系可能分子轨道和UNO的差异更大，如果届时从分子轨道中不好找出合适的用于CASSCF活性空间的轨道，就可以考虑用UNO。
接下来用#P CAS(2,2)/6-31G* guess=read做计算即可，这些UNO就被用作初猜轨道了，且UNO 16、17都处在活性空间。结果和上一节基于MO来做是一致的。
上面用UNO做CAS计算共经历了三步：
#P UHF/6-31G* guess=mix
#P HF/6-31G* guess=(naturalorbitals,read,save,only)
#P CAS(2,2)/6-31G* guess=read
实际上可以化简为两步
#P UHF/6-31G* guess=mix
#P CAS(2,2,UNO)/6-31G* guess=read
这代表，做CASSCF的时候，先从check文件中读取上一步的UHF波函数，自动转化为UNO，然后再用它作为CASSCF的初猜。
这里都是假定与HOMO、LUMO编号对应的那两个UNO是要被纳入活性空间的，但如果和实际情况不符，要取的是别的UNO，那么应该用guess(read,alter)来调换UNO轨道顺序，就像平时调换MO顺序做CASSCF那样。（无论是两步方法还是三步方法，都是这么用alter调换UNO顺序）
大家若有兴趣，可以在CASSCF做完之后观看一下chk里的轨道，会发现CASSCF优化出的16、17号轨道相较HOMO、LUMO，和UNO 16、17号在形状上明显更为接近，故此例使用UNO比用MO作为CASSCF初猜轨道可能更有优势。不过，如果基于MO做CASSCF就能正常收敛，如C4H8，其实没必要非得用UNO。
3 UHF波函数的双自由基特征
这里只讨论UHF的情况，UDFT的情况在处理上是完全一致的。
由于双自由基是单重态的，直接用UHF的结果和RHF是相同的，必须用guess=mix等关键词获得对称破缺解，这在《谈谈片段组合波函数与自旋极化单重态》（）已经谈了很多了，这里不再累述，只着重说说怎么计算UHF波函数的双自由基特征。
UHF波函数下，双自由基可以基于三种方法来计算：
(1)基于alpha和beta轨道重叠来计算
(2)基于&S^2&来计算
(3)基于UNO占据数来计算
对于H2这种理想体系，没有其它轨道影响，这三种形式是等价的，而实际情况由于其它轨道的影响，三种方式计算结果会有一定的差异。下面就分别谈谈这三种计算方式，并且用C4H8作为实例。
3.1 基于alpha和beta轨道重叠来计算
UHF波函数中alpha轨道彼此间是正交归一的，beta轨道彼此间也是正交归一的，但是alpha和beta轨道之间并非总是正交归一。对于双自由基体系，我们可以看成RHF波函数原本的一个双占据轨道（通常是HOMO）发生了分裂，形成空间分布并不重合的一个alpha轨道和一个beta轨道。二者空间分布分离越大，显然双自由基特征越强，这可以用它们重叠积分来衡量。
基于自旋投影UHF理论，Chem. Phys. Lett., 33, 330 (1975)中给出了UHF波函数双自由基特征表达式
y=1-(2*T)/(1+T^2)
其中T是空间发生分离的那一对儿alpha和beta轨道间的重叠积分。这个表达式和前面的y=2*c2^2形式是内在对应的。
用Multiwfn可以计算这个重叠积分。我们先用#P UHF/6-31G* guess=mix计算体系得到对称破缺波函数，并获得fch文件，载入进Multiwfn后进入主功能100，再进功能5，选选项2，我们就得到了每个alpha轨道和序号对应的beta轨道的重叠积分：
Overlap between the& & 13th alpha and beta orbitals:& & 0.986728
Overlap between the& & 14th alpha and beta orbitals:& &-0.961990
Overlap between the& & 15th alpha and beta orbitals:& &-0.995562
Overlap between the& & 16th alpha and beta orbitals:& & 0.042710
Overlap between the& & 17th alpha and beta orbitals:& & 0.149644
Overlap between the& & 18th alpha and beta orbitals:& & 0.939765
第16号alpha和beta轨道的重叠积分0.042710就是上式中的T。因为这两个轨道，如下图所示，正是分别描述那两个单电子的轨道：
（注：Multiwfn里beta序号在alpha序号之后。载入文件后屏幕提示从77号轨道开始是beta轨道，因此在主功能0看第16号beta轨道应该选77+16-1=92号轨道）
3.png (88.02 KB, 下载次数: 2)
16:37 上传
把T=0.042710代入y表达式，得到双自由基特征为0.915，比2.2节基于CAS(2,2)波函数得到的0.732大不少。
3.2 基于&S^2&来计算
波函数的自旋平方算符期望值&S^2&和精确值S^2_Exact=S(S+1)之间有如下关系（见Szabo的Modern quantum chemistry p107）
4.png (3.25 KB, 下载次数: 2)
16:36 上传
其中Nα和Nβ分别是alpha和beta电子数，S_ij是i号alpha和j号beta轨道重叠积分。由此式可见，如果每个alpha占据轨道都能和一个beta占据轨道唯一地匹配（积分值为1）而和其它的正交（积分值为0），那么&S^2&就是精确值S^2_Exact，RHF、ROHF波函数就属于这种情况。而对于UHF计算，只要由于自旋极化导致alpha和beta占据轨道分布发生分离，那么&S^2&肯定会偏离精确值，这就是常说的UHF波函数的自旋污染问题。
假设只考虑对应单电子的那一对儿alpha和beta轨道而忽略其它的，上式可简化为
ΔS^2=&S^2&-S^2_Exact=1-T^2
值得一提的是，双自由基体系S^2_Exact=0，因此对于理想双自由基，即T=0的情况，&S^2&应为1。
将T=√(1-ΔS^2)代入3.1节y的表达式，就可以得到基于&S^2&与精确值偏差量ΔS^2的双自由基特征表达式
y=1-2*√(1-ΔS^2)/(2-ΔS^2)
不过这个式子往往没法用。比如C4H8这个体系在#P UHF/6-31G* guess=mix计算后输出的&S^2&=ΔS^2=1.0142，代入y表达式的话，根号里就成负值了。这是由于其它alpha和beta轨道间也出现了一定的分离，明显影响了体系的&S^2&所致，只考虑一对儿alpha和beta轨道不合适。
3.3 基于UNO占据数来计算
重叠积分T和UNO占据数也有直接联系。设只考虑对应单电子的那一对儿alpha和beta轨道而忽略其它的，则编号对应HOMO和LUMO的UNO的占据数（记为n_HOMO和n_LUMO）将分别为1+T和1-T，故T=(n_HOMO-n_LUMO)/2。因此，我们可以通过UNO占据数获得T再获得y。对于理想双自由基，n_HOMO和n_LUMO都将为1。
2.3节中我们已经看到了C4H8的UNO占据数了，n_HOMO=1.123073，n_LUMO=0.876927，故T=0.123，y=0.758。和3.1节直接通过重叠积分算出来的0.915有不小差别，这是因为这里得到自然轨道占据数的时候把其它轨道的影响也纳入进去了，结果倒是比较接近于CAS(2,2)得到的0.732。因此，如果不想做CAS，而想基于方便快速的UHF波函数来定量考察双自由基特征的话，通过UNO占据数计算是比较理想的选择。C4H8的0.7几的双自由基特征说明此体系是较明显的双自由基，尽管alpha和beta单电子分布还没有完全分离，离理想双自由基还有点距离。
实际上，只需要UNO占据数的话，并不需要2.3节那么多关键词而且分两步，只需要用#P UHF/6-31G* guess=mix pop=NO算一次就够了，会看到这样的输出
& && && && && && && && && && && && & 16& && && &&&17& && && &&&18& && && &&&19& && && &&&20
& &&&Eigenvalues --& &&&1.12307& &0.87693& &0.00395& &0.00347& &0.00317
& &1 1& &C&&1S& && && && &&&0.01911& &0.00845&&-0.07955& &0.10559&&-0.08133
& &2& && &&&2S& && && && &&&-0.04188&&-0.01951& &0.15207&&-0.21148& &0.16843
Eigenvalues便是这些UNO的占据数，即UHF无自旋密度矩阵的本征值。
相对于用CASSCF计算双自由基，UHF虽然不如它准确，但好处是快、方便，而且有自旋密度，通过自旋密度能直接看到未成对儿alpha和beta电子分布，只要在Multiwfn里进入主功能5后选自旋密度然后按提示操作就能得到自旋密度等值面图，单电子在什么位置一目了然。而CI波函数虽然精确、完美表现了两个单电子处于自旋相反的状态，但是却得不到自旋密度。
如3.2节的&S^2&公式可见，只要用UHF描述双自由基，就一定会出现自旋污染，而且对于理想双自由基偏差值高达1.0，于是就总有人盲目地批判UHF/UDFT有自旋污染，好像偏差值这么大就根本不能用似的，这是大错特错！&S^2&与理想值偏差大，虽然意味着波函数与真实波函数偏差大，但并不意味着能量不合理！对于H2完全解离这样的理想双自由基，UHF和FCI结果是一致的，尽管前者表面上看自旋污染达到1.0已经大得要命了。自旋污染对于DFT波函数更是难以讨论，本来DFT引入波函数的初衷就是为了获得可靠的动能密度而不是构造真实的波函数，故不应当根据波函数的合理性来讨论DFT结果的合理性。从实效上来讲，UDFT计算双自由基还是不错的。
北京科音自然科学研究中心：&&不定期开办各层次量子化学、分子动力学、Multiwfn程序培训
思想家公社的门口Blog：
Multiwfn量子化学波函数分析程序主页：
计算化学公社论坛：
思想家公社QQ群1号：号：。用于讨论理论、计算化学，两个群讨论范畴相同，可加入任意其一但不可都加入，申请信息必须注明具体研究方向，否则一概不批。
此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!
在线时间67 小时
帖子威望eV
Level 3 能力者, 积分 308, 距离下一级还需 192 积分
谢谢Sob这么努力地分享经验！！
在线时间3 小时
帖子威望eV
Level 2 能力者, 积分 60, 距离下一级还需 90 积分
谢谢，收获很大。
在线时间217 小时
帖子威望eV
Level 5 (御坂), 积分 1681, 距离下一级还需 2319 积分
做过不少双自由基体系的计算，发现DFT方法很有效，结果与实验也吻合地很好。一直想学习如何使用CASSCF呢
PS: Masayoshi Nakano也做了不少有关y值的工作
在线时间36 小时
帖子威望eV
Level 3 能力者, 积分 283, 距离下一级还需 217 积分
我之前看到有人用qcisd(t)/ccsd(t)中的T1检测来判断是不是要考虑multi-reference的效应，感觉这个T1和文中T好像很相似，他们是一样的么？或者具有正相关的关系？
在线时间2514 小时
帖子威望eV
我之前看到有人用qcisd(t)/ccsd(t)中的T1检测来判断是不是要考虑multi-reference的效应，感觉这个T1和文中T ...
虽然有一点关系，但形式还是很不一样。T1诊断类似于判断CI中的单激发组态所占比例，其数值越大，体现出HF基态波函数的主导性越差，越不适合作为参考态。
北京科音自然科学研究中心：&&不定期开办各层次量子化学、分子动力学、Multiwfn程序培训
思想家公社的门口Blog：
Multiwfn量子化学波函数分析程序主页：
计算化学公社论坛：
思想家公社QQ群1号：号：。用于讨论理论、计算化学，两个群讨论范畴相同，可加入任意其一但不可都加入，申请信息必须注明具体研究方向，否则一概不批。
此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!
在线时间36 小时
帖子威望eV
Level 3 能力者, 积分 283, 距离下一级还需 217 积分
虽然有一点关系，但形式还是很不一样。T1诊断类似于判断CI中的单激发组态所占比例，其数值越大，体现出HF ...
好的，谢谢sob！
在线时间2 小时
帖子威望eV
Level 1 能力者, 积分 14, 距离下一级还需 11 积分
T1诊断值如果不是很高，0.02左右的需要考虑多参考效应么
在线时间2514 小时
帖子威望eV
T1诊断值如果不是很高，0.02左右的需要考虑多参考效应么
属于模棱两可。对于DFT，T1越高，建议用HF成份越低的泛函，如果明显大于0.02则应当用纯泛函。对于CCSD(T)，0.02左右一般也能得到合理结果。
北京科音自然科学研究中心：&&不定期开办各层次量子化学、分子动力学、Multiwfn程序培训
思想家公社的门口Blog：
Multiwfn量子化学波函数分析程序主页：
计算化学公社论坛：
思想家公社QQ群1号：号：。用于讨论理论、计算化学，两个群讨论范畴相同，可加入任意其一但不可都加入，申请信息必须注明具体研究方向，否则一概不批。
此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!
在线时间452 小时
帖子威望eV
Level 4 (黑子), 积分 1349, 距离下一级还需 151 积分
sob老师，现在有没有软件能做 多参考-DFT呢（理论上 密度泛函理论用多参考好像有点奇怪。。。）？
在线时间2514 小时
帖子威望eV
sob老师，现在有没有软件能做 多参考-DFT呢（理论上 密度泛函理论用多参考好像有点奇怪。。。）？
但Grimme有个DFT/MRCI方法（JCP,111,5645），也有相应的程序，但可能不是你说的那个意思
北京科音自然科学研究中心：&&不定期开办各层次量子化学、分子动力学、Multiwfn程序培训
思想家公社的门口Blog：
Multiwfn量子化学波函数分析程序主页：
计算化学公社论坛：
思想家公社QQ群1号：号：。用于讨论理论、计算化学，两个群讨论范畴相同，可加入任意其一但不可都加入，申请信息必须注明具体研究方向，否则一概不批。
此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!
在线时间192 小时
帖子威望eV
Level 4 (黑子), 积分 802, 距离下一级还需 698 积分
请教sob老师一下
Yamaguchi 这个计算耦合常数的方程& &J= (EHS-EBS)/(S^2)HS-(S^2)BS 中的S^2，是用S^2理想值还是用计算结果中输出S^2值？
在线时间2514 小时
帖子威望eV
请教sob老师一下
Yamaguchi 这个计算耦合常数的方程& &J= (EHS-EBS)/(S^2)HS-(S^2)BS 中的S^2，是用S^2 ...
计算结果中的
北京科音自然科学研究中心：&&不定期开办各层次量子化学、分子动力学、Multiwfn程序培训
思想家公社的门口Blog：
Multiwfn量子化学波函数分析程序主页：
计算化学公社论坛：
思想家公社QQ群1号：号：。用于讨论理论、计算化学，两个群讨论范畴相同，可加入任意其一但不可都加入，申请信息必须注明具体研究方向，否则一概不批。
此账号为诸Sobereva共用
Money and papers are rubbish, get a real life!
在线时间452 小时
帖子威望eV
Level 4 (黑子), 积分 1349, 距离下一级还需 151 积分
但Grimme有个DFT/MRCI方法（JCP,111,5645），也有相应的程序http://bbs.keinsci.com/thread-20 ...
在线时间11 小时
帖子威望eV
Level 2 能力者, 积分 65, 距离下一级还需 85 积分
sob老师，麻烦问一下，我用您给的C4H8的坐标就可以算出同样的结果，为什么用我自己建模得到的坐标跟您的值差好多啊}