使用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文件

喵星大佬

这玩意的cas和nevpt2的(88)是怎么挑出来的。。。。感觉要14，14
14，14的话cas倒是问题不大但我觉得nevpt2应该做不动

如果要算三重态改成n_singlets=5就是算5个单重态

应该是如果要算单重态

sobereva

最近有文章声称范围分离双杂化结合SCS-ADC(2)算各类激发能都很理想：https://doi.org/10.1021/acs.jctc.1c01100，MRCC能做，可以考虑测一测

biogon

喵星大佬 发表于 2022-3-22 03:01
这玩意的cas和nevpt2的(88)是怎么挑出来的。。。。感觉要14，14
14，14的话cas倒是问题不大但我觉得nevpt2 ...

14，14能mrpt2算动的

biogon

之前看的ADC2对一般性体系测得benchmark，感觉精度还是很拉胯

ggdh

喵星大佬 发表于 2022-3-22 03:01
这玩意的cas和nevpt2的(88)是怎么挑出来的。。。。感觉要14，14
14，14的话cas倒是问题不大但我觉得nevpt2 ...

取的mp2的自然轨道 占据0.05 和1.95 之内
简单测试, 就没有搞更大的空间了...

喵星大佬

ggdh 发表于 2022-3-22 09:22
取的mp2的自然轨道 占据0.05 和1.95 之内
简单测试, 就没有搞更大的空间了...

感觉这个办法并不那么靠谱
难道不是要把s1和t1激发的主要轨道都扔进去嘛
照理说NEVPT2应该是这对里面最可靠的

ggdh

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

请问怎么确定s1和t1激发的主要轨道有哪些?

喵星大佬

ggdh 发表于 2022-3-22 09:32
请问怎么确定s1和t1激发的主要轨道有哪些?

看TDDFT的然后做HF的轨道里找成分差不多的？
@zjxitcc 专家来回答

ggdh

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

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

Accelerator

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

ggdh

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

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

chrinide

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

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

exity

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

心向暖阳

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