计算化学公社

标题: 使用adcc和pyscf计算ADC级别的单三线态能级(for MR-TADF) [打印本页]
作者
Author: ggdh    时间: 2022-3-22 01:34
标题: 使用adcc和pyscf计算ADC级别的单三线态能级(for MR-TADF)
本帖最后由 ggdh 于 2022-3-22 18:54 编辑

类似的部分程序(待补充):
Turbomole: 支持ri加速的scs-cc2, 但是要收费
eT: 有multilevel-cc加速, 但不支持三重态计算
orca: ADC不支持计算三重态, DLPNO-STEOM-CCSD 速度很快, 但容易出现收敛问题
MRCC: 待测试

1. 程序安装
1.1 安装anaconda或者miniconda
1.2 安装pyscf和adcc
  1. conda create -n pyscf -c pyscf pyscf
  2. conda activate pyscf
  3. conda install -c conda-forge adcc
复制代码
另外也可以选择别的方法比如pip, 参考:
https://adc-connect.org/v0.15.13/installation.html
https://pyscf.org/install.html

2. 输入文件的准备
2.1 准备一个优化好的xyz文件, 手边没有现成的可以用附件中的xyz文件, 其他格式的结构文件可以用Multiwfn (主功能100里的子功能2)转成xyz
(, 下载次数 Times of downloads: 35)
结构图所示,该分子具有一定的MR-TADF结构特征

2.2 创建pyscf的输入文件BN.py, 其中内容如下
  1. from pyscf import gto, scf, lib
  2. import adcc

  4. lib.num_threads(12)
  5. # Run SCF in pyscf
  6. mol = gto.M(atom=gto.mole.fromfile("BN.xyz"), basis='cc-pvdz')
  7. scfres = scf.RHF(mol)
  8. scfres.kernel()

  10. adcc.set_n_threads(12)
  11. state = adcc.adc2(scfres, n_triplets=1)
  12. print(state.describe())
  13. print(state.describe_amplitudes())
复制代码
其中需要注意的是:


2.3运行程序
  1. export PYSCF_TMPDIR=/tmp
  2. conda activate pyscf
  3. python BN.py | tee BN.out
复制代码
其中第一行是设缓存目录, 大体系很耗硬盘, 所以如果这里尽量设一个体积大的硬盘位置
然后部分输出结果如下:
  1. converged SCF energy = -539.461545133644
  2. Starting adc2 triplet Jacobi-Davidson ...
  3. Niter n_ss  max_residual  time  Ritz values
  4.   1     4       0.32167  23.3s  [0.29368883]
  5.   2     8      0.039105  12.7s  [0.13900202]
  6.   3    12     0.0017295  21.4s  [0.1199706]
  7.   4    16    0.00093827  24.2s  [0.11816422]
  8.   5    20    0.00016175  24.8s  [0.11744624]
  9. === Converged ===
  10.     Number of matrix applies:    20
  11.     Total solver time:           106s 519ms
  12. +------------------------------------------------------------------------+
  13. | adc2                                    triplet ,  converged               |
  14. +------------------------------------------------------------------------+
  15. |  #        excitation energy     osc str    |v1|^2    |v2|^2       |
  16. |          (au)           (eV)                                                     |
  17. |  0     0.1174462      3.195875   0.0000    0.9119   0.08811  |
  18. +------------------------------------------------------------------------+
复制代码
可以看到其中有adc2的收敛情况, 激发态能级, 振子强度, 单激发态和双激发的比例

