这是本文档旧的修订版！
计算动态超极化率、双光子吸收截面TPA Cross Section
本文参考文献资料：
这里我们介绍双光子吸收截面σTPA的计算过程。根据文献，
- 精细结构常数α=e$^2$/4πε$_0$cħ，其中e电子电荷，ε$_0$真空介电常数，c真空光速，ħ约化普朗克常数。
- 整数N取决于实验设置，很多文章中都取N=4
- γ约化形式的二阶超极化率γ(-ω,ω,ω,-ω)，在ADF软件中，超极化率选项选择γ TPA，并设置ω1，即可计算得到γ(-ω1,ω1,ω1,-ω1)
- ω频率
- Im[]表示取虚部
- σ$_{TPA}$单位为Göppert-Mayer (1 GM = 10$^{−50}$ cm$^4$ s photon$^{−1}$)
原子单位下:
- α=1/137.035999679 a.u.
- e=1 a.u.
- ħ=1 a.u.
- ADF得到的γijkl值单位为a.u.
如果我们得到原子单位的γTPA、ω，那么GM单位的值 = (4*3.1415927$^3$*(0.529177$^4$)*2.418884*10/(137$^2$*15))*γTPA*ω$^2$ = 0.0008343253*γTPA*ω$^2$
我们仅仅选取其中1-NMe2为例计算TPA截面。
计算背景
分子结构
文章附录材料提供了几何结构优化完成之后的分子结构（实际上使用ADF也同样可以优化，参考：geoopt）
1-NMe2分子结构：（点击后，复制文字，在Input窗口粘贴即可）
参数设置
关于泛函：
SAOP是ADF独有的一种模型势，响应性质计算方面，一般优于其他LDA、GGA、metaGGA泛函，因此本文采用SAOP。
关于基组：
原则上而言，基组越大越精确，但是根据文献中介绍，基组对结果影响不大，因此本文采用了很小的DZP基组。如果体系中包含较重的元素，还是建议为重元素设置TZP，甚至TZ2P、QZ4P基组。设置方法参考：diffbasisfordiffelement。SAOP不允许使用Frozen Core，而且TPA计算也不建议使用冻芯近似，因此设置为None。
数值精度：
这里为了快速计算得到结果，因此使用了默认值Normal，实际上建议选择Good。每提高一个等级，计算量都会数倍增长。
TPA计算，选择γ TPA，并设置ω1，后面的单位，用户可以点击后自由选择，例如这里选择了eV，并设置ω1=1.56 eV；lifetime这个值最好通过拟合分子的吸收数据来获得，不过这个值在相似的分子之间变化不大，因此估计值并不困难。文献“Journal of Chemical Physics 123, 174110 , 2005”中使用了0.004 Hatree，本文使用0.0034 Hatree。
γ TPA选项，实际上使用的是基于TDDFT的damped cubic response theory来计算。
保存任务并提交作业。
结果查看
γTPA的结果，主要在*.out文件中。搜索“IMAGINARY SECOND HYPERPOLARIZABILITY”即看到γ的虚部（更前面则是实部数据）：
----------------------------------------------------
IMAGINARY SECOND HYPERPOLARIZABILITY
----------------------------------------------------
X Y Z
X X X 107766780.3550 3358308.0941 401555.9119
Y 3562512.2519 29233046.0292 1337070.5306
Z 387007.5189 1345541.2048 81403.7813
X Y X 3562512.2519 29233046.0292 1337070.5306
Y 42404008.8720 -294151.9274 1021197.6465
Z 1470347.0140 610561.4349 49380.5430
X Z X 387007.5189 1345541.2048 81403.7813
Y 1470347.0140 610561.4349 49380.5430
Z -472299.5636 53831.2906 -2282.2650
Y X X 3358315.1989 42850778.7646 1474115.9282
Y 29233049.0338 -87861.6272 615867.8827
Z 1345561.1655 1039131.8245 49700.3983
Y Y X 29233049.0338 -87861.6272 615867.8827
Y -294150.6247 100180120.4322 2161451.0943
Z 610581.8039 2177757.0074 70576.2211
Y Z X 1345561.1655 1039131.8245 49700.3983
Y 610581.8039 2177757.0074 70576.2211
Z 53832.5230 -520809.3156 -11313.5051
Z X X 401555.7471 1474104.7497 -479494.6257
Y 1337078.5881 615854.4081 54685.0684
Z 81451.5589 49699.1777 -2183.0447
Z Y X 1337078.5881 615854.4081 54685.0684
Y 1021203.3117 2161444.4591 -527540.2199
Z 49380.1514 70625.0554 -11490.4655
Z Z X 81451.5589 49699.1777 -2183.0447
Y 49380.1514 70625.0554 -11490.4655
Z -2281.5329 -11310.6297 4031.7741
根据公式，我们需要虚部中γijkl，ijkl=xxyy、xyyx、xyxy；yyxx、yxxy、yxyx；yyzz、yzzy、yzyz；zzyy、zyyz、zyzy；zzxx、zxxz、zxzx；xxzz、xzzx、xzxz共18个数据，以及α=β时，各出现3次的xxxx、yyyy、zzzz。γTPA，这27项之和（用户可以借助Excel导入文本形成数据表格，借助其公式功能较为高效、不易出错地整理完毕所有计算结果）。
上面数据中，Im[γXXXX] = 107766780.3550 a.u.，Im[γXXYY] = 29233046.0292 a.u.，读者可以类似找到其他数据。
本例中，27个数据之和$γ^{TPA}$=824647744.9551 a.u.，ω = 1.56 eV = 1.56/27.2113956 = 0.0573289 a.u.（在*.run文件中也可以找到这个数值），带入前面的公式中，得到2261 GM。
分别计算1.46~1.76 eV之间的结果（计算输入输出文件下载），得到$σ^{TPA}$曲线：
峰的位置与文献完全符合，但是峰值差一倍（实际上我们这里的计算得到的峰的位置、强度，比文献中的结果更接近实验结果）。γ值非常敏感，积分精度（Good与Normal之间大约对γ带来10%左右的差异）、键长的极其细微变化，都会导致很大的扰动，因此可以认为与文献结果是一致的（下图蓝色曲线）：
实验结果1-NMe2峰的位置在1.72eV，我们的计算结果峰位置在1.58eV附近，文献中为1.55eV；峰值实验结果360 GM，我们的计算结果是2376 GM，文献中为4246 GM。
我们类似计算了1-NH2，定性结果与文献一致（上图红色曲线）。其中实验峰值130 GM，峰的位置1.77 eV，我们计算峰值1097 GM，峰位置1.76 eV（只计算了4个点），文献计算峰值1071 GM，峰的位置1.77 eV。
致谢
感谢该文献作者Abdou Boucekkine教授，帮助确定公式中Im[γXXXX]、Im[γYYYY]、Im[γZZZZ]需要乘以3。