Spin-flip方法中RODFT难收敛解决办法

zjxitcc

最常见的Spin-flip方法为SF-CIS和SF-TDDFT，它们一般以高自旋三重态为参考态波函数，将一个alpha电子翻转为beta电子，产生满足<Sz>=0的大量行列式，随后构造CI哈密顿矩阵进行对角化获得单重态、三重态、五重态、...，每个电子态都由多个行列式线性组合而成。


(copied from ORCA 5.0.4 manual P794)

由于这类方法通常可以正确描述具有多组态特征的S0基态、正确描述S0/S1圆锥交叉点，且（部分方法）有解析梯度可进行结构优化、寻找过渡态和极小能量交叉点等任务，是多参考方法外的最优选择。在平时的计算中高自旋三重态可以用ROHF/RODFT算，也可以用UHF/UDFT算；但在Spin-flip这类方法中，为尽可能减小自旋污染或严格保持自旋纯态，经常只能选择ROHF/RODFT。因为我们的目的是计算多个电子态，自旋污染会导致我们无法根据自旋多重度辨认各个电子态。如果连S0, S1, T1, T2都说不清楚，那研究光化学/光物理/双自由基等问题只能是一堆浆糊。本文中RODFT也称为ROKS，UDFT也称为UKS。

在五、六年甚至更早以前，如果做ROHF/RODFT计算碰到SCF不收敛，基本上没有特别高效的解决手段，只能是尝试一个个关键词，且通常无法检验波函数稳定性，很难判断侥幸收敛的开壳层波函数是否靠谱。这几年各个量化软件的发展让笔者觉得这个问题已经有系统性的解决办法，且成功率在99%以上，不需要盲目尝试，因此有必要单独写一篇分享一下计算技巧。

1. RODFT计算
RODFT计算与Spin-flip方法的计算可以分别使用不同的量化程序，不是非得在一个程序中完成。Gaussian软件目前仍然是SCF收敛性最稳健的程序，我们可以用它做完RODFT计算，然后传轨道给其他程序进行后续的Spin-flip方法计算。这里展示一个基态O2的单点计算文件

1. %chk=O2_T.chk
   
2. %mem=4GB
   
3. %nprocshared=2
   
4. #p ROBHandHLYP/def2TZVP nosymm int=nobasistransform
   
5. 
   
6. title
   
7. 
   
8. 0 3
   
9. O   0.0   0.0   0.0
   
10. O   0.0   0.0   1.1616
    

复制代码

这个例子SCF收敛没有什么难度，仅作为示例方便阅读。做完Gaussian计算后有O2_T.chk文件，运行

1. formchk O2_T.chk O2_T.fch

复制代码

产生O2_T.fch文件。

开源程序PySCF支持ROHF/RODFT的轨道优化二阶算法，且支持波函数稳定性检验，是目前几乎最强大的RO方法计算程序。这里的描述为笔者个人经过大量计算经验总结得出，不代表任何读者或软件官方的观点。可能很多读者不熟悉用PySCF做这种计算，这里展示一个详细的输入文件脚本

1. from pyscf import gto, dft, lib
   
2. from mokit.lib.py2fch_direct import fchk
   
3. 
   
4. lib.num_threads(64)
   
5. mol = gto.M()
   
6. mol.atom = '''
   
7. C   -0.12264859        0.00599759       -0.75931820
   
8. C   -0.13393227        0.05086047        0.75572332
   
9. O    1.29017598        0.03740604       -0.91156021
   
10. O    1.12879647       -0.57344495        0.94545004
    
11. H   -0.54539878       -0.92543383       -1.16182063
    
12. H   -0.91791476       -0.55372886        1.23567184
    
13. H   -0.13429527        1.07605238        1.15229831
    
14. H   -0.56478278        0.88229116       -1.25644448
    
15. '''
    
16. 
    
17. mol.basis = 'aug-cc-pVDZ'
    
18. mol.charge = 0
    
19. mol.spin = 2
    
20. mol.verbose = 4
    
21. mol.build(parse_arg=False)
    
22. 
    
23. mf = dft.ROKS(mol)
    