3. 加速的技巧
加速通常意味着损失精度或者有产生误差风险, 但是postHF往往会面对体系太大无法计算的情况(用这个程序必须保证内存必须要大)
所以往往不得不放弃部分精度来完成计算. 这其中涉及到的道德,IF,良心,电费等因素还需用户自行权衡.
首先是可以改变的是基组的大小和体系中非共轭部分的裁剪, 这就不多说了.
然后影响速度的参数主要是如下几个参数:
state = adcc.adc2(scfres, n_triplets=1, frozen_core=0,frozen_virtual=0, n_guesses=4, conv_tol=1e-6)
其中:
n_triplets/n_singlets: 所需计算的态的数量
frozen_core: 冻结内层轨道的数量, 从低往高数
frozen_virtual: 冻结空轨道的数量, 从高往低数
n_guesses: 计算的猜的数量, 默认是激发态的数量的2倍, 但不会小于4
conv_tol: 收敛标准, 默认是SCF收敛标准的10倍但不会小于1e-6, 上面的例子中这个值大概是在1e-4
分别对这几个参数做了测试, 耗时和精度下降的情况入下图所示 (, 下载次数 Times of downloads: 17)
其中time列是耗时(秒), eV列是激发态跃迁能, ΔeV列代表精度的损耗,
假设道德底线是0.01 eV, 假设其他体系也有类似的规律, 那么从图中可得得结论如下
收敛标准conv_tol对耗时影响较大, 但对精度影响不大1e-4(默认)都能接受,保险的话设1e-5(如果要提高到1e-5, 注意必须得另外把scf的精度调高到1e-6)
猜数量n_guesses对耗时影响大, 但对精度影响不大(这个应该要看体系, 改小这个参数一定要小心)
冻结内层电子frozen_core对耗时影响大, 如果只冻结内层电子则对精度影响不大(这个体系的内层轨道--1S轨道--正好是14个)
冻结高层空轨道frozen_virtual对耗时影响没那么大, 但由于空轨道数量多, 总体对速度影响还是大, 但是对精度影响较大(0.01eV级别)
综上所述,conv_tol用默认就够, n_guesses可以不放心的稍微改小, frozen_core可以比较放心的冻结内层轨道, frozen_virtual调大则需要耗费道德值.

4. 自动设加速参数的输入文件
根据第三节的测试, 我们发现conv_tol和n_guesses比较好设, 但是frozen_core和frozen_virtual比较让人纠结
所以我写了根据占据轨道和空轨道能量的高低来自动设frozen_core和frozen_virtual的输入文件如下
  1. from pyscf import gto, scf, lib
  2. import adcc

  4. occ_threash=3 # orbitals 3 Ha lower than HOMO will be frozen
  5. virtual_threash=5 # orbitals 3 Ha higher than LUMO will be frozen
  6. lib.num_threads(12)
  7. # Run SCF in pyscf
  8. mol = gto.M(atom=gto.mole.fromfile("BN.xyz"), basis='cc-pvdz')
  9. scfres = scf.RHF(mol)
  10. scfres.kernel()
  11. for i,o in enumerate(scfres.mo_occ):
  12.     if o == 0:
  13.         LUMO_idx = i
  14.         HOMO_idx = i-1
  15.         break
  16. occ_threash = scfres.mo_energy[HOMO_idx] - occ_threash
  17. virtual_threash = scfres.mo_energy[LUMO_idx] + virtual_threash
  18. fo = len([i for i in scfres.mo_energy if i < occ_threash])
  19. fv = len([i for i in scfres.mo_energy if i > virtual_threash])
  20. print("{:d} occ orbitals lower than {:.3f} eV will be frozen".format(fo,occ_threash*27.212))
  21. print("{:d} virtual orbitals higher than {:.3f} eV will be frozen".format(fv,virtual_threash*27.212))
  22. # Run adcc
  23. adcc.set_n_threads(12)
  24. state = adcc.adc2(scfres, n_singlets=1, n_guesses=2, conv_tol=1e-4, frozen_core=fo, frozen_virtual=fv)
  25. print(state.describe())
  26. print(state.describe_amplitudes())
复制代码
大家可以自行修改文件开头的occ_threash和virtual_threash值来设置根据能量来固定多少轨道, 如果道德水平在0.01eV, 我估计可以设4,4 (比HOMO低4Ha 比LUMO高4Ha的轨道纳入冻结), 如果道德水平在0.1eV,我估计可以设3,1. 大家可以弄个小体系自行尝试.

5. 一个小的benchmark
(, 下载次数 Times of downloads: 26)
一些讨论:

附件中包含了excel数据,以及分子结构xyz文件







作者
Author: 喵星大佬    时间: 2022-3-22 03:01
本帖最后由 喵星大佬 于 2022-3-22 03:03 编辑

