TDDFT-ris 2.0 更新，200原子体系无压力, 支持Multiwfn

JohnCase

ris方法原文：Zhou, Zehao, Fabio Della Sala, and Shane M. Parker.
                   Minimal auxiliary basis set approach for the electronic excitation spectra
                   of organic molecules. The Journal of Physical Chemistry Letters
                   14, no. 7 (2023): 1968-1976.
                    https://pubs.acs.org/doi/10.1021/acs.jpclett.2c03698

最近得空优化了TDDFT-ris的 toy code： John-zzh/pyscf_TDDFT_ris: efficient TDDFT-ris based on MOKIT and PySCF

优化方法：
1. 允许J项和K项使用不同的辅助基，默认J项使用sp辅助基提高精度，K项使用s辅助基，JK辅助基的轨道指数仍是一样的，都是0.2/r**2。
2. 对K项做了轨道截断处理，降低计算成本且不影响精度，默认删除比HOMO高>40eV的虚轨道，以及比LUMO低>40eV的占据轨道。
3. 完美支持Multiwfn做后续的激发态分析以及UV-vis光谱绘制。
4. 目前程序的效率支持~300原子大体系RKS-TDDFT-ris激发态计算，支持杂化泛函和范围分离杂化泛函。（TDA以及纯泛函的代码还没优化，有需求的朋友可以给我一些动力）
5. 支持单精度计算。

Jan15更新
1）已优化TDA和纯泛函。2）现在自动输出UV、ECD光谱到单独的文件。

有机体系速度： 117原子 wb97x-d/def2-TZVP/10 states ,      4核16G内存， 2min。
                       243原子 wb97x-d3bj/def2-TZVP/10 states， 4核 64G内存，12min。
                       362原子  wb97x-d3bj/def2-TZVP/10 states，4核 64G内存，1h 10min。
普通杂化泛函会比范围分离泛函便宜。


杂化泛函TDDFT-ris是N^4标度，在大体系上，AO->MO积分变换占了很可观的时间，其次是TDDFT矩阵向量乘。（当然，态越多矩阵向量乘也越贵）

纯泛函最便宜，TDDFT-ris方法的前身TDDFT-as方法，就是纯泛函，适合金属团簇体系，是N^3标度。
并行效率不高，主要是吃内存。因为基于pyscf的三中心积分变换理论上来说给的内存越大越快。当然，程序也允许小内存分批次转换。
以后有空继续优化代码，还有很多油水可以压榨。

========================================= =======
                       详细安装使用流程，及Multiwfn分析的步骤：
================================================
1. 首先下载程序到本地

1. git clone git@github.com:John-zzh/pyscf_TDDFT_ris.git

复制代码

2. 设置Python环境变量(可以放进~/.bashrc)

1. export PYTHONPATH=absolue_path_to_ris_repo:$PYTHONPATH

复制代码

absolute_path_to_ris_repo是pyscf_TDDFT_ris的绝对目录

3. 安装PySCF，以及负责读取分子基态信息/格式转换的MOKIT, 这一步方法较多

-- 使用conda安装（推荐）

1. conda create -n ris-mokit-pyscf-py39 python=3.9
   
2. conda activate ris-mokit-pyscf-py39
   
3. conda install mokit pyscf -c mokit/label/cf -c conda-forge

复制代码

--下载MOKIT预编译版，并且使用pip安装Pyscf
安装MOKIT预编译版的详细描述
pip install pyscf

================== 此时安装完成 ==================

这里我用一个小例子做计算流程的演示，49原子的视黄醛分子。
先用g16做PBE0/def2-SVP的DFT基态计算，输入文件gau.gjf：

1. %chk=gau.chk
   
2. %nprocs=8
   
3. %mem=16GB
   
4. #p PBE1PBE/def2SVP  
   
5. 
   
6. retinal
   
7. 
   
8. 0 1
   
9. 
   
10. xyz略...
    
11. 
    
12. g16 gau.gjf   

复制代码

耗时44秒。

然后formchk gau.chk gau.fch，拿到fch文件。

==================== 接下来开始TDDFT-ris计算  =======================

使用conda安装MOKIT的要用conda激活当时的环境。如果是下载预编译的MOKIT， 则需要设置MOKIT的环境变量。
我这里用conda演示完整的环境变量设置，方便批量作业提交。

1. conda activate ris-mokit-pyscf-py39
   
2. export PYTHONPATH=absolue_path_to_ris_repo:$PYTHONPATH
   
3. export MKL_NUM_THREADS=4
   
4. export OMP_NUM_THREADS=4
   
