使用adcc和pyscf计算ADC级别的单三线态能级(for MR-TADF)

ggdh

类似的部分程序(待补充):
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

结构图所示,该分子具有一定的MR-TADF结构特征

2.2 创建pyscf的输入文件BN.py, 其中内容如下

1. from pyscf import gto, scf, lib
   
2. import adcc
   
3. 
   
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()
   
9. 
   
10. adcc.set_n_threads(12)
    
11. state = adcc.adc2(scfres, n_triplets=1)
    
12. print(state.describe())
    
13. print(state.describe_amplitudes())

复制代码

其中需要注意的是:

• 12为openmp并行核数, 这里并不支持mpi并行, 要分别为pyscf和adcc设
• BN.xyz为xyz文件名
• cc-pvdz为基组名
• n_triplets=1的意思是算一个三重态, 如果要算单重态改成n_singlets=5就是算5个单重态
  

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
分别对这几个参数做了测试, 耗时和精度下降的情况入下图所示
其中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
   
3. 
   
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

一些讨论:

• adc2得到的结果和用orca的STEOM-DLPNO-CCSD得到的结果近似
• 纯泛函和HF%小的泛函(B3LYP)会因为低估T1而产生大的ΔEST
• 高HF%的杂化泛函/范围分离泛函/双杂化泛函对T1的估算比较准, 但是会因为高估S1而产生大的ΔEST
• 双杂化泛函表现并不好, 但部分泛函加了SOS后情况得到了好转其中表现最好的是SOS-PBE-QIDH和SOS-RSX-QIDH
• 注意上述结论仅仅针对当前一个分子得到, 对其他体系是否适用还需大家自行测试
  

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

ggdh

心向暖阳 发表于 2022-4-9 20:45
我测试了几个分子，觉得双杂化算出来的Est为负，可能跟波函数不稳定有关，导致结果乱七八糟的

还有就是 ...

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

心向暖阳

wzkchem5 发表于 2022-4-11 14:31
对，如果只算N个激发态，有很大概率得到的不是前N个激发态，而是得到了一些高的却漏掉了低的。非双杂化泛 ...

谢谢爱心苹果老师

wzkchem5

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

所以CIS(D) ...

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

心向暖阳

wzkchem5 发表于 2022-4-10 23:08
激发态排序是乱的，是因为激发态是按doubles correction之前的能量排的，做完doubles correction以后为了 ...

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

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

心向暖阳

ggdh 发表于 2022-4-10 22:06
马上有mrcc的攻略了, 那个才是真的C, 敬请期待

tql

wzkchem5

心向暖阳 发表于 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

ggdh

心向暖阳 发表于 2022-4-9 20:45
我测试了几个分子，觉得双杂化算出来的Est为负，可能跟波函数不稳定有关，导致结果乱七八糟的

还有就是 ...

马上有mrcc的攻略了, 那个才是真的C, 敬请期待

心向暖阳

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

还有就是我的激发能居然算出来排序结果都是乱的，奇怪。

心向暖阳

有个问题哈，就是ADC是不是高估激发能，虽然detEST算得挺准得

exity

"道德水平在0.1eV"
哈哈哈哈哈哈....
我的道德水准不允许我用nm来比较激发能。

chrinide

Accelerator 发表于 2022-3-22 17:12
分享一个去年做的常见计算方法在一个有机荧光分子（结构忘了，大概是个磷杂芳环）的S1, T1能级的测评：
其 ...

个人认为（不负责任滴）大概率实验测得的“实验值”是一个不那么准确的“统计平均”值；而理论计算的结果只是（真空下？）单个最低能量构象下的理论值。除非你的实验极端高精尖--BEC（玻色爱因斯坦凝聚）下单分子束缚下测量；或者激光冷却/光腔束缚单分子测量；分子束缚势阱（把小分子注入C60-C180等巴基球分子中）下测得的

ggdh

Accelerator 发表于 2022-3-22 17:12
分享一个去年做的常见计算方法在一个有机荧光分子（结构忘了，大概是个磷杂芳环）的S1, T1能级的测评：
其 ...

DSD-PBEP86这个离谱啊, 我刚算了算上面的体系, deltaEST都是负的
另外你这个里面CC的方法只有第一个, 但是S1和T1用的方法不完全一样, 结果让人有点怀疑.
其他泛函的趋势好像和我测试的差不多.

Accelerator

分享一个去年做的常见计算方法在一个有机荧光分子（结构忘了，大概是个磷杂芳环）的S1, T1能级的测评：
其中除了CIS(D)与实验值十分好地符合外，其他都表现很差。

ggdh

sobereva 发表于 2022-3-22 09:04
最近有文章声称范围分离双杂化结合SCS-ADC(2)算各类激发能都很理想：https://doi.org/10.1021/acs.jctc.1c0 ...

这是上周5才放出的MRCC新版才加入的新功能把, sob的消息好灵通啊