Calculating singlet-triplet gap and open-shell singlet character

Question

1. I am trying to calculate the singlet-triplet energy gap ($S_0-T_1$) of a conjugated system (18 atoms). For the singlet $S_0$ state, assumed to be closed shell, the Gaussian options were

   opt=tight b3lyp/6-311g(d,p) integral(ultrafinegrid)
   ...
   0 1
   ...
   

   and for the energy of the triplet $T_1$ state

   opt=tight rob3lyp/6-311g(d,p) integral(ultrafinegrid)
   ...
   0 3
   ...
   

   Does the energy difference of those two calculations have any physical significance? (i.e. should have I used ...td-ub3lyp...0 3 ... instead?)

2. I believe a system I am currently working on may have some open-shell character for the singlet ground state so I performed the following with the optimized geometry I got from a TD calculation of a triplet state:

   #p ub3lyp/6-311+g(d,p) scf=(tight) guess=(mix,always) gfinput gfprint iop(6/7=3)
   ...
   0 1
   ...
   

   For the stability check, do I have to use the guess=(mix,always) on they keyword route as well?

Answer 1

There is a related discussion on which method to use for a singlet triplet guess, which I encourage you to also read: U- or RO-method for Singlet-Triplet Gap?

This post is going to be a little longer, so I'll present the key points upfront. It is reasonable to calculate (and/or optimise) the wave function with a restricted method, then test it for instabilities. The same applies to to triplet (etc.) wave functions, where I would also try to start with a restricted open shell approach. To answer 1, it is a reasonable first approach for a singlet-triplet gap. If you are more interested in the electronic excitations a time-dependent treatment is necessary.

As a general recommendation I would advise to split the calculations up and read the respective wave functions from the checkpoint file to perform additional calculations and checks, so that you always have a fallback.
Hence the answer to 2, once you have created and converged a solution to the wave function you can read it in again. Do not generate a new initial guess.

You should also always keep an eye on spin contamination, because this is essentially telling you if the method you are using is suitable in the first place. When it comes to broken symmetry solutions, you should always consider to run multi-reference calculations, too.

(Note that guess(always) applies only to geometry optimisations; see below.)

---

I'm going to exercise a couple of calculations on dioxygen, because it's a nice study case. I will be using the BP86 functional with the def2-SVP basis set and density fitting. I have done the calculations with Gaussian 09 rev. D.01.
It is somewhat an expansion on the example in Exploring Chemistry with Electronic Structure Methods, 3rd Edition (p.71).

I am going to present the full input of the calculations and at the end the absolute energy to compare to, in case you'd like to run the calculations yourself. Please note that I have not run any geometry optimisations and assumed a equilibrium bond distance of 120 pm. I think that is a reasonable choice since this is just an example and the optimisations can be considered as additional exercise.

I labelled the the used methods RHF, ROHF, and UHF even though I am using density functional approximations.

Because the main focus is the use of the stable keyword, we start with a wave function that is already know to be not stable. The RHF solution of dioxygen, filename bp86svp.rhf.sing.gjf:

%chk=bp86svp.rhf.sing.chk
%nproc=2
%mem=10000MB
#p RBP86/Def2SVP/W06
DenFit
gfinput gfoldprint iop(6/7=3)
symmetry(loose)
pop=full pop=nbo
scf(xqc)

title card unused

0 1
O  0.6 0.0 0.0
O -0.6 0.0 0.0

This input uses a few keywords, that are not strictly necessary to complete the calculations. For example the route section #p requests detailed output. DenFit is already implied with the choice of W06 as an additional basis set. I like to write it in nevertheless, to exactly know what I am doing. The same goes for symmetry(loose); occasionally I like to turn symmetry of, so I know this time I really wanted it. The line gfinput gfoldprint iop(6/7=3) is included so that molden (❤) can read (and display) the orbitals. There is no shame in always requesting population analyses.

In this very case even scf(MaxCycle=9999) will not lead to convergence; with scf(xqc) we are allowing the algorithm to switch to quadratic convergence if conventionally the convergence cannot be achieved. Note that this may converge to a different solution than originally attempted.

For publication type calculations I would also include int(ultrafine) and scf(tight); but I forgot it here.

From here I went on to calculate the stability of the wave function, i.e. bp86svp.rhf.sing.stab.gjf:

%chk=bp86svp.rhf.sing.stab.chk
%nproc=2
%mem=10000MB
%oldchk=bp86svp.rhf.sing.chk
#p RBP86/Def2SVP/W06
DenFit
geom=check guess=read
gfinput gfoldprint iop(6/7=3)
scf(xqc)
stable

title card unused

0 1

Note that in this case the SCF will converge within the conventional cycles and hence the energy is different. You can additionally check this by reading in the checkpoint file and perform a new SCF cycle.

However, in the stability analysis you will find the following output:

 ***********************************************************************
 Stability analysis using  singles matrix:
 ***********************************************************************

 Eigenvectors of the stability matrix:

 Excited state symmetry could not be determined.
 Eigenvector   1:      Triplet-?Sym  Eigenvalue=-0.0560505  =2.000
       8 ->  9         0.70645

 Excited state symmetry could not be determined.
 Eigenvector   2:      Singlet-?Sym  Eigenvalue= 0.0372059  =0.000
       8 ->  9         0.70689

 Excited state symmetry could not be determined.
 Eigenvector   3:      Triplet-?Sym  Eigenvalue= 0.1233031  =2.000
       7 ->  9         0.70493

 Excited state symmetry could not be determined.
 Eigenvector   4:      Triplet-?Sym  Eigenvalue= 0.2035093  =2.000
       5 ->  9         0.70373

 Excited state symmetry could not be determined.
 Eigenvector   5:      Singlet-?Sym  Eigenvalue= 0.2902800  =0.000
       6 ->  9         0.70685

 Excited state symmetry could not be determined.
 Eigenvector   6:      Singlet-?Sym  Eigenvalue= 0.3353297  =0.000
       5 ->  9         0.70254
 SavETr:  write IOETrn=   770 NScale= 10 NData=  16 NLR=1 NState=    6 LETran=     118.
 The wavefunction has an RHF -> UHF instability.

As expected we see a couple of instabilities. At this point we can run another calculation with stable=opt to find such a wave function. For reasons I don't quite understand why it finds the the broken symmetry solution instead of the triplet solution.

For exercise purposes we'll manually generate the broken symmetry solution; i.e. bp86svp.uhf.sing.gjf:

%chk=bp86svp.uhf.sing.chk
%nproc=2
%mem=10000MB
%oldchk=bp86svp.rhf.sing.chk
#p UBP86/Def2SVP/W06
DenFit
scf(xqc)
gfinput gfoldprint iop(6/7=3)
symmetry(loose)
pop=full pop=nbo
geom=check guess=(read,mix)

title card unused

0 1

The keyword guess(mix) generates a UHF guess by mixing the HOMO with the LUMO. You only need to do this once. Note that the keyword guess(always) only applies to optimisations and will perform a new initial guess at every point along the optimisation. I personally have no idea why this option exists or what it is used for.

Now we can go ahead and test for stability; i.e. bp86svp.uhf.sing.stab:

%chk=bp86svp.uhf.sing.stab.chk
%nproc=2
%mem=10000MB
%oldchk=bp86svp.uhf.sing.chk
#p UBP86/Def2SVP/W06
DenFit
scf(xqc)
gfinput gfoldprint iop(6/7=3)
symmetry(loose)
pop=full pop=nbo
geom=check guess=read
stable

title card unused

0 1

Note that we do not need to mix the orbitals in the wave function again. The stability section will now read like this:

 ***********************************************************************
 Stability analysis using  singles matrix:
 ***********************************************************************

 Eigenvectors of the stability matrix:

 Excited state symmetry could not be determined.
 Eigenvector   1:  1.012-?Sym  Eigenvalue=-0.0000016  =0.006
      8A ->  9A        0.70640
      8B ->  9B       -0.70640

 Excited state symmetry could not be determined.
 Eigenvector   2:  1.006-?Sym  Eigenvalue= 0.0930085  =0.003
      8A ->  9A        0.70688
      8B ->  9B        0.70688

 Excited state symmetry could not be determined.
 Eigenvector   3:  2.245-?Sym  Eigenvalue= 0.2031898  =1.010
      7A ->  9A       -0.70381
      7B ->  9B        0.70381
 SavETr:  write IOETrn=   770 NScale= 10 NData=  16 NLR=1 NState=    3 LETran=      64.
 The wavefunction is stable under the perturbations considered.