5. python absolute_path_to_ris_repo/main.py -f gau.fch -func pbe0 -n 5 -M 4000 -pt 0.0001

复制代码

（export PYTHONPATH=absolue_path_to_ris_repo:$PYTHONPATH这一步如果已经做过了可以省略）
MKL_NUM_THREADS OMP_NUM_THREADS是OMP并行的咒语。MOKIT和Numpy都认这个。
-f gau.fch是刚刚拿到的fch文件。
-func pbe0 泛函名称需要提供，或者提供HF成份的系数-ax 0.25也可以，此处和基态计算的泛函参数一致即可。
-n 5 计算5个激发态
-M 8000 8G内存。其实程序不会严格限制内存，只是按照这个数字给积分变换分批次，否则大体系有可能撑爆节点内存。
-pt 0.0001 打印跃迁系数的阈值。做电子激发分析的话需要设置为0.0001。 默认是0.05。

注：如果使用其他软件，比如ORCA，做DFT基态计算，那么预先使用MOKIT把输出文件转换为fch文件即可
（假设你的gbw文件名为molecule.gbw）

1. orca_2mkl molecule -mkl
   
2. mkl2fch molecule.mkl molecule.fch

复制代码

TDDFT-ris计算2.17秒完成，

1. args.nroots 5
   
2. Woring directory: /home/jojo/ris/gaussian_multiwfn/TDDFT
   
3. nbf = 434
   
4. nif = 434
   
5. mo_coeff.shape (434, 434)
   
6. mo_energy.shape (434,)
   
7. Restricted Kohn-Sham
   
8. functional:  pbe0
   
9. get_mf_from_fch_time: 0.73 seconds
   
10. This job can use: 24CPUs
    
11. use new code
    
12. use different J and K fitting basis: J with sp and K with s
    
13. self.nroots 5
    
14. use single precision? True
    
15. loading default XC functional paramters from parameter.py
    
16. use hybrid XC functional
    
17. self.a_x 0.25
    
18. cartesian or spherical electron integral = _sph
    
19. n_occ = 78
    
20. n_vir = 356
    
21. MO truncation in K with threshold 40 eV above HOMO and below LUMO
    
22. rest_occ = 57
    
23. rest_vir = 169
    
24. ==================== RIJ ====================
    
25. Asigning minimal auxiliary basis set: sp
    
26. The exponent alpha: 0.2/R^2
    
27.    auxmol.cart = False
    
28.     Full eri3c in shape (np.int64(434), np.int64(434), np.int64(112)), will take 161 MB
    
29.     max_mem_for_one_batch 2000.0
    
30.     small 3c2e, just incore
    
31.     pre_T_pq time: 0.1 seconds
    
32. pre_T_pq_to_Tpq one batch time 0.0 seconds
    
33. T_ia_J time 0.5 seconds
    
34. ==================== RIK ====================
    
35. Asigning minimal auxiliary basis set: s
    
36. The exponent alpha: 0.2/R^2
    
37.    auxmol.cart = False
    
38.     Full eri3c in shape (np.int64(434), np.int64(434), np.int64(49)), will take 70 MB
    
39.     max_mem_for_one_batch 2000.0
    
40.     small 3c2e, just incore
    
41.     pre_T_pq time: 0.1 seconds
    
42. pre_T_pq_to_Tpq one batch time 0.0 seconds
    
43.     pre_T_pq time: 0.0 seconds
    
44. pre_T_pq_to_Tpq one batch time 0.0 seconds
    
45.     pre_T_pq time: 0.0 seconds
    
46. pre_T_pq_to_Tpq one batch time 0.0 seconds
    
47. T_ia_K T_ij_K T_ab_K time 0.3 seconds
    
48. ======= TDDFT Eigen Solver Statrs =======
    
49. size of A matrix = 27768
    
50. Using non-orthogonalized Krylov subspace (nKs) method.
    
51. 
    
52.         Furche, Filipp, Brandon T. Krull, Brian D. Nguyen, and Jake Kwon.
    
53.         Accelerating molecular property calculations with nonorthonormal Krylov space methods.
    
54.         The Journal of Chemical Physics 144, no. 17 (2016).
    
55. 
    
56. iter: 1    max|R|: 1.27e-01   new_vectors: 10  MVP: 0.1 seconds
    
57. iter: 2    max|R|: 3.37e-02   new_vectors: 5  MVP: 0.1 seconds
    
58. iter: 3    max|R|: 1.28e-02   new_vectors: 5  MVP: 0.1 seconds
    
59. iter: 4    max|R|: 5.19e-03   new_vectors: 5  MVP: 0.1 seconds
    
