SCF convergence error in the determination of the excited state energy

Makhvela

Hello, I faced SCF convergence error in the determination of the excited state energy. How can I solve this convergence error? thank you

sobereva

Please upload output file so that current convergence status can be understood by others.

Makhvela

The output file is heavy. I am trying to upload it, but it didn't upload.

zjxitcc

SF-TDA, MRSF-TDA, SA-SF-TDA methods are usually based on the ROKS (restricted open-shell Kohn-Sham, also called RODFT) reference wave function. But ROKS is often very difficult to be converged. The SCF algorithm in GAMESS is not as robust as that in mainstream/popular quantum chemistry (QC) programs, which makes the ROKS calculation even harder.

I wrote a tutorial (in Chinese) which aims to solve this problem thoroughly. The idea is to separate the target calculation into two steps: (i) ROKS calculation, (ii) SF-/MRSF-/SA-SF-/TDA calculation. We can use one QC program (which has very robust ROKS code, e.g. PySCF/ORCA/Gaussian) to perform the ROKS calculation, and use another program (GAMESS/OpenQP/Q-Chem, etc) to perform the desired spin-flip calculation. A combination of two programs makes the whole calculation more efficient and fewer SCF convergence faliures.

I'll take your molecule as an example. Firstly, we write/create the following Python script `Cz2BP_MRSFt_PBE0.py`

1. from pyscf import gto,  dft, lib
   
2. from mokit.lib.py2gms import py2gms
   
3. 
   
4. lib.num_threads(64) # 64 CPU cores
   
5. mol = gto.M()
   
6. # 82 atom(s)
   
7. mol.atom = '''
   
8. C               3.76399300   -0.50818100    0.38545400
   
9. C              2.60277900    0.14659700   -0.06114300
   
10. (not shown)
    
11. '''
    
12. mol.basis = {
    
13. (def-TZVP basis set data, not shown)
    
14. }
    
15. mol.charge = 0
    
16. mol.spin = 2
    
17. mol.verbose = 4
    
18. mol.build(parse_arg=False)
    
19. 
    
20. # routine ROKS calculation
    
21. mf = dft.ROKS(mol)
    
22. mf.xc = 'pbe0'
    
23. mf.grids.atom_grid = (99,590)
    
24. mf.verbose = 4
    
25. mf.max_cycle = 50
    
26. mf.max_memory = 160000 # 160 GB
    
27. mf.level_shift = 0.1
    
28. mf.kernel()
    
29. if mf.converged is False:
    
30.   mf = mf.newton()
    
31.   mf.kernel()
    
32. 
    
33. # stable=opt for ROKS wave function
    
34. mo, _, stable, _ = mf.stability(tol=1e-7,return_status=True)
    
35. if stable is False:
    
36.   dm = mf.make_rdm1(mo, mf.mo_occ)
    
37.   mf.kernel(dm0=dm)
    
38.   mo, _, stable, _ = mf.stability(tol=1e-7,return_status=True)
    
39. 
    
40. # export GAMESS .inp file
    
41. py2gms(mf, 'Cz2BP_MRSFt_PBE0.inp', mrsf=True)
    

复制代码

This is also an input file of the open source QC package PySCF. You can find this script attached Cz2BP_MRSFt_PBE0.py (7.83 KB, 下载次数 Times of downloads: 0)

2025-11-11 13:54 上传 Uploaded

点击下载
Click to download

You can submit the PySCF job, e.g.

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

复制代码

This python script says that firstly the ROKS PBE0 calculation would be performed, secondly check the ROKS wave function stability and try to find a lower solution (in an iterative way) if the ROKS wave function is not stable. Thirdly, export/punch stable ROKS wave function and create GAMESS MRSF-TDA input file called `Cz2BP_MRSFt_PBE0.inp`. Here py2gms is a function from MOKIT, which is able to (correctly) transform MOs from PySCF to GAMESS. You can find Cartesian coordinates, basis set data (def-TZVP), MRSF keywords, MO coefficients are well written in .inp

1. $CONTRL SCFTYP=ROHF RUNTYP=ENERGY ICHARG=0 MULT=3 NOSYM=1 ICUT=11
   
2.   MAXIT=200 ISPHER=1 DFTTYP=PBE0 TDDFT=MRSF $END
   
3. $SYSTEM MWORDS=300 $END
   
4. $SCF DIRSCF=.T. FDIFF=.F. MOM=.T. $END
   
5. $GUESS GUESS=MOREAD NORB=1090 $END
   
