计算化学公社

标题: SOS方法计算NLO性质，优选基组，以及对CPDFT的一些怨念 [打印本页]

作者
Author: yjcmwgk 时间: 2015-9-1 15:40
标题: SOS方法计算NLO性质，优选基组，以及对CPDFT的一些怨念
本帖最后由 yjcmwgk 于 2015-9-5 10:03 编辑

SOS方法计算NLO性质，优选基组，以及对CPDFT的一些怨念(Update 1)

最近我用CAM-B3LYP泛函，对NH3分子的NLO响应系数，借助Multiwfn，用SOS方法进行了一些测试，主要目的是想测试出该项计算中，极化函数该用多少。
在sob对原帖进行了评价后，我重新开贴吧，这次叫update 1。这update 1不但有原来的内容，也有sob提醒的内容，也有我新加入的内容。
此次测试用5种基组，Ahlrichs基组、带有aug弥散函数的Dunning基组、带有d-aug弥散函数的Dunning基组、Pople基组、sadlej基组。
介绍完毕，开始正文。

先列出测试结果：

(, 下载次数 Times of downloads: 69)

(, 下载次数 Times of downloads: 80)
注1：在用aug-cc-pv5z、daug-cc-pvqz、daug-cc-pv5z基组进行计算时，对氨分子进行TD计算，总激发态数目绝不止1000个。但是考虑到计算量，我只截取了前一千个激发态。
注2：def2-tzvpd、def2-tzvppd、sadlej-pvtz和sadlej+基组来自Basis Exchange，其他基组是高斯自带的。
注3：该测试结果仅对氨分子有效。我觉得以后我会去测试几个简单的共轭体系，苯甲酸（偶极矩体系）和三聚氰胺（八极矩体系）是两个不错的测试体系。
注4：动态超极化率和非线性光学响应系数的计算，用lambda=1064nm的原因是，这是最常用的激光器之一。

我的目的是：在泛函不变的情况下，增加极化函数的数量，看看NLO计算到底需要多少极化函数？增加到什么程度，NLO计算结果就稳定下来，不会大变了？
我必须说，我对Pople基组的测试结果是绝望的、对Dunning基组是失望的、对Ahlrichs基组是满意的、对sadlej基组是欣喜的。
看表格吧。
（1）对于Pople基组，增加极化函数的数量，无论是alpha，beta还是gamma，我增加极化函数直到(3df,3pd)，也仍旧没有稳定到某个数值的迹象。Pople基组没有d及以上角动量弥散函数，在超极化率计算上式完败的。
（2）对于Ahlrichs基组，尽管def2-tzvppd比def2-tzvpd多了近百个激发态，但是无论是alpha，beta还是gamma都没有大的变化。而def2-tzvppd比def2-tzvpd的gauss函数数量一般来说要大上15-25%呢。所以我到现在仍旧觉得def2-tzvpd够用了啊。此处@ChemiAndy @jiewei @sobereva @14xlliu @卡开发发 @我本是个娃娃  @小范范1989 @superrice
（3）对于Dunning基组来说，单说“增加极化函数数量，使得结果稳定到某个数值”这唯一目的，对于alpha的稳定来说，是aug-cc-pvdz就达到了稳定。但是beta值的稳定，要增加到aug-cc-pvqz。而对于gamma来说，似乎稳值遥遥无期！况且还要考虑到Dunning基组远大于上述两个基组的gauss函数数量。sob说由于超极化率对弥散函数很敏感，也可以同时测试下每个角动量都含两层弥散的d-aug-cc-pVnZ系列，测试发现对于alpha来说daug-cc-pVnZ系列并不比aug-cc-pVnZ系列更好，两者收敛至同一个稳定值；对于beta来说d-aug-cc-pVnZ系列意义也不大，对于gamma来说，d-aug-cc-pV5Z要比aug-cc-pV5Z的数值偏离了大概十分之一。同时，对于gamma来说，aug-cc-pV5Z与aug-cc-pVQZ差距明显，但d-aug-cc-pV5Z与d-aug-cc-pVQZ差距已经很小，似乎已经达到稳定。这说明每个角动量都含两层弥散的d-aug-cc-pVnZ系列在gamma值的计算上明显好于aug-cc-pVnZ系列。然而这并没有什么用。仅仅是一个氨分子，采用d-aug-cc-pvqz和d-aug-cc-pv5z基组就分别需要近400个和近700个Gauss Func了！这在实际研究中，体系稍微变大，就无法计算了。
（4）sadlej基组的测试结果让人欣喜！作为专门为超计划率计算而设计的基组，sadlej+基组仅用130+个Gauss Func就达到了表中aug-cc-pv5z的计算水平。与aug-cc-pvtz基组相比，sadlej+基组大大简化了对s和p函数的描述，细化了对d型函数的描述，但是删除了f型函数的描述。这是否表明，在SOS(TDDFT)的NLO计算中，aug-cc-pVnZ中对f极化函数的描述事实上没有太大必要？同时我们也看到，对gamma的计算，sadlej+基组仍不及aug-cc-pvtz基组，是否说明我们可以通过在sadlej+基组中通过某些替换来实现双弥散函数的描述？当然这又不是特别必要，毕竟在实际研究中很少遇到氨分子这么简单的小分子。sadlej+基组的Gauss Func比def2-tzvpd要多一些，但是准确度却明显要好。这是个取舍问题。
（5）当我进行过d-aug-cc-pv5z基组的测试后，我的直观感受是：与其追求双弥散函数来提高这十分之一的精度，不如关注计算量是不是电脑能承受的，以及好几千个Gauss Func在一起进行SCF迭代能不能收敛的问题。此处也@ChemiAndy @jiewei @sobereva @14xlliu @卡开发发 @我本是个娃娃  @小范范1989 @superrice
以上是我的一个粗糙测试。还请各位keinsciers批评指正

