使用Multiwfn对第一超极化率做双能级和三能级模型分析
使用Multiwfn对第一超极化率做双能级和三能级模型分析
Using Multiwfn to perform two-level and three-level model analyses for first hyperpolarizability
文/Sobereva@北京科音
First release: 2019-Sep-6 Last update: 2023-Aug-11
1 前言
研究第一超极化率(beta)的文章里经常做一种叫双能级模型的分析用于阐明影响beta的关键因素、便于讨论不同体系beta大小差异的原因，比如在笔者的J. Comput. Chem., 38, 1574 (2017)中就用了这种方式讨论。之前笔者在《谈谈计算第一超极化率的双能级公式》（http://sobereva.com/361）里专们对此方法进行了介绍。本质上，双能级公式就是从《使用Multiwfn基于完全态求和(SOS)方法计算极化率和超极化率》（http://sobereva.com/232）文中介绍的计算beta的完整的SOS公式中简化来的。
使用双能级公式的前提之一是必须有一个激发态是关键态(crucial state)，也就是它对beta的贡献远大于其它激发态，因此可以只用基态和这个激发态来讨论beta。然而很多情况下关键态难以确定、模棱两可，甚至有时候明显能看出就是需要同时考虑两个激发态才能定性正确描述实际的beta，比如第一激发态和第二激发态振子强度都较大而它们的能量差很小，此时就明显不能强行用双能级分析了，得同时把这两个激发态一起考虑，笔者管这叫三能级模型，推导过程见下一节。
为了便于大家实现双能级和三能级分析，从2019-Sep-6更新的Multiwfn开始，在完全态求和(SOS)模块里新增了子功能20，专门进行双/三能及分析，使用相当容易，下面就进行介绍。首先笔者把方法的原理介绍一下，然后给出一个具体分析例子。Multiwfn可以在其主页http://sobereva.com/multiwfn免费下载。不了解Multiwfn的话看《Multiwfn FAQ》（http://sobereva.com/452）。
2 原理
Multiwfn手册的3.27.2节里对双能级、三能级公式的推导做了非常详细的说明，十分建议一看。这里只是简单提一下关键点。
SOS方法计算beta张量的ABC分量公式如下
对于静态极限（ω=0）的情况，ZZZ分量的公式可以写为以下形式
其中Δi是第i激发态的激发能，μ0i_Z是基态到第i激发态的跃迁偶极矩的Z分量，μij_Z是第i激发态到第j激发态的跃迁偶极矩的Z分量，Δμi_Z是第i激发态的（非弛豫密度对应的）偶极矩矢量与基态偶极矩矢量之差的Z分量。
由上式可见，静态beta的ZZZ分量可以写为每个激发态独自的贡献以及不同激发态之间耦合所产生的贡献之和。如果只有一个激发态起到主导地位，那就可以用下面的双能级模型讨论
如果有两个激发态同时起到关键地位，哪个都不能简单地忽略，那就得用下面的三能级模型讨论，包括两个激发态自身的贡献和二者的耦合项
再强调一下，用双/三能级模型讨论影响beta大小的因素的时候必须满足以下两条：
(1)beta的总值必须能被beta_XXX、beta_YYY、beta_ZZZ其一所主导，这样我们才能针对起到主导作用的分量来进行双/三能级模型的分析来解释影响beta总大小的主要因素（注：为此，有的时候可能需要在计算前对体系做恰当的旋转来使得以上三个分量之一最大化，从而便于双/三能级分析）
(2)打算用双能级分析时，必须有一个激发态对被考察的beta分量的贡献远大于其它激发态；打算用三能级分析时，必须有两个激发态的贡献远大于其它激发态。
通常，Donor-pi-acceptor型体系可以满足以上要求。分析前应当令Donor-acceptor连线的方向（或者更广义地说，电子激发时电荷转移的方向）平行于某个笛卡尔轴，以满足上面第(1)条。
3 实例
下面就用一个典型的Donor-pi-acceptor型体系“氨基-联苯-硝基”作为例子说一下怎么用Multiwfn做双能级和三能级分析。用到的输入文件满足Multiwfn手册3.21节开头里说的对输入文件的要求即可。简单来说，对于Gaussian用户，一般需要做TDDFT计算，并且写上IOp(9/40=4)关键词，计算后产生的fch和输出文件都需要保留。本例用到的文件D-pi-A.fchk和D-pi-A.out已经提供在了程序包的examples\excit目录下（注：这个计算任务是在CAM-B3LYP/6-31g(d)级别下完成的，对于考察超极化率来说，用这种不带弥散函数的基组不可能得到准确的结果，但本文仅作为演示，所以暂且忽略这点）。
首先我们产生Multiwfn的SOS模块的输入文件，这需要依赖D-pi-A.fchk和D-pi-A.out里的信息计算出各个态的偶极矩和态之间的跃迁偶极矩。启动Multiwfn后依次输入
examples\excit\D-pi-A.fchk
18 //电子激发分析
5 //计算各个态的偶极矩和态之间的跃迁偶极矩
examples\excit\D-pi-A.out
3 //将结果导出为当前目录下的SOS.txt文件
然后就可以做双/三能级分析了。在此之前，返回主菜单，进入主功能0看一眼体系的朝向，如下所示。可见体系是顺着X方向的，因此我们之后要考察的是beta_XXX。而且如果你有计算beta的常识的话，也会知道这种体系beta_XXX主宰了总beta，即考察beta_XXX就等于近似考察了总beta。
启动Multiwfn，输入
SOS.txt //载入当前目录下的SOS.txt
24 //(超)极化率分析
2 //通过完全态求和(SOS)方式研究（超）极化率
20 //双、三能级分析
此时程序将激发态的基本信息输出了出来，便于用户选择被考察的态
Excitation energy (E,eV) and transition dipole moment (X,Y,Z,total) between gro
und state to excited states (a.u.)
State 1 E= 3.9069 X= 0.44441 Y= -0.00044 Z= -0.00182 Tot= 0.44442
State 2 E= 4.0624 X= 2.52778 Y= -0.00309 Z= -0.00733 Tot= 2.52779
State 3 E= 4.4166 X= -0.00263 Y= -0.00167 Z= -0.02343 Tot= 0.02364
State 4 E= 4.7912 X= 0.00069 Y= 0.34055 Z= -0.01274 Tot= 0.34079
State 5 E= 4.8872 X= -0.00391 Y= 0.19205 Z= -0.17054 Tot= 0.25687
由上可见，第2激发态的跃迁偶极矩的模（Tot值）明显大于其它的态（相应地振子强度也大），而且它的激发能又比较低，因此这就是典型的关键态(crucial state)。接着输入
1 //考察X方向，即被分析的是beta_XXX
2 //将第2激发态用于双能级分析
马上输出了以下信息
Excited state 2
Excitation energy 4.0624 eV
Transition dipole moment component: 2.5278 a.u.
Oscillator strength component: 0.635943 a.u.
Variation of dipole moment component: 6.5629 a.u.
beta evaluated by the two-level model: 11289.099582 a.u.
前面介绍双能级公式中涉及的激发能、跃迁偶极矩、偶极矩相对于基态的变化在上面都给出了，对于后两者指的是我们选择的X分量。另外，由这些数据折算的振子强度的X分量也给出了（具体公式看手册3.27.2节）。利用这些信息，我们就可以讨论这个体系与与它类似的体系的beta的差异主要来自于哪些因素。
PS：大家也可以尝试把被分析的方向设为Y或Z，会发现算出的数值都非常接近于零，这也主要在于第2激发态的跃迁偶极矩的Y和Z分量都近乎为0。
接下来我们用三能级模型进行考察。在Multiwfn窗口里输入
20 //再次做双/三能级分析
1 //再次考察beta_XXX
1,2 //将激发态1、2用于三能级模型分析
结果如下
Excited state 1
Excitation energy 3.9069 eV
Transition dipole moment component: 0.4444 a.u.
Oscillator strength component: 0.018904 a.u.
Variation of dipole moment component: -1.0098 a.u.
Excited state 2
Excitation energy 4.0624 eV
Transition dipole moment component: 2.5278 a.u.
Oscillator strength component: 0.635943 a.u.
Variation of dipole moment component: 6.5629 a.u.
Transition dipole moment between states 1 to 2: 1.475397 a.u.
Contribution of excited state 1: -58.048845 a.u.
Contribution of excited state 2: 11289.099582 a.u.
Contribution of cross term: 927.898744 a.u.
beta evaluated by the three-level model: 11251.731412 a.u.
上面给出的第2激发态的信息和之前我们做双能级分析的时候得到的相同。上面的信息包含了三能级公式里面涉及到的所有量。当必须用三能级模型才能说明问题的情况，大家就可以对不同体系考察上述这些量对beta_XXX的影响。
从上面的输出也可以看到第1激发态自身的贡献（-58.0），以及第1激发态与第2激发态之间耦合产生的贡献（927.9）相对于总值11251.7来说都很小，体现出用双能级模型其实就足够了，把第1激发态纳入作为三能级模型考察也不会得到更多信息。也很容易理解为什么第1激发态的独自贡献（-58.0）那么小，因为双能级模型指出某个激发态对beta的贡献正比于其偶极矩相对于基态的变化量以及跃迁偶极矩的平方。
原理上来说，还可以做四能级、五能级...分析，但是由于基本不可能用得到，所以Multiwfn也就没留出相应功能。由于三能级分析的时候已经给出两个激发态自身以及激发态之间的耦合项了，因此若想同时考虑更多能级，靠多次使用三能级公式就可以等效实现（比如想用2、3、4激发态做四能级分析，通过做2-3、2-4、3-4三次三能级分析就可以等效实现）。
顺带一提，有时候我们的电子激发计算的输出文件里可能有很多态，比如50个乃至上百个，但使用双/三能级公式考察的时候一般只会牵扯到最低的几个，显然此时对输出文件里所有激发态都做偶极矩、跃迁偶极矩的计算完全没必要，会造成大量无意义的耗时。为此，可以恰当设定Multiwfn的settings.ini里的maxloadexc，比如设为5，那么主功能18的子功能5就只会读入前5个激发态，计算偶极矩和跃迁偶极矩也只对它们计算，比起算所有激发态会省时甚巨，此时产生的SOS.txt就足矣用于做双/三能级分析了。