We already know that the triplet state should be more stable. As I hinted, I like to start with a ROHF calculation. In principle these calculations are always instable with respect to an UHF solution. However, usually this leads to faster results, especially if you have already a converged RHF solution. In this case we could readily work from the broken symmetry solution and directly calculate the triplet from there. Choose whatever is more convenient for you.
For this illustrative purpose I chose to generate a new initial guess; i.e. bp86svp.rohf.trip.com:

%chk=bp86svp.rohf.trip.chk
%nproc=2
%mem=10000MB
#p ROBP86/Def2SVP/W06
DenFit
gfinput gfoldprint iop(6/7=3)
symmetry(loose)
pop=full pop=nbo

title card unused

0 3
O  0.6 0.0 0.0
O -0.6 0.0 0.0

And get the UHF triplet from there; i.e. bp86svp.uhf.trip.gjf:

%chk=bp86svp.uhf.trip.chk
%nproc=2
%mem=10000MB
%oldchk=bp86svp.rohf.trip.chk
#p UBP86/Def2SVP/W06
DenFit
gfinput gfoldprint iop(6/7=3)
symmetry(loose)
pop=full pop=nbo
scf(xqc)
geom=check guess=read

title card unused

0 3

And then obviously check for instabilities; i.e. bp86svp.uhf.trip.stab.gjf:

%chk=bp86svp.uhf.trip.stab.chk
%nproc=2
%mem=10000MB
%oldchk=bp86svp.uhf.trip.chk
#p UBP86/Def2SVP/W06
DenFit
gfinput gfoldprint iop(6/7=3)
symmetry(loose)
pop=full pop=nbo
scf(xqc)
geom=check guess=read
stable

title card unused

0 3

We will find the following stability section:

 ***********************************************************************
 Stability analysis using  singles matrix:
 ***********************************************************************

 Eigenvectors of the stability matrix:

 Excited state symmetry could not be determined.
 Eigenvector   1:  3.003-?Sym  Eigenvalue= 0.2313384  =2.004
      6B ->  8B        0.70567
      7B ->  9B        0.70567

 Excited state symmetry could not be determined.
 Eigenvector   2:  3.003-?Sym  Eigenvalue= 0.2313396  =2.004
      6B ->  9B       -0.70568
      7B ->  8B        0.70568

 LASSW sees vector with all elements less than  0.100000000000E-03
 LASSW sees vector with all elements less than  0.100000000000E-03
 Excited state symmetry could not be determined.
 Eigenvector   3:  3.002-?Sym  Eigenvalue= 0.2704717  =2.003
      6B ->  8B        0.70711
      7B ->  9B       -0.70711
 SavETr:  write IOETrn=   770 NScale= 10 NData=  16 NLR=1 NState=    3 LETran=      64.
 The wavefunction is stable under the perturbations considered.

A time-dependent calculation will look for excitations based on the reference wave function. Therefore the solution will not be different. If you don't need to perform this calculation because you are interested in for example the UV/Vis spectrum, then there is no reason to do it.

I guess we have run the necessary scenarios now, so here are the results:

Command file              Functional            Energy / Hartree ( cycles )
bp86svp.rhf.sing          E(RB-P86)       =       -150.155033509 (      3 )
bp86svp.rhf.sing.maxc     No energy statement found.
bp86svp.rhf.sing.scf      E(RB-P86)       =       -150.154053425 (     16 )
bp86svp.rhf.sing.stab     E(RB-P86)       =       -150.154053425 (     16 )
bp86svp.rhf.sing.stab.opt E(UB-P86)       =       -150.201668045 (      4 )
bp86svp.rohf.trip         E(ROB-P86)      =       -150.213780181 (      9 )
bp86svp.uhf.quin.stable   E(UB-P86)       =       -149.651924155 (      8 )
bp86svp.uhf.sing          E(UB-P86)       =       -150.201668044 (      8 )
bp86svp.uhf.sing.stab     E(UB-P86)       =       -150.201668044 (      1 )
bp86svp.uhf.trip          E(UB-P86)       =       -150.215687576 (      8 )
bp86svp.uhf.trip.stab     E(UB-P86)       =       -150.215687574 (      1 )
bp86svp.uhf.trip.td       E(UB-P86)       =       -150.215687575 (      9 )

Lastly and maybe most important: with DFA it is never enough to just use one method. For more on this: DFT Functional Selection Criteria. Winter is coming!