----------以上是基于测试结果给出的一些东西，下面是一些前期工作的总（怨）结（念）。------------
SOS大法好！可惜对基组要求实在太高！以前我折腾CPDFT的时候，Cam-B3LYP或者LC-BLYP泛函，外加6-31+G(d)包打天下的。
然而CPDFT那是撞大运、“霰弹枪式”模拟。当年俺们优选了N的N次方个泛函也没找到较为通用的泛函(Cam-B3LYP和LC-BLYP勉强)。后来弄那个阻尼震荡，到死也找不到通用的阻尼系数计算方式来，连估算方式都欠奉。我们组做NLO模拟（CP-DFT）的人很少，基本上就是我师兄(Y. Zhang)，我爱人(L. Zhang)和我自己。师兄出国后也在NLO做了一些工作，现在不做了。我有段时间没做了，现在又拾起来了。我爱人不做了。背景介绍完毕。
我们遇到了4个问题
第一，以J. Phys. Chem. C 2013, 117, 1833的Fig. 7和J. Phys. Chem. C 2012, 116, 21973的Fig. 3为例，beta的模拟，莫说数值，连趋势都做不好！superrice也说“关于beta和gamma，一般来说确实不容易模拟。碰到有的体系你会发现连正负号都是错的，根本不要说趋势了。有时候要撞上大运了，你会发现DFT居然比CCSD的结果还要好。你可以试试functional tuning的办法来算beta：10.1002/cphc.201300256  还有 10.1021/jp507226v  运气好的话能跟接近mp2的水平。”看来他的怨念似乎也挺深的。事实上superrice也多次提到所谓的functional tuning，其实还是挑泛函的代名词，仍然是“霰弹枪式”模拟。
第二，以J. Phys. Chem. A 2012, 116, 10249的Fig. 5为例，Harmonic light intensity as a function of the polarization angle Ψ by polar representation的模拟，与实验值差距较大，只能给个趋势。当然，能给趋势已经不错了，superrice说“碰到有的体系你会发现连正负号都是错的，根本不要说趋势了”，这事儿我们也遇到过多次。看起来，倒是最为粗糙的FF方法反而安全？？
第三，关于NLO分子材料与外加电磁场发生共振震荡。如果不添加阻尼系数，则会产生过度震荡的问题，如J. Phys. Chem. A 2012, 116, 10249的Fig. 2(b)。如果添加阻尼系数，则无法找到通用的阻尼系数计算方式来，连估算方式都欠奉。我们在ChemPhysChem 2015, 16, 1889的Eq. 6中的阻尼系数，事实上是通过拟合实验值得到的。如何通过理论方式比较准确地获得该值，尚未找到可靠方法。
第四，在ChemPhysChem 2015, 16, 1889的Fig. 6-7以及Turk. J. Chem. 2014, 38, 1046的Fig. 6（这是一篇会议论文，无论是否付费都可以点击直接下载）中，我们计算了Periodic vibration of electron density driven by an external electromagnetic field，动图在此。但这是用FF方法计算得到的。我们对自己很不满。请问有什么好方法么？况且我通过外加电场以FF方法来模拟电子云密度的转移，是在用静态模拟动态过程，实质上还是多个静态过程。这种过程最让我自己诟病的是，这办法等于默认了一个前提，即“无论外加电场如何变化，分子内部的电子云密度总能在一瞬间弛豫完毕，根本不需要弛豫时间”。这并不总是合理的！比如，对于1064nm的激发光源来说，分子外的外加电磁场方向的转换周期是3.5乘10的负15次方秒，外加电磁场维持同一方向的时间不超过1.8乘10的负15次方秒，再考虑到正弦函数在相当多的时候强度接近零，能明显对分子内电子有单一方向推力的时间估算为1乘10的负15次方秒算是合适的，而电子运动的速度不超过10的7次方米每秒，外层电子的运动速度甚至在10的6次方米每秒数量级，在10的负15次方秒内最多只能走一个纳米的距离。我们再考虑单个电子的运动速度必然远快于电子云整体的运动速度，电子云整体能移动半个纳米都算他移动得快了（想想气体理论中的平均自由程以及各种衍生概念，还有，每一个原子核都是电子的陷阱，电子哪儿能真的自由移动呢？此处还可以想象宏观中的铜导线中，单个电子运动速度上亿米每秒，电场移动速度就是光速，但是电子在导线中的整体移动速度只有几厘米每秒）半个纳米有多长？苯环的直径就接近半个纳米了吧（此处用电子云密度等值面的方法确定分子直径而不是H到H原子核的距离）。所以“无论外加电场如何变化，分子内部的电子云密度总能在一瞬间弛豫完毕，根本不需要弛豫时间”并不总是合理的！它的合理区间应该是激发光的波长超过2000nm甚至更长的波长！但是我并不知道如何去解决它！

