KirariNto：通过比较NTO重叠来实现激发态跟踪的程序

Accelerator

KirariNto是对Marco Campetella和Juan Sanz Garcia在10.1002/jcc.26162中提出的算法的实现，但相比而言使用的近似程度更低。原文作者虽然也开发了实现这一算法的程序，但并不公开。

KirariNto用于在激发态构型优化过程中尽可能维持激发态成分不变。其思想是在除第一步外的每一步，都先进行一个TD-DFT单点计算，对每个激发态分别进行NTO分析，与前一步选定的感兴趣的激发态的NTO轨道计算重叠积分，选取重叠程度最大的，作为下一步感兴趣的激发态并得到受力。在KirariNto中，这些过程均借助Multiwfn完成。为了适配Multiwfn的功能，计算两个电子态A和B之间的NTO重叠积分的具体做法如下：
1. 对连续两步的分子的几何结构进行对齐；
2. 对A和B分别得到NTO轨道和本征值；
3. 取出A和B各自的NTO-hole轨道和NTO-electron轨道，设置占据数为本征值的四次方根，共计保存成4个fchk文件；
4. 分别对这4个fchk文件计算电子密度，然后以hole-hole，electron-electron的方式得到重叠积分并相加。

KirariNto调用PySCF和geomeTRIC作为优化器，通过对Gaussian或ORCA的输出信息进行处理，得到每一步的能量和梯度并输入到优化器中。因此它需要如下依赖：
Linux操作系统；
PySCF；
geomeTRIC（可通过pip3 install geometric安装）
Numpy。
在运行前，需要指定与Multiwfn、ORCA、Gaussian相关的环境变量：

1. export KIRARINTO_MULTIWFN=~/Multiwfn_3.8_dev_bin_Linux_noGUI/Multiwfn
   
2.     export KIRARINTO_GAUSSIAN=~/g16/g16
   
3.     export KIRARINTO_FORMCHK=~/g16/formchk
   
4.     export KIRARINTO_ORCA=/opt/orca5/orca
   
5.     export KIRARINTO_ORCA2MKL=/opt/orca5/orca_2mkl

复制代码

运行方式如下：

1. usage: kirariNto.py [-h] [-c CONSTRAINTS] [-p PROGRAM] [-d DELETE] [-g GRID] input
   
2. 
   
3. positional arguments:
   
4.   input                 An input file for the TD-DFT calculation of your initial geometry.
   
5. 
   
6. optional arguments:
   
7.   -h, --help            show this help message and exit
   
8.   -c CONSTRAINTS, --constraints CONSTRAINTS (optional) An external geometry constraint file that follows the
   
9.                         format requirement of PySCF.
   
10.   -p PROGRAM, --program PROGRAM The program you are using (Gaussian or ORCA; case insensitive). Default: Gaussian if input file contains .gjf, and ORCA if input file contains .inp
    
11.   -d DELETE, --delete DELETE (optional) If the .chk or .gbw files to be deleted after each step. yes (default) or no.
    
12.   -g GRID, --grid GRID  (optional) 1, 2 (default) or 3. The accuracy of integral grid for overlap evaluation. Denser for a large value.
    

复制代码

需要首先准备一个用于对初始结构进行TD-DFT单点计算的输入文件，并指定nstates和nroot。例如（以下仅以展示格式为目的，请忽视其中计算水平）：
Gaussian:

1. %mem=8GB
   
2. %nprocshared=8
   
3. %chk=u_td_1.chk
   
4. # td=(<b>nstates=3,root=1</b>) b3lyp em=gd3bj iop(9/40=4) gen
   
5. 
   
6. TC
   
7. 
   
8. 0 1
   
9. C                  1.49653099   -0.03970163    0.25798558
   
10. C                  0.96661422    1.10990407   -0.32719382
    
11. C                 -0.29560541    1.04203217   -0.91629061
    
12. H                  2.46280152   -0.01555637    0.71693905
    
13. H                  1.51812645    2.02681494   -0.32424024
    
14. C                 -0.42910744   -1.21052831   -0.33465934
    
15. O                 -1.07773001   -2.28888230   -0.33812694
    
16. N                 -0.96370843   -0.11703076   -0.90600943
    
17. N                  0.78578354   -1.17302439    0.24032710
    
18. H                  1.16155785   -1.99994675    0.65864886
    
19. N                 -0.88665177    2.23323444   -1.54286970
    
20. H                 -1.88298815    2.18859514   -1.46992353
    
21. H                 -0.55480592    3.05391114   -1.07771587
    
22. 
    
23. -c -h -o -n 0
    
24. sto-3g
    
25. ****
    
26. 
    
27. 
    

复制代码

ORCA:

1. ! blyp def2-sv(p) def2/j
   