60. iter: 5    max|R|: 1.70e-03   new_vectors: 4  MVP: 0.1 seconds
    
61. iter: 6    max|R|: 9.69e-04   new_vectors: 1  MVP: 0.0 seconds
    
62. Finished in 6 steps, 0.72 seconds
    
63. final subspace = 30
    
64. max_norm = 9.69e-04
    
65. MVcost     0.4901s 67.86%
    
66. fill_holder_cost 0.0892s 12.35%
    
67. subgencost 0.0266s 3.69%
    
68. subcost    0.0267s 3.69%
    
69. full_cost  0.0362s 5.01%
    
70. ======= TDDFT Eigen Solver Done =======
    
71. check norm of X^TX - Y^YY - I = 4.38e-08
    
72. print RKS transition coefficients larger than 1.00e-01
    
73. index of HOMO: 78
    
74. index of LUMO: 79
    
75. Excited State     1:      Singlet-A      2.9053 eV  427.04 nm  f=0.0001   <S**2>=0.000
    
76.              76 -> 79               0.65193
    
77.              76 -> 80              -0.23655
    
78. Excited State     2:      Singlet-A      2.9834 eV  415.87 nm  f=1.3705   <S**2>=0.000
    
79.              77 -> 79              -0.14597
    
80.              78 -> 79               0.69388
    
81. Excited State     3:      Singlet-A      3.8900 eV  318.95 nm  f=0.4812   <S**2>=0.000
    
82.              77 -> 79              -0.64082
    
83.              78 -> 79              -0.15515
    
84.              78 -> 80              -0.23941
    
85. Excited State     4:      Singlet-A      4.3099 eV  287.87 nm  f=0.0396   <S**2>=0.000
    
86.              75 -> 79               0.35318
    
87.              77 -> 79               0.19885
    
88.              78 -> 80              -0.57193
    
89. Excited State     5:      Singlet-A      4.9689 eV  249.69 nm  f=0.0002   <S**2>=0.000
    
90.              76 -> 79               0.27008
    
91.              76 -> 80               0.59980
    
92.              76 -> 81              -0.19345
    
93.              76 -> 83              -0.11782
    
94. transition coefficient data also written to gau_TDDFT_ris_coeff_Multiwfn.txt
    
95. ================================================
    
96. eV       nm       cm^-1    oscillator strength
    
97. 2.905    427      23433    0.00012678
    
98. 2.983    416      24063    1.37051890
    
99. 3.890    319      31375    0.48124768
    
100. 4.310    288      34762    0.03955060
     
101. 4.969    250      40077    0.00018934
     
102. spectra data also written to gau_TDDFT_ris_eV_os_Multiwfn.txt
     
103. total ris time: 2.48 seconds
     
104. 
     
105.     Please cite the TDDFT-ris method:
     
106. 
     
107.         1.  Zhou, Zehao, Fabio Della Sala, and Shane M. Parker.
     
108.             Minimal auxiliary basis set approach for the electronic excitation spectra
     
109.             of organic molecules. The Journal of Physical Chemistry Letters
     
110.             14, no. 7 (2023): 1968-1976.
     
111.             (must cite)
     
112. 
     
113.         2.  Zhou, Zehao, and Shane M. Parker.
     
114.             Converging Time-Dependent Density Functional Theory Calculations in Five Iterations
     
115.             with Minimal Auxiliary Preconditioning. Journal of Chemical Theory and Computation
     
116.             20, no. 15 (2024): 6738-6746.
     
117.             (for efficient orbital truncation technique)
     
118. 
     
119.         2.  Giannone, Giulia, and Fabio Della Sala.
     
120.             Minimal auxiliary basis set for time-dependent density functional theory and
     
121.             comparison with tight-binding approximations: Application to silver nanoparticles.
     
122.             The Journal of Chemical Physics 153, no. 8 (2020).
     
123.             (The idea of TDDFT-ris originates from TDDFT-as)
     
124. 
     
125.     And cite the pyscf-TDDFT-ris package:
     
126. 
     