24. mf.xc = 'bhandhlyp'
    
25. mf.grids.atom_grid = (99,590)
    
26. mf.max_cycle = 128
    
27. mf.max_memory = 128000 #MB
    
28. old_e = mf.kernel()
    
29. 
    
30. mo = mf.stability()[0]
    
31. dm = mf.make_rdm1(mo, mf.mo_occ)
    
32. mf.kernel(dm0=dm)
    
33. 
    
34. mf.stability()
    
35. fchk(mf, 'high_spin.fch', density=True)

复制代码

文件含义：先进行正常RODFT计算，SCF收敛后电子能量为-228.760185 a.u.，检验波函数稳定性，发现不稳定，mf.stability()会返回新的轨道，使用该轨道再次进行SCF计算，收敛后能量为-228.787950 a.u.，再检验波函数稳定性，发现稳定，导出RODFT轨道进fch文件。fchk是MOKIT的一个模块，用于在PySCF计算结束后导出fch文件。

若读者尚未进行任何Spin-flip计算，或计算中途遇到RODFT SCF不收敛，建议采用上述任一程序进行计算。另外，在平时的DFT计算中不同情景下常用泛函不同，但在Spin-flip方法中最常用的泛函是杂化泛函BHHLYP（Gaussian输入文件中写作BHandHLYP），不宜使用纯泛函。

2. RODFT波函数稳定性
对于开壳层体系的SCF结果，笔者建议检验一下波函数稳定性（除非是你已经算过很多次、非常了解的体系）然后再进行Spin-flip方法的计算。上面已经展示过如何使用PySCF进行RODFT计算并检验波函数稳定性，但读者碰到的情况可能更多是手上已有一个波函数文件，但未检验稳定性。假设文件名叫high_spin.fch，我们可以写一个Python脚本stab.py，内容如下

1. from pyscf import dft, lib
   
2. from mokit.lib.gaussian import load_mol_from_fch
   
3. from mokit.lib.fch2py import fch2py
   
4. from mokit.lib.rwwfn import read_nbf_and_nif_from_fch, read_na_and_nb_from_fch, \
   
5. get_occ_from_na_nb
   
6. 
   
7. lib.num_threads(64) # 64 CPU cores
   
8. fchname = 'high_spin.fch'
   
9. 
   
10. mol = load_mol_from_fch(fchname)
    
11. mf = dft.ROKS(mol)
    
12. mf.xc = 'bhandhlyp'
    
13. mf.grids.atom_grid = (99,590)
    
14. mf.max_memory = 128000 #MB
    
15. nbf, nif = read_nbf_and_nif_from_fch(fchname)
    
16. na, nb = read_na_and_nb_from_fch(fchname)
    
17. mo = fch2py(fchname, nbf, nif, 'a')
    
18. mo_occ = get_occ_from_na_nb(nif, na, nb)
    
19. dm = mf.make_rdm1(mo, mo_occ)
    
20. mf.kernel(dm0=dm)
    
21. mf.stability()

复制代码

提交PySCF任务

1. python stab.py >stab.out 2>&1

复制代码

在输出文件中可以看到SCF 2圈收敛，且波函数是稳定的

1. init E= -228.787949822635
   
2.   HOMO = -0.205358439066859  LUMO = 0.00469108326441472
   
3. cycle= 1 E= -228.787949757274  delta_E= 6.54e-08  |g|= 3.9e-06  |ddm|= 5.78e-06
   
4.   HOMO = -0.205356878045115  LUMO = 0.00469095686765871
   
5. cycle= 2 E= -228.78794975727  delta_E= 3.98e-12  |g|= 5.37e-06  |ddm|= 6.31e-06
   
6.   HOMO = -0.20535887115258  LUMO = 0.0046911310210829
   
7. Extra cycle  E= -228.787949757264  delta_E= 6.59e-12  |g|= 7.77e-06  |ddm|= 9.14e-06
   
8. converged SCF energy = -228.787949757264
   
9. <class 'pyscf.dft.roks.ROKS'> wavefunction is stable in the internal stability analysis

复制代码

