计算化学公社

 找回密码 Forget password
 注册 Register
Views: 6145|回复 Reply: 22

[综合交流] 使用adcc和pyscf计算ADC级别的单三线态能级(for MR-TADF)

[复制链接 Copy URL]

877

帖子

36

威望

4803

eV
积分
6400

Level 6 (一方通行)

发表于 Post on 2022-3-22 01:34:52 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 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
bn.png
结构图所示,该分子具有一定的MR-TADF结构特征

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

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

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

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


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







BN.xyz

1.25 KB, 下载次数 Times of downloads: 9

mrtadf.xlsx

27.07 KB, 下载次数 Times of downloads: 4

评分 Rate

参与人数
Participants 6
威望 +1 eV +25 收起 理由
Reason
ene + 5
lonemen + 5 钟老师大法好!
丁越 + 5 赞!
hebrewsnabla + 5 精品内容
sobereva + 1
卡开发发 + 5 精品内容

查看全部评分 View all ratings

1464

帖子

1

威望

2591

eV
积分
4075

Level 6 (一方通行)

喵星人

发表于 Post on 2022-3-22 03:01:10 | 显示全部楼层 Show all
本帖最后由 喵星大佬 于 2022-3-22 03:03 编辑

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

应该是如果要算单重态

4万

帖子

99

威望

4万

eV
积分
89888

管理员

公社社长+计算化学玩家

发表于 Post on 2022-3-22 09:04:28 | 显示全部楼层 Show all
最近有文章声称范围分离双杂化结合SCS-ADC(2)算各类激发能都很理想:https://doi.org/10.1021/acs.jctc.1c01100,MRCC能做,可以考虑测一测
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班,内容介绍以及往届资料购买请点击链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的最佳途径。培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

1158

帖子

1

威望

2797

eV
积分
3975

Level 5 (御坂)

发表于 Post on 2022-3-22 09:14:25 | 显示全部楼层 Show all
喵星大佬 发表于 2022-3-22 03:01
这玩意的cas和nevpt2的(88)是怎么挑出来的。。。。感觉要14,14
14,14的话cas倒是问题不大但我觉得nevpt2 ...

14,14能mrpt2算动的

1158

帖子

1

威望

2797

eV
积分
3975

Level 5 (御坂)

发表于 Post on 2022-3-22 09:16:29 | 显示全部楼层 Show all
之前看的ADC2对一般性体系测得benchmark,感觉精度还是很拉胯

877

帖子

36

威望

4803

eV
积分
6400

Level 6 (一方通行)

 楼主 Author| 发表于 Post on 2022-3-22 09:22:59 | 显示全部楼层 Show all
本帖最后由 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 之内
简单测试, 就没有搞更大的空间了...

1464

帖子

1

威望

2591

eV
积分
4075

Level 6 (一方通行)

喵星人

发表于 Post on 2022-3-22 09:28:54 | 显示全部楼层 Show all
本帖最后由 喵星大佬 于 2022-3-22 09:32 编辑
ggdh 发表于 2022-3-22 09:22
取的mp2的自然轨道 占据0.05 和1.95 之内
简单测试, 就没有搞更大的空间了...

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

877

帖子

36

威望

4803

eV
积分
6400

Level 6 (一方通行)

 楼主 Author| 发表于 Post on 2022-3-22 09:32:56 | 显示全部楼层 Show all
喵星大佬 发表于 2022-3-22 09:28
感觉这个办法并不那么靠谱
难道不是要把s1和t1激发的主要轨道都扔进去嘛
照理说NEVPT2应该是这对里 ...

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

1464

帖子

1

威望

2591

eV
积分
4075

Level 6 (一方通行)

喵星人

发表于 Post on 2022-3-22 09:37:26 | 显示全部楼层 Show all
ggdh 发表于 2022-3-22 09:32
请问怎么确定s1和t1激发的主要轨道有哪些?

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

877

帖子

36

威望

4803

eV
积分
6400

Level 6 (一方通行)

 楼主 Author| 发表于 Post on 2022-3-22 09:51:15 | 显示全部楼层 Show all
sobereva 发表于 2022-3-22 09:04
最近有文章声称范围分离双杂化结合SCS-ADC(2)算各类激发能都很理想:https://doi.org/10.1021/acs.jctc.1c0 ...

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

422

帖子

7

威望

4764

eV
积分
5326

Level 6 (一方通行)

BSJ Institute

发表于 Post on 2022-3-22 17:12:17 | 显示全部楼层 Show all
分享一个去年做的常见计算方法在一个有机荧光分子(结构忘了,大概是个磷杂芳环)的S1, T1能级的测评:
QQ图片20220322170928.png 其中除了CIS(D)与实验值十分好地符合外,其他都表现很差。

评分 Rate

参与人数
Participants 1
eV +5 收起 理由
Reason
ggdh + 5 谢谢分享

查看全部评分 View all ratings

877

帖子

36

威望

4803

eV
积分
6400

Level 6 (一方通行)

 楼主 Author| 发表于 Post on 2022-3-22 18:03:34 | 显示全部楼层 Show all
Accelerator 发表于 2022-3-22 17:12
分享一个去年做的常见计算方法在一个有机荧光分子(结构忘了,大概是个磷杂芳环)的S1, T1能级的测评:
其 ...

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

272

帖子

0

威望

3937

eV
积分
4209

Level 6 (一方通行)

发表于 Post on 2022-3-22 20:27:43 | 显示全部楼层 Show all
Accelerator 发表于 2022-3-22 17:12
分享一个去年做的常见计算方法在一个有机荧光分子(结构忘了,大概是个磷杂芳环)的S1, T1能级的测评:
其 ...

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

312

帖子

1

威望

3106

eV
积分
3438

Level 5 (御坂)

发表于 Post on 2022-3-22 20:57:17 | 显示全部楼层 Show all
"道德水平在0.1eV"
哈哈哈哈哈哈....
我的道德水准不允许我用nm来比较激发能。

68

帖子

0

威望

1605

eV
积分
1673

Level 5 (御坂)

发表于 Post on 2022-3-27 21:11:29 | 显示全部楼层 Show all
有个问题哈,就是ADC是不是高估激发能,虽然detEST算得挺准得

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-2 23:16 , Processed in 0.689484 second(s), 26 queries .

快速回复 返回顶部 返回列表 Return to list