Calculating singlettriplet gap and openshell singlet character
There is a related discussion on which method to use for a singlet triplet guess, which I encourage you to also read: U or ROmethod for SingletTriplet 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 singlettriplet gap. If you are more interested in the electronic excitations a timedependent 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 multireference 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 def2SVP 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.100000000000E03 LASSW sees vector with all elements less than 0.100000000000E03 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 timedependent 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(RBP86) = 150.155033509 ( 3 ) bp86svp.rhf.sing.maxc No energy statement found. bp86svp.rhf.sing.scf E(RBP86) = 150.154053425 ( 16 ) bp86svp.rhf.sing.stab E(RBP86) = 150.154053425 ( 16 ) bp86svp.rhf.sing.stab.opt E(UBP86) = 150.201668045 ( 4 ) bp86svp.rohf.trip E(ROBP86) = 150.213780181 ( 9 ) bp86svp.uhf.quin.stable E(UBP86) = 149.651924155 ( 8 ) bp86svp.uhf.sing E(UBP86) = 150.201668044 ( 8 ) bp86svp.uhf.sing.stab E(UBP86) = 150.201668044 ( 1 ) bp86svp.uhf.trip E(UBP86) = 150.215687576 ( 8 ) bp86svp.uhf.trip.stab E(UBP86) = 150.215687574 ( 1 ) bp86svp.uhf.trip.td E(UBP86) = 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!
Related videos on Youtube
kimdazar
Updated on October 18, 2020Comments

kimdazar about 3 years
I am trying to calculate the singlettriplet energy gap ($S_0T_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/6311g(d,p) integral(ultrafinegrid) ... 0 1 ...
and for the energy of the triplet $T_1$ state
opt=tight rob3lyp/6311g(d,p) integral(ultrafinegrid) ... 0 3 ...
Does the energy difference of those two calculations have any physical significance? (i.e. should have I used
...tdub3lyp...0 3 ...
instead?)I believe a system I am currently working on may have some openshell 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/6311+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?

permeakra over 7 years1.a you need to use same level of theory (open shell or unrestricted) for both S0 and T1. 1.b you should consider using nondft, as spinpolarized 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/6311++g** should be used, at least for final energies.

Martin  マーチン about 7 years@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.

kimdazar almost 7 years@Martinマーチン Hey, They are piconjugated systems containing phosphorus, nitrogen, sulfur, carbon, and hydrogen. I have been reading about some methods such as broken symmetry DFT, and spinflip TDDFT, which help in cases of potential openshell singlet ground states.

kimdazar almost 7 yearsMartin!! . 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 pisystem I would naturally need to first optimize then,can I use that optimized geometry for all the wavefunction stab calculations?

kimdazar almost 7 yearsTo 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.

Martin  マーチン almost 7 years@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.