# Calculating singlet-triplet gap and open-shell singlet character

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?

• 1.a you need to use same level of theory (open shell or unrestricted) for both S0 and T1. 1.b you should consider using non-dft, as spin-polarized dft has somewhat questionably fundation. It is OK for big systems when there is no choice, but for systems of 18 atoms something like ump2/6-311++g** should be used, at least for final energies. – permeakra Jul 29 '16 at 18:43
• @permeakra 1.a not necessarily. If either formalism is a good description it does not matter. If not then both are inappropriate. See also chemistry.stackexchange.com/q/49174/4945 1.b UMP2 will be worse than any DFT (and that is not considering that there could be a transitionmetal in there). Since OP did not specify what kind of atoms these 18 are, it is hard to judge if a higher method is suitable. – Martin - マーチン Oct 14 '16 at 4:58
• @Martin-マーチン Hey, They are pi-conjugated systems containing phosphorus, nitrogen, sulfur, carbon, and hydrogen. I have been reading about some methods such as broken symmetry DFT, and spin-flip TD-DFT, which help in cases of potential open-shell singlet ground states. – kimdazar Jan 13 '17 at 21:46

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
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

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
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)

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)
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!

• Martin!! . Thank you for such a detailed and illustrative explanation!!! I am still confused as of why you specify scf(xqc) together with guess=read in the stability calcs. 1) I naturally understand you exemplified a typical HF procedure in which the SCF calculates/creates/guesses a wavefunction until the energy and density (orbitals!) criteria are met, plus the fact that it is a very simple diatomic molecule (thus no optimization); if I was doing a DFT on the pi-system I would naturally need to first optimize then,can I use that optimized geometry for all the wavefunction stab calculations? – kimdazar Jan 14 '17 at 0:54
• To be clear, the opt keyword on the stable=opt portion of a stability input job does not refer to a geometrical optimization, right? Thanks a lot. – kimdazar Jan 14 '17 at 1:01
• @kimdazar I think there was no particular reason to specify xqc after reading in the orbitals other than laziness, or maybe the stability analysis needed it. I don't remember; it shouldn't do any harm though. You should use the optimised geometry for stability calculations. If you find instabilities you need to optimise again with the new wave function. Stable=opt is only for single point calculations; yes. – Martin - マーチン Jan 14 '17 at 1:25