话说我这么坦诚地袒露自己当前工作的不足，不给个威望值啥的鼓励下

如果点进来的人不幸发现不知道我在说什么，你就当我是一本正经的胡说八道好啦

作者
Author: sobereva 时间: 2015-9-1 19:33
daug-cc-pV5Z虽然算实际分子没什么实用性，但作为极小体系测试时的金标准不错。（说来Phys. Chem. Chem. Phys., 2001, 3, 4661-4666这篇用FF算稀有气体gamma，甚至用到骇人的t-aug-cc-pV6Z）

从静态的结果来看，def2+D算alpha和beta都很好，但gamma和金标准比差得太多，毕竟那也只是专门给计算极化率优化的，用到gamma勉为其难，还是sadlej+更可靠。

d-aug-cc-pVDZ其实也已经挺好了，计算量还没sadlaj+大，无论是静态还是动态，无论是alpha还是gamma，都已经达到了aug-cc-pV5Z的水准，和金标准daug-cc-pV5Z差异也不是太大，也在一定程度上表明高角动量极化函数用并不很大，而弥散函数的影响更为主要。另外可见d-aug-cc-pVDZ和sadlaj+结果相差很小，由于基组尺寸还成，所以体系大点的话，都可以考虑用d-aug-cc-pVDZ或sadlej+。但毕竟这只是一个体系的测试，还不足以完全说明问题，以后可以测试更多的。（之后测试可以顺便记录下实际计算耗时，因为不同角动量的GTF的积分计算时间相差不少，GTF数和实际耗时不那么好挂钩）。