这玩意的cas和nevpt2的(88)是怎么挑出来的。。。。感觉要14，14
14，14的话cas倒是问题不大但我觉得nevpt2应该做不动
如果要算三重态改成n_singlets=5就是算5个单重态

应该是如果要算单重态

作者
Author: sobereva    时间: 2022-3-22 09:04
最近有文章声称范围分离双杂化结合SCS-ADC(2)算各类激发能都很理想：https://doi.org/10.1021/acs.jctc.1c01100，MRCC能做，可以考虑测一测
作者
Author: biogon    时间: 2022-3-22 09:14
喵星大佬 发表于 2022-3-22 03:01
这玩意的cas和nevpt2的(88)是怎么挑出来的。。。。感觉要14，14
14，14的话cas倒是问题不大但我觉得nevpt2 ...

14，14能mrpt2算动的
作者
Author: biogon    时间: 2022-3-22 09:16
之前看的ADC2对一般性体系测得benchmark，感觉精度还是很拉胯
作者
Author: ggdh    时间: 2022-3-22 09:22
本帖最后由 ggdh 于 2022-3-22 09:25 编辑
喵星大佬 发表于 2022-3-22 03:01
这玩意的cas和nevpt2的(88)是怎么挑出来的。。。。感觉要14，14
14，14的话cas倒是问题不大但我觉得nevpt2 ...

取的mp2的自然轨道 占据0.05 和1.95 之内
简单测试, 就没有搞更大的空间了...
作者
Author: 喵星大佬    时间: 2022-3-22 09:28
本帖最后由 喵星大佬 于 2022-3-22 09:32 编辑
ggdh 发表于 2022-3-22 09:22
取的mp2的自然轨道 占据0.05 和1.95 之内
简单测试, 就没有搞更大的空间了...

感觉这个办法并不那么靠谱
难道不是要把s1和t1激发的主要轨道都扔进去嘛
照理说NEVPT2应该是这对里面最可靠的
作者
Author: ggdh    时间: 2022-3-22 09:32
喵星大佬 发表于 2022-3-22 09:28
感觉这个办法并不那么靠谱
难道不是要把s1和t1激发的主要轨道都扔进去嘛
照理说NEVPT2应该是这对里 ...

请问怎么确定s1和t1激发的主要轨道有哪些?
作者
Author: 喵星大佬    时间: 2022-3-22 09:37
ggdh 发表于 2022-3-22 09:32
请问怎么确定s1和t1激发的主要轨道有哪些?

看TDDFT的然后做HF的轨道里找成分差不多的？
@zjxitcc 专家来回答
作者
Author: ggdh    时间: 2022-3-22 09:51
sobereva 发表于 2022-3-22 09:04
最近有文章声称范围分离双杂化结合SCS-ADC(2)算各类激发能都很理想：https://doi.org/10.1021/acs.jctc.1c0 ...

这是上周5才放出的MRCC新版才加入的新功能把, sob的消息好灵通啊
作者
Author: Accelerator    时间: 2022-3-22 17:12
分享一个去年做的常见计算方法在一个有机荧光分子（结构忘了，大概是个磷杂芳环）的S1, T1能级的测评：
(, 下载次数 Times of downloads: 21) 其中除了CIS(D)与实验值十分好地符合外，其他都表现很差。
作者
Author: ggdh    时间: 2022-3-22 18:03
Accelerator 发表于 2022-3-22 17:12
分享一个去年做的常见计算方法在一个有机荧光分子（结构忘了，大概是个磷杂芳环）的S1, T1能级的测评：
其 ...

DSD-PBEP86这个离谱啊, 我刚算了算上面的体系, deltaEST都是负的
另外你这个里面CC的方法只有第一个, 但是S1和T1用的方法不完全一样, 结果让人有点怀疑.
其他泛函的趋势好像和我测试的差不多.
作者
Author: chrinide    时间: 2022-3-22 20:27
Accelerator 发表于 2022-3-22 17:12
分享一个去年做的常见计算方法在一个有机荧光分子（结构忘了，大概是个磷杂芳环）的S1, T1能级的测评：
其 ...