若读者的RODFT是Q-Chem算的，Q-Chem目前不支持RODFT波函数稳定性检验。不过Q-Chem算完就有fch文件（根据版本的不同，有可能后缀是.FChk），启动Python，运行

1. from mokit.lib.qchem import standardize_fch
   
2. standardize_fch('high_spin.fch')

复制代码

可获得文件high_spin_std.fch，后面的步骤就与本小节一开始写的步骤一样了。

ORCA的SF-TDDFT目前只能基于三重态UHF/UDFT参考态，不能基于ROHF/RODFT，因此上文所示技巧没有直接帮助。若好奇仍想检验RODFT波函数稳定性（到ORCA 5.0.4为止仍不支持），可以运行

1. orca_2mkl high_spin -molden
   
2. molden2fch high_spin.molden -orca

复制代码

产生high_spin.fch文件，检验波函数稳定性的步骤不再赘述。

3. Spin-flip类型计算
在获得fch文件后便可直接产生各种Spin-flip计算的输入文件，同时传收敛的DFT轨道，具体示例见下。

3.1 GAMESS的SF-TDDFT
运行

1. fch2inp high_spin.fch -sf

复制代码

产生high_spin.inp文件，坐标、基组信息、轨道数据和方法关键词全都写好了，可直接提交SF-TDDFT计算任务

1. $GMS high_spin.inp 00 16 >high_spin.gms 2>&1 &

复制代码

其中16表示16核并行。输出文件中可以看到SCF立即收敛

1. ITER EX      TOTAL ENERGY        E CHANGE  DENSITY CHANGE    DIIS ERROR      INTEGRALS    SKIPPED
   
2.           * * *   INITIATING DIIS PROCEDURE   * * *
   
3.    1  0     -228.7879496773  -228.7879496773   0.005965249   0.000148515       43190509       1148
   
4.    2  1     -228.7879496421     0.0000000352   0.000540755   0.000005344       42883334      14932
   
5.    3  2     -228.7879496401     0.0000000020   0.000095039   0.000001233       42421386      33461
   
6.    4  3     -228.7879496398     0.0000000002   0.000029150   0.000001298       41913822      51524
   
7.    5  4     -228.7879496398     0.0000000001   0.000001886   0.000000547       41554248      65582
   
8.    6  5     -228.7879496398    -0.0000000000   0.000007179   0.000000066       40613185     101140
   
9.    7  6     -228.7879496398    -0.0000000000   0.000000597   0.000000012       40398109     106392
   
10.    8  7     -228.7879496398    -0.0000000000   0.000000234   0.000000015       38760838     164687
    
11. 
    
12.           -----------------
    
13.           DENSITY CONVERGED
    
14.           -----------------

复制代码

立即进入SF-TDDFT计算环节，默认计算5个电子态。GAMESS的SCF收敛性比较差，我们通过传轨道让它迅速完成SCF计算，可以节约很多人力物力。这里的文件名high_spin.fch仅做为示例，表示三重态（或更高）自旋多重度的一个波函数fch文件，读者可使用任意文件名。这里也可以提供三重态UHF/UDFT的fch文件，但会增大目标电子态的自旋污染，一般不建议使用。

3.2 GAMESS的MRSF-TDDFT
运行

1. fch2inp O2_T.fch -mrsf

复制代码

产生O2_T.inp文件，可直接提交MRSF-TDDFT计算任务，默认计算5个单重态。注意MRSF只能基于限制性开壳层，勿提供UHF/UDFT的fch文件。由于MRSF几乎无自旋污染，且已经实现了解析导数，因此可以取代SF-TDDFT。

3.3 ORCA的SF-TDDFT
运行

1. fch2mkl O2_T.fch -sf

复制代码

产生ORCA输入文件O2_T_o.inp和轨道文件O2_T_o.mkl。将mkl转化为gbw文件

1. orca_2mkl O2_T_o -gbw

复制代码

即可提交SF-TDDFT计算任务。inp文件中的核数和内存可以根据自己的需求修改。上文提到过ORCA不支持从RODFT出发做SF计算，请注意提供三重态UDFT的fch文件给fch2mkl小程序。