127.         1. Zehao Zhou, pyscf-TDDFT-ris, (https://github.com/John-zzh/pyscf_TDDFT_ris)
     
128. 
     
129.         2. Jingxiang Zou, Molecular Orbital Kit (MOKIT), (https://gitlab.com/jxzou/mokit)
     
130. 
     
131.         3. PySCF: the Python-based simulations of chemistry framework,
     
132.            Q. Sun, et. al., and G. K.-L. Chan, WIREs Comput. Mol. Sci. 8, e1340 (2018)
     
133.            (https://pyscf.org/about.html)

复制代码

除了屏幕上打印的信息，在当前目录还生成了3个文件，
gau_TDDFT_ris_eV_os_Multiwfn.txt，  激发能/振子强度信息，用于Multiwfn绘制UV-vis激发光谱
gau_TDDFT_ris_eV_rs_Multiwfn.txt，  激发能/转子强度信息，用于Multiwfn绘制ECD光谱
gau_TDDFT_ris_coeff_Multiwfn.txt，    跃迁系数信息，在Multiwfn做电子激发分析的步骤中代替高斯的log/out文件即可


============================ Multiwfn 绘制TDDFT-ris 输出的UV-vis、ECD光谱  ===========================
打开Multiwfn，直接读取gau_TDDFT_ris_eV_os_Multiwfn.txt文件，
选择
11 Plot IR/Raman/UV-Vis/ECD/VCD/ROA/NMR spectrum，
3:UV-Vis
0 Plot spectrum
即可绘制激发光谱。
我这里和g16 ab-initio PBE0/def-2SVP的计算结果做对比，
可以写个multiple.txt文件，让Multiwfn第一步直接读取multiple.txt，即可同时绘制两者对比图
gau.log                                            pbe0-def2SVP
gau_TDDFT_ris_eV_os_Multiwfn.txt   pbe0-ris-def2SVP

可见TDDFT-ris和精确的TDDFT光谱非常接近！这打破了半经验方法的传统印象。

同理，读取gau_TDDFT_ris_eV_rs_Multiwfn.txt文件即可绘制ECD光谱。

============================ Multiwfn 做NTO分析  ===========================

打开Multiwfn，先读取gau.fch文件，
选择
18 Electron excitation analysis
6 Generate natural transition orbitals (NTOs)
输入gau_TDDFT_ris_coeff_Multiwfn.txt文件路径
选择感兴趣的态 2
3 Output NTO orbitals to .mwfn file (recommended)
D:\tmp\test.wfn
此时退出并重新打开Multiwfn，加载D:\tmp\test.wfn即可观看NTO轨道

（空穴-电子等分析方法同理可得）

========================== 演示结束，使用愉快 ==========================
有任何使用问题/开发建议随时联系，也可以合作计算大批量激发态任务。

yzh

请问老师发射光谱（PL）的结果怎么样，对比性如何？谢谢

JohnCase

yzh 发表于 2024-12-24 22:45
请问老师发射光谱（PL）的结果怎么样，对比性如何？谢谢

荧光吗？目前还没有开发梯度，因此没有这个功能（以后有空再做

wzkchem5

对于足够大的体系应该还是O(N^4)标度吧，因为RI-K是O(N^4)的scaling？有没有考虑用domain fitting之类的手段降低scaling？

JohnCase

wzkchem5 发表于 2024-12-25 00:42
对于足够大的体系应该还是O(N^4)标度吧，因为RI-K是O(N^4)的scaling？有没有考虑用domain fitting之类的手 ...

确实，我看了下详细的时间，即便截断了虚轨道，RIK也还是比RIJ更贵。

JohnCase

wzkchem5 发表于 2024-12-25 00:42
对于足够大的体系应该还是O(N^4)标度吧，因为RI-K是O(N^4)的scaling？有没有考虑用domain fitting之类的手 ...

我有空去瞅瞅domain fitting。

yzh

JohnCase 发表于 2024-12-24 23:01
荧光吗？目前还没有开发梯度，因此没有这个功能（以后有空再做

谢谢，期待老师的荧光

smooth85

按照了流程测试了一下，代码很好用，计算速度很快。感谢作者的付出~

有两个问题：
1. 虽然设置了MKL_NUM_THREADS 和 OMP_NUM_THREADS，但是CPU依然只用了不到一半在跑。

2. 是否能做ECD ？

另外，期待TDA代码的加入

再次感谢开发者~

JohnCase

smooth85 发表于 2024-12-26 08:56
按照了流程测试了一下，代码很好用，计算速度很快。感谢作者的付出~

有两个问题：

我自己测试2-4核以后可能就没有更多收获了，毕竟主要部分是numpy计算。
可以做ECD，但是相关的benchmark工作还没有，精度是不确定的。
TDA代码其实是有的，只是这次改进的时候这部分代码没翻新。。。如果只是看激发光谱的话TDDFT比TDA的振子强度要好不少。我可能元旦以后才有空把TDA这部分也优化一下~

smooth85

感谢回复，期待更新

JohnCase

smooth85 发表于 2024-12-26 08:56
按照了流程测试了一下，代码很好用，计算速度很快。感谢作者的付出~

有两个问题：

现在支持ECD了