个人认为（不负责任滴）大概率实验测得的“实验值”是一个不那么准确的“统计平均”值；而理论计算的结果只是（真空下？）单个最低能量构象下的理论值。除非你的实验极端高精尖--BEC（玻色爱因斯坦凝聚）下单分子束缚下测量；或者激光冷却/光腔束缚单分子测量；分子束缚势阱（把小分子注入C60-C180等巴基球分子中）下测得的
作者
Author: exity    时间: 2022-3-22 20:57
"道德水平在0.1eV"
哈哈哈哈哈哈....
我的道德水准不允许我用nm来比较激发能。
作者
Author: 心向暖阳    时间: 2022-3-27 21:11
有个问题哈，就是ADC是不是高估激发能，虽然detEST算得挺准得
作者
Author: 心向暖阳    时间: 2022-4-9 20:45
本帖最后由 心向暖阳 于 2022-4-9 20:53 编辑

我测试了几个分子，觉得双杂化算出来的Est为负，可能跟波函数不稳定有关，导致结果乱七八糟的

还有就是我的激发能居然算出来排序结果都是乱的，奇怪。
作者
Author: ggdh    时间: 2022-4-10 22:06
心向暖阳 发表于 2022-4-9 20:45
我测试了几个分子，觉得双杂化算出来的Est为负，可能跟波函数不稳定有关，导致结果乱七八糟的

还有就是 ...

马上有mrcc的攻略了, 那个才是真的C, 敬请期待
作者
Author: wzkchem5    时间: 2022-4-10 23:08
心向暖阳 发表于 2022-4-9 13:45
我测试了几个分子，觉得双杂化算出来的Est为负，可能跟波函数不稳定有关，导致结果乱七八糟的

还有就是 ...

激发态排序是乱的，是因为激发态是按doubles correction之前的能量排的，做完doubles correction以后为了方便和doubles correction之前的结果比较，没有重新排序。
Est为负不一定是错的，有的分子本来就是T1高于S1，这种分子最近两年挺热门的。可以参见https://pubs.acs.org/doi/10.1021/acs.jpca.1c10492
作者
Author: 心向暖阳    时间: 2022-4-11 08:51
ggdh 发表于 2022-4-10 22:06
马上有mrcc的攻略了, 那个才是真的C, 敬请期待

tql
作者
Author: 心向暖阳    时间: 2022-4-11 09:13
wzkchem5 发表于 2022-4-10 23:08
激发态排序是乱的，是因为激发态是按doubles correction之前的能量排的，做完doubles correction以后为了 ...

对于第一个问题哈，是不是需要多算一些激发态呢，然后以校正后的激发能自己确定激发态次序

所以CIS(D)是这么的香，

作者
Author: wzkchem5    时间: 2022-4-11 14:31
心向暖阳 发表于 2022-4-11 02:13
对于第一个问题哈，是不是需要多算一些激发态呢，然后以校正后的激发能自己确定激发态次序

所以CIS(D) ...

对，如果只算N个激发态，有很大概率得到的不是前N个激发态，而是得到了一些高的却漏掉了低的。非双杂化泛函的TDDFT也会有这个问题，但是远不如双杂化泛函严重，而且出现这个问题的原因也不一样（对于非双杂化泛函主要是因为Davidson方法不能保证收敛到最低的N个激发态；对于双杂化泛函这个因素也存在，但更主要的因素是doubles correction改变了激发态的顺序），所以非双杂化泛函的TDDFT一般只需要算大概N+3个根就可以比较好地保证得到了最低的N个激发态，但双杂化泛函经常需要多于这个数。
作者
Author: 心向暖阳    时间: 2022-4-11 15:32
wzkchem5 发表于 2022-4-11 14:31
对，如果只算N个激发态，有很大概率得到的不是前N个激发态，而是得到了一些高的却漏掉了低的。非双杂化泛 ...

谢谢爱心苹果老师
作者
Author: ggdh    时间: 2022-4-11 16:01
心向暖阳 发表于 2022-4-9 20:45
我测试了几个分子，觉得双杂化算出来的Est为负，可能跟波函数不稳定有关，导致结果乱七八糟的

还有就是 ...

能不能把你测的分子发出来看看(如果涉及到保密就算了)



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3