3.4 Q-Chem的SA-SF-DFT
运行

1. fch2qchem high_spin.fch -sasf

复制代码

即可提交SA-SF-DFT计算任务。假设我们用16线程并行

1. qchem -nt 16 -np 1 high_spin.in high_spin.out high_spin

复制代码

注意SA-SF-DFT是自旋纯态方法，只能从RODFT出发，勿使用UDFT的fch文件。

在这4个例子中，fch2inp，fch2mkl，fch2qchem是MOKIT的小程序，分别用于传轨道给GAMESS/ORCA/Q-Chem软件，并且已经帮用户写好了体系信息和计算关键词。使用这些技巧完成计算，发表文章时记得需要引用RODFT计算程序、Spin-flip方法计算程序和MOKIT。

4. 支持Spin-flip方法的程序
根据笔者的经验，目前Spin-flip方法总结如下（这里只讨论可获取的免费程序或可购买的商业程序，不讨论不开放的程序）

程序	方法	解析导数	从何种自旋出发	指定激发态自旋	结果自旋纯态
GAMESS	SF-TDDFT	有	三重态或更高	否	否
ORCA	SF-TDDFT	有	仅三重态	否	否
GAMESS	MRSF-TDDFT	有	仅三重态	支持	几乎纯态，无法看出自旋污染
Q-Chem	SA-SF-TDDFT	暂无	三重态或更高	暂无	自旋纯态

第4列"从何种自旋出发"的意思是程序允许从哪一种自旋多重度的RODFT出发做spin-flip计算，例如从三重态RODFT出发翻转1个电子可以得到单重态和三重态（还有五重态等，需求较少），从四重态出发翻转1个电子可以得到二重态和四重态，也就是说GAMESS的SF-TDDFT理论上是可以寻找二重态-二重态极小能量交叉点和二重态-四重态极小能量交叉点的。Q-Chem的SA-SF-DFT目前没有解析梯度，用于结构优化和寻找交叉点还是太费劲了，一般只能用来做单点计算获得自旋纯态的各个电子态。

第5列"指定激发态自旋"意思是，是否可以计算某个特定自旋多重度的多个电子态。这很容易理解，在传统TDDFT计算中也有类似的需求：若参考态波函数是闭壳层RDFT，则可通过TD(singlet)和TD(triplet)关键词来指定计算单重态激发态或三重态激发态。在这里也是类似的意思。GAMESS的MRSF-TDDFT支持指定单重态、三重态甚至五重态；而不论哪个程序的SF-TDDFT都不能指定目标电子态的自旋，只能是算出什么结果就接受什么。SA-SF-DFT目前无法同时计算单重态和三重态：假设从三重态出发，只能获得多个单重态，没法获得任何三重态的电子态；从四重态出发，只能获得多个二重态，没法给出任何四重态的电子态。

目前这些方法都没有实现计算旋轨耦合的功能，将来只有在SA-SF-DFT中有希望实现，因为其他方法都不满足自旋纯态。可以看到，Spin-flip方法还有巨大的发展空间：解析二阶导数、非绝热耦合矢量、指定目标自旋子空间、虚轨道空间压缩、态追踪算法、超越TDA近似等许多方向都急需高效的程序实现。不过，即使是自旋纯态的Spin-flip方法也不是十全十美的，也存在理论局限性，并且同时还有自旋纯态RO-TDDFT和X-TDDFT等其他开壳层TD方法参与竞争。
2024年5月1日注：MRSF-TDDFT开发者在2023年发表了计算旋轨耦合的文章，估计不久后可以在GAMESS中用到。由于MRSF-TDDFT不严格满足自旋纯态，在计算SOC时会扔掉一小部分引起自旋污染的组态，用剩下的纯态进行SOC计算。

免责声明：笔者不是从事Spin-flip方法开发的专家，上文若有错误或争议，欢迎留言讨论。

相关阅读
1. 离线安装PySCF-2.x
2. 利用MOKIT从Gaussian向其他量化程序传轨道
3. 利用MOKIT从PySCF向其他量化程序传轨道