这个测试很有意义，不妨之后扩展一下发个文章。

pople是根本没法用，别说极化率、超极化率了，即便加上弥散函数后连偶极矩都算不太准。这里顺便把之前做过的水的气相偶极矩测试数据列在这（实验值是1.8546 Debye），算是对alpha/beta/gamma这一系列不同阶数的响应性质的补充：
B3LYP/aug-cc-pVTZ：1.8471
MP2/aug-cc-pVTZ：1.8591
B3LYP/cc-pVTZ：1.9191
HF/aug-cc-pVTZ：1.94
B3LYP/sadlej-pVTZ：1.8618
B3LYP/def2-TZVPD：1.8508
B3LYP/6-311++G**：2.1591 （居然误差比6-31G*还大！虽然不免巧合）
B3LYP/6-31G*：2.0952

GTF数：
aug-cc-pVTZ：126 GTFs
sadlej-pVTZ：88 GTFs
def2-TZVPD：83 GTFs
6-311++G**：54 GTFs
可见sadlej-pVTZ和def2-TZVPD算低阶响应属性都很出色，比aug-cc-pVnZ划算得多。

作者
Author: superrice 时间: 2015-9-1 22:28
一般来说不都假设在激光作用下，分子来不及弛豫吗。。

作者
Author: yjcmwgk 时间: 2015-9-2 08:40

superrice 发表于 2015-9-1 22:28
一般来说不都假设在激光作用下，分子来不及弛豫吗。。

原子核当然不允许移动位置，电子云当然必须移动位置啊

作者
Author: superrice 时间: 2015-9-2 22:36

yjcmwgk 发表于 2015-9-2 08:40
原子核当然不允许移动位置，电子云当然必须移动位置啊

一般都加个damping factor吧，这个值一般都是实验拟合的吧，我貌似没看到过文献计算这个值的。有时候不用太纠结我觉得，有的体系确实不好算。我们旁边有个组，做的超高精度方法的开发，看他们发的JCP的文章就像天书一样，但实用性几乎为0，根本不能应用实际体系，感觉现在计算机硬件还差得远。

作者
Author: yjcmwgk 时间: 2015-9-3 10:14

superrice 发表于 2015-9-2 22:36
一般都加个damping factor吧，这个值一般都是实验拟合的吧，我貌似没看到过文献计算这个值的。有时候不用 ...

我在第三条说的就是这个damping factor呀。
真是没有好办法么

作者
Author: superrice 时间: 2015-9-3 23:07

yjcmwgk 发表于 2015-9-3 10:14
我在第三条说的就是这个damping factor呀。
真是没有好办法么

个人认为没啥好办法=。=

作者
Author: 我本是个娃娃 时间: 2015-9-16 14:09

sobereva 发表于 2015-9-1 19:33
daug-cc-pV5Z虽然算实际分子没什么实用性，但作为极小体系测试时的金标准不错。（说来Phys. Chem. Chem. Ph ...

建议老师把超极化率这一块儿的知识单独写一个帖子，详详细细地讨论一下，或者在明年初的书里好好地写上一笔

作者
Author: stecue 时间: 2015-9-27 22:59
我对这块不熟，但是觉得用DFT算非线性光学原理上就有点奇怪。即使是TDDFT本身就是线性响应框架内的，为什么可以计算非线性光学响应呢？

欢迎光临计算化学公社 (http://bbs.keinsci.com/)