6. $DFT NRAD0=99 NLEB0=590 NRAD=99 NLEB=590 $END
   
7. $TDDFT NSTATE=5 NRAD=99 NLEB=590 $END
   

复制代码

Note that `DFTTYP=BHHLYP` is written in the .inp file by default since the BHandHLYP is the most frequently used functional in spin-flip calculations. Since you use PBE0, we need to modify it into `DFTTYP=PBE0`. Usually the SCF iteration in GAMESS would be converged in 1~2 cycles when we apply this trick. And the program immediately moves forward to the MRSF calculation. There are some rare cases that the energy goes higher or oscillating even if the initial ROKS energy is already the lowest. We can add the keyword `MOM=.T.` if such case happens. Finally we can submit the GAMESS job, e.g.

1. rungms Cz2BP_MRSFt_PBE0.inp 00 64 >Cz2BP_MRSFt_PBE0.gms 2>&1 &
   

复制代码

A smooth convergence can be found in GAMESS output file (.gms)

1.   ITER EX      TOTAL ENERGY        E CHANGE  DENSITY CHANGE    DIIS ERROR      INTEGRALS    SKIPPED
   
2.           * * *   INITIATING DIIS PROCEDURE   * * *
   
3.    1  0    -1843.2831356194 -1843.2831356194   0.044601763   0.000078659    2194165061210776939510
   
4.    2  1    -1843.2831355007     0.0000001186   0.001964492   0.000009821    2193945982510777418946
   
5.    3  2    -1843.2831193900     0.0000161107   0.001395311   0.000574187    2193943440710777425512
   
6.    4  3    -1843.2831354921    -0.0000161021   0.000055589   0.000011031    2193944894210777422167
   
7.    5  4    -1843.2831354991    -0.0000000071   0.000038054   0.000006210    2193945115510777421894
   
8.    6  5    -1843.2831355005    -0.0000000014   0.000043029   0.000004214    2193944894310777422381
   
9.    7  6    -1843.2831355012    -0.0000000007   0.000013974   0.000001654    2193944769910777422655
   
10.    8  7    -1843.2831355008     0.0000000004   0.000010141   0.000003306    2193944743110777422770
    
11.    9  8    -1843.2831355012    -0.0000000004   0.000014131   0.000001722    2193944703110777422832
    
12.   10  9    -1843.2831355012    -0.0000000001   0.000013286   0.000001411    2193944612910777422949
    
13.   11 10    -1843.2831355013    -0.0000000001   0.000004101   0.000000793    2193944530810777423034
    
14.   12 11    -1843.2831355013     0.0000000000   0.000007330   0.000000671    2193944506410777423072
    
15.   13 12    -1843.2831355013    -0.0000000000   0.000001489   0.000000591    2193944470610777423118
    
16.   14 13    -1843.2831355013     0.0000000000   0.000005537   0.000000944    2193944503710777423108
    
17.   15 14    -1843.2831355012     0.0000000001   0.000006524   0.000001680    2193944509210777423084
    
18.   16 15    -1843.2831355013    -0.0000000001   0.000000745   0.000000180    2193944511910777423067
    
19.   17 16    -1843.2831355013    -0.0000000000   0.000000127   0.000000014    2193944505710777423069
    
20. 
    
21.           -----------------
    
22.           DENSITY CONVERGED
    
23.           -----------------
    

复制代码

Many people who want to perform spin-flip calculations encountered such problem and solve it successfully by using this trick. Maybe you can, or you have solved this problem by using other techniques, but I want to remind you that this molecule has two different ROKS solutions

1. converged SCF energy = -1843.27474917186
   
2. converged SCF energy = -1843.28312805381
   

复制代码

This 1st solution is energetically higher and not stable (which means it is a saddle point in the wave function parameter space), while the 2nd solution is stable and energetically favorable.

By the way, if you did not use MOKIT and PySCF before, but now are williing to have a try, you can simply run

1. conda create -n mokit-py39 python=3.9
   
2. conda activate mokit-py39
   
3. conda install mokit -c mokit
   
4. pip install pyscf
   
5. conda deactivate
   

复制代码

to install both MOKIT and PySCF. For further installation guide, please click here. For more utilities/modules to transform MOs among QC packages, you can read the following tutorials
《利用MOKIT从PySCF向其他量化程序传轨道》，《利用MOKIT从Gaussian向其他量化程序传轨道》，《利用MOKIT从ORCA向其他量化程序传轨道》