2. %pal nprocs 8 end
   
3. %tddft nroots 5 iroot 1 end
   
4. * xyz 0 1
   
5. C                  1.49653098   -0.03970163    0.25798558
   
6. C                  0.96661422    1.10990407   -0.32719382
   
7. C                 -0.29560541    1.04203217   -0.91629061
   
8. H                  2.46280151   -0.01555637    0.71693905
   
9. H                  1.51812644    2.02681493   -0.32424024
   
10. C                 -0.42910744   -1.21052830   -0.33465934
    
11. O                 -1.07773001   -2.28888229   -0.33812694
    
12. N                 -0.96370843   -0.11703076   -0.90600943
    
13. N                  0.78578354   -1.17302438    0.24032710
    
14. H                  1.16155784   -1.99994674    0.65864886
    
15. N                 -0.88665177    2.23323443   -1.54286969
    
16. H                 -1.88298814    2.18859513   -1.46992352
    
17. H                 -0.55480592    3.05391113   -1.07771587
    
18. *
    

复制代码

假如叫做task.gjf，接下来执行：

1. python3 kirariNto.py task.gjf

复制代码

程序会自动在当前目录创建一个JOBS目录，构型优化每一步的文件都会保存在其中。运行过程中，会不断输出类似如下的信息，可以看到每个态与前一个态之间的NTO重叠积分，以及最终选择了哪个态来走下一步：

1. Step    0 : Gradient = 9.646e-02/1.700e-01 (rms/max) Energy = -389.6216215640
   
2. Hessian Eigenvalues: 2.30000e-02 2.30000e-02 2.30000e-02 ... 5.66572e-01 5.76846e-01 8.02088e-01
   
3. NTO overlap between the state 1 and the previous TES is 0.10911
   
4. NTO overlap between the state 2 and the previous TES is 0.08976
   
5. NTO overlap between the state 3 and the previous TES is 0.10217
   
6. Now the state with largest overlap is determined to be 1
   

复制代码

以上述gjf文件为例，走10步后到达收敛（收敛限与Gaussian默认收敛限相同）：

1. Step   10 : Displace = 6.239e-04/1.256e-03 (rms/max) Trust = 1.141e-02 (=) Grad = 6.950e-05/1.542e-04 (rms/max) E (change) = -389.6852923380 (-8.480e-07) Quality = 1.256
   
2. Hessian Eigenvalues: 1.50842e-02 2.29339e-02 2.30019e-02 ... 6.85113e-01 7.31987e-01 7.87066e-01
   
3. Converged! =D
   

复制代码

程序地址：
https://github.com/RimoAccelerator/KirariNto/tree/main

qczgzly

看起来KirariNto在激发态优化过程中也需要计算受力，得到的是能量最小的激发态结构。那么维持激发态成分不变与不做这种考虑，所得结果一样吗，KirariNto方式优势是什么呢？

Accelerator

qczgzly 发表于 2024-1-13 18:18
看起来KirariNto在激发态优化过程中也需要计算受力，得到的是能量最小的激发态结构。那么维持激发态成分不 ...

用来处理优化过程中激发态发生交叉的情况。原文作者举了一个例子，胞嘧啶的S0最优结构垂直激发后有一个pi-pi*态和一个n-pi*态，分别是S1和S2，但S1的全局最优结构中n-pi*态变成了S1，假如提前不知道这件事，想要研究pi-pi*态，就会首先沮丧地发现收敛到了一个n-pi*态，然后再转而去优化S2。更复杂的体系，可能还会发现电子态的顺序来回反复改变，根本不知道该优化S几。使用KirariNto，就可以最大化一次优化到想研究的电子态的概率。以胞嘧啶为例，在记录初始构型的输入文件中分别写root=1和root=2，优化过程中就会自动探测并最终收敛到S2和S1势能面上的最优结构。

qczgzly

Accelerator 发表于 2024-1-13 18:36
用来处理优化过程中激发态发生交叉的情况。原文作者举了一个例子，胞嘧啶的S0最优结构垂直激发后有一个pi ...

原来如此！非常感谢您的解答~

wzkchem5

我感觉跟踪NTO可能不如跟踪跃迁密度矩阵更robust。
考虑一个四能级体系，有两个近简并的占据轨道i、j，两个近简并的空轨道a、b，其他轨道的贡献可以忽略。那么它可能会有四个激发态，(1/sqrt(2))(|i→a>+|j→b>)、(1/sqrt(2))(|i→a>-|j→b>)、(1/sqrt(2))(|i→b>+|j→a>)、(1/sqrt(2))(|i→b>-|j→a>)。那么问题来了，这四个激发态的NTO，都是由i、j按任意角度混合而成的两个占据轨道，外加由a、b按任意角度混合成的两个虚轨道，而且四个轨道的本征值都是0.5。也就是这四个激发态的NTO的区别，完全取决于对角化时引入的简并本征态之间的随机酉变换，结构优化时前一步的某一个态被指认成后一步的四个激发态里的任意一个的概率都是相等的。
但是另外一方面，这四个态又是货真价实的四个不同的态，它们的跃迁密度矩阵X是完全不同的，分别是1/sqrt(2)乘以以下四个矩阵

1. 1 0
   
2. 0 1
   
3. 
   
4. 1 0
   
5. 0 -1
   
6. 
   
7. 0 1
   
8. 1 0
   
9. 
   
10. 0 1
    
11. -1 0

复制代码

只不过计算NTO时对角化的是X*X^T、X^T*X，而这四个态的X*X^T、X^T*X矩阵是一样的，都是

1. 1/2 0
   
2. 0 1/2

复制代码

所以也就是说，存在靠NTO无法区分，但靠跃迁密度矩阵能区分的情况。另一方面，如果两个激发态的NTO不一样，或NTO的本征值不一样，那么就意味着X*X^T不一样，因此X肯定不一样。所以不存在NTO能区分、但跃迁密度矩阵区分不了的情形。
不仅如此，以上的例子并不是凭空编造的，而是有现实原型，也就是卟啉配合物。卟啉的HOMO和HOMO-1接近简并，LUMO和LUMO+1严格简并，因此确实会出现以上的现象，最低的4个单重态激发态分裂成两组，2eV左右有一组简并激发态，叫做Q带；3eV左右有另一组简并激发态，叫做B带或Soret带（https://pubs.aip.org/aip/jcp/art ... Substitution-on-the）。虽然真实的卟啉配合物的HOMO和HOMO-1不会严格简并，但是不好说有没有某种取代卟啉的轨道简并性足够好，导致激发态跟踪错了。
考虑到直接计算两个跃迁密度矩阵的overlap的计算量可能比较大，实际操作可以把跃迁密度矩阵做奇异值分解（SVD），易知得到的奇异矢量就是NTO，但相比一般计算NTO的步骤，由跃迁密度矩阵SVD得到NTO有两个优点：（1）即使有两个NTO的本征值偶然简并，也能分清楚哪个占据NTO跃迁到了哪个虚NTO，也就是能分清【50%的i→a，50%的j→b】和【50%的i→b，50%的j→a】这两种情况；（2）由奇异值的正负号，可以知道比如说这50%的i→a和50%的j→b之间是同相位还是反相位。所以如果手头已经有一个基于NTO做激发态跟踪的程序，那么程序核心部分可以不变，把NTO的计算改为由跃迁密度矩阵的SVD得到，在比较两个结构下的NTO时把占据NTO/虚NTO的对应关系，以及奇异值也考虑上，这样应该能得到更robust的结果。

zjxitcc

推荐一篇很相关的文章《Transition orbital projection approach for excited state tracking》https://doi.org/10.1063/5.0081207
里面也有介绍一些仅凭NTO无法分辨的例子。当初文章作者想算法时问了wzkchem5一些这方面的问题。

Accelerator

wzkchem5 发表于 2024-1-13 22:28
我感觉跟踪NTO可能不如跟踪跃迁密度矩阵更robust。
考虑一个四能级体系，有两个近简并的占据轨道i、j，两 ...

我对此有一个问题。
对于两两简并的四能级系统，(1/sqrt(2))(|i→a>+|j→b>)和(1/sqrt(2))(|i→a>-|j→b>)能量不同非常好理解（与此同时它们对应的Multiwfn中定义的电子/空穴密度也是不同的，因为定义式中的交叉项不同），但(1/sqrt(2))(|i→a>+|j→b>)和(1/sqrt(2))(|i→b>+|j→a>)两个态（在楼上推荐的文章中叫做bijectivity导致的两个态），是否有可能在大多数情况下（特别是因为对高称性导致轨道简并的情况下）是简并的？目前是否有已知的完全由于bijectivity导致激发态能量不同的例子？
上文给出了C18环的S2和S12态作为此案例，但这两个态对应的NTO本征值不同，不能作为严格的反例。

wzkchem5

Accelerator 发表于 2024-1-14 14:28
我对此有一个问题。
对于两两简并的四能级系统，(1/sqrt(2))(|i→a>+|j→b>)和(1/sqrt(2))(|i→a>-|j→b ...

有的，比如i和a局域在一个片段上，j和b局域在另一个片段上，两个片段相距很远。那么(1/sqrt(2))(|i→a>+|j→b>)的激子结合能很负，(1/sqrt(2))(|i→b>+|j→a>)的激子结合能接近0，所以前者能量更低。