[darkdevil]

Hi!
I havn't been posting questions here for quite a while, hope everyone is safe in this pandemic :O

I was trained as an organic chemists and am still new to DFT using Gaussian.
I just did an open shell unrestricted DFT calculation of an organic radical cation (positively charged, doublet state) , and generated some really nice electron spin density pictures. What I am not sure about is the difference between the positive and negative ISO values of the lobes after generating the cube files. They turned out to be displayed as blue (positive) and green (negative) lobes wrapped around the molecule in GaussView. I am wondering what does the physical meaning of it?

Do the "balloons" mean where the extra "radical" electron will be most likely to locate? but why there are two different colors? I am guessing if blue represents the "up" spin state and green represents the "down" spin state?

Image is attached below:



Thanks and stay safe!

[pm133]

It looks to me like you might have some spin contamination here. Your output file will tell you that.
The colours are for spin up and spin down.
If there's no spin contamination then this might be showing that your unpaired electrons are delocalised.

[darkdevil]

It looks to me like you might have some spin contamination here. Your output file will tell you that.
The colours are for spin up and spin down.
If there's no spin contamination then this might be showing that your unpaired electrons are delocalised.

Thanks for your reply.
What exactly is a spin contamination and how do you know there is such contamination by looking at the ESD images? and how do I check it in the output file? I googled about it and have a rough idea that it is related to a different in spin multiplicity than the input spin state, but what causes this to happen?

[pm133]

Your question worries me slightly because it suggests you might be trying to use Gaussian as a black box tool. Is this the case? I'll assume for now that you know the maths behind this stuff and try to explain. If you don't know, you might need someone to help you.

Your spin state is not a pure doublet because in unrestricted shell calculations your determinants are not eigenfunctions of the S^2 operator.
Instead, it's an approximate doublet formed from a linear combination of doublet, quartet, sextet etc. spin states.
These higher spin states remain in the calculation and hence the name "spin contamination". As I remember, it forces a higher S^2 value than you would expect for a pure state and the energy increases as well.
Visually, you can get a clue to this potentially happening when you see your doublet smeared across a molecule when you don't expect it and specifically when you see phases of up and down spin split up and spread across the molecule as you see in your diagram.

Firstly, check your Gaussian output file to see the result of the calculation and make sure it is an unrestricted solution (which I'm sure it will be). It will tell you next to the final energy value.
Secondly there should be a value of S^2 alongside the result I think. That SHOULD be what you would expect for a pure doublet (0.7500) but it will be different and if I remember always higher than this because of the inclusion of higher spin states which you have not removed. I think Gaussian might try to use projection operators to remove as much of it as it can but I seem to remember this not being trivial at all and isn't guaranteed to work.

If you do have spin contamination you need to determine how bad it is by seeing how far from the pure doublet value you get. If it goes higher than 10% above this figure then you can't use the result because it's contaminated beyond repair. Either way, you need to check for wavefunction instabilities to check if a lower energy solution can be found. Use the Gaussian docs to see how to do that as I've been out the game for 3 years and I've forgotten. I think there's a keyword "...something...=stable".

My PhD was in open shell calculations using transition metals and I can assure you it isn't necessarily a trivial exercise to get to results you know you can trust. It's nice to see you being ambitious about using the tools but Gaussian isn't a block box for open shell stuff. You'll maybe need to start diving into the theoretical stuff as described in Szabo and Ostlund (Modern Quantum Chemistry). They do a great description of spin contamination on pages 104 to 107.

I'm afraid you've not asked a trivial question here so good luck with it. Let me know how you get on.

[darkdevil]

Your question worries me slightly because it suggests you might be trying to use Gaussian as a black box tool. Is this the case? I'll assume for now that you know the maths behind this stuff and try to explain. If you don't know, you might need someone to help you.

Your spin state is not a pure doublet because in unrestricted shell calculations your determinants are not eigenfunctions of the S^2 operator.
Instead, it's an approximate doublet formed from a linear combination of doublet, quartet, sextet etc. spin states.
These higher spin states remain in the calculation and hence the name "spin contamination". As I remember, it forces a higher S^2 value than you would expect for a pure state and the energy increases as well.
Visually, you can get a clue to this potentially happening when you see your doublet smeared across a molecule when you don't expect it and specifically when you see phases of up and down spin split up and spread across the molecule as you see in your diagram.

Firstly, check your Gaussian output file to see the result of the calculation and make sure it is an unrestricted solution (which I'm sure it will be). It will tell you next to the final energy value.
Secondly there should be a value of S^2 alongside the result I think. That SHOULD be what you would expect for a pure doublet (0.7500) but it will be different and if I remember always higher than this because of the inclusion of higher spin states which you have not removed. I think Gaussian might try to use projection operators to remove as much of it as it can but I seem to remember this not being trivial at all and isn't guaranteed to work.

If you do have spin contamination you need to determine how bad it is by seeing how far from the pure doublet value you get. If it goes higher than 10% above this figure then you can't use the result because it's contaminated beyond repair. Either way, you need to check for wavefunction instabilities to check if a lower energy solution can be found. Use the Gaussian docs to see how to do that as I've been out the game for 3 years and I've forgotten. I think there's a keyword "...something...=stable".

My PhD was in open shell calculations using transition metals and I can assure you it isn't necessarily a trivial exercise to get to results you know you can trust. It's nice to see you being ambitious about using the tools but Gaussian isn't a block box for open shell stuff. You'll maybe need to start diving into the theoretical stuff as described in Szabo and Ostlund (Modern Quantum Chemistry). They do a great description of spin contamination on pages 104 to 107.

I'm afraid you've not asked a trivial question here so good luck with it. Let me know how you get on.

Hi, I am doing organic electronics research. I am trying to compare the substituent effect using some model compounds in DFT. We observed a difference in stability for organic radical cation in terms of electrical conductivity. For example, compound A has basic thiophene unit and its conductivity would decay over a time scale of several days while we assum the radical cationic of compound A may somehow be destabilized by its electronic structure especially at the C-H bond on thiophene . For compound B, it gives out a stable conductivity for months, we assume that the methoxy groups can somhow stabilize the structure.

From the electron spin density diagram, we can see the spin is located on the C-H bond of the thiophene in compound A. I wonder if I can use this as an evidence to suggest the following radical cation dissociation at the C-H bond to give:

1) a hydrogen radical and a parent cation ( H^• + P⁺) ; or
2) a hydrogen cation (proton) and a parent radical ( H⁺ + P^•)

If the above mechanism really happen, a conductive compound can no longer have a delocalized polaron (or radical cation) along the pi-conjugated backbone, the electrical conductivity would then decrease.

The other image shows the spin density on compound B which is totally different from compound A that the methoxy group on the right hand side somehow act as a protecting group to prevent radical cation dissociation to happen.



but this is just my hypothesis right now.

[pm133]

Your work sounds interesting but for the reasons I described above, I am not sure you can be confident about that spin density calculation.
I would be especially careful about saying the spin is in the middle of the C-H bond. Have you performed an NBO analysis to confirm that?

[darkdevil]

Your work sounds interesting but for the reasons I described above, I am not sure you can be confident about that spin density calculation.
I would be especially careful about saying the spin is in the middle of the C-H bond. Have you performed an NBO analysis to confirm that?

I've also done an NBO analysis for the two compounds.
I was wondering when we oxidize a compound, do we first remove the electron from the Beta MO first?
This is the NBO of the C-H bond in the radical cationic compound A and compound B. Red lobes are UP spin and Green lobes are DOWN spin.

If we use the open shell model, the first electron oxidation (removing 1 electron from the HOMO or valence shell) would be removing the electron from the Beta MO which has the higher energy. Is this correct?

If so, the electron removed in compound A would be located at the single C-H bond (the green lobe), and the resulting product would be a weak C-H bond or a broken C-H bond.
For compound B, the NBO looks totally different, and the electron being removed would be from either the conjugated backbone or the lone pair electrons of oxygen instead of a single bond and prevent the radical dissociation reactions to occur.

Please correct me if I am wrong. Much appreciated!

[pm133]

Any NBO or other analysis is totally dependent on your calculation being free of spin contamination and wavefunction instabilities.

Is your wavefunction stable?
What is your S^2 value?

BTW the NBO analysis will tell you where the electrons are located and in what quantities. For that, you'll need to read the logfile. Orbital pictures are notoriously difficult to interpret physically and can be misleading.

[darkdevil]

Any NBO or other analysis is totally dependent on your calculation being free of spin contamination and wavefunction instabilities.

Is your wavefunction stable?
What is your S^2 value?

BTW the NBO analysis will tell you where the electrons are located and in what quantities. For that, you'll need to read the logfile. Orbital pictures are notoriously difficult to interpret physically and can be misleading.

Thanks for replying. This is from the output file after doing the NBO analysis for compound A

The electronic state of the initial guess is 2-A.
 Initial guess <Sx>= 0.0000 <Sy>= 0.0000 <Sz>= 0.5000 <S**2>= 0.7500 S= 0.5000
 Requested convergence on RMS density matrix=1.00D-08 within 128 cycles.
 Requested convergence on MAX density matrix=1.00D-06.
 Requested convergence on             energy=1.00D-06.
 No special actions if energy rises.
 Integral accuracy reduced to 1.0D-05 until final iterations.
 EnCoef did     8 forward-backward iterations
 EnCoef did     6 forward-backward iterations
 EnCoef did     3 forward-backward iterations
 Initial convergence to 1.0D-05 achieved.  Increase integral accuracy.
 SCF Done:  E(UB3LYP) =  -2361.43404646     A.U. after   28 cycles
             Convg  =    0.5853D-08             -V/T =  2.0064
 <Sx>= 0.0000 <Sy>= 0.0000 <Sz>= 0.5000 <S**2>= 0.7860 S= 0.5178
 <L.S>= 0.000000000000E+00
 Annihilation of the first spin contaminant:
 S**2 before annihilation     0.7860,   after     0.7508

What does the number 0.7860 and 0.7508 mean?

[pm133]

Those are your S^2 values before and after Gaussian attempted projection methods to remove the spin contamination.
You can see the contamination in these numbers but projection seems to have removed it to an excellent degree. This is good news as you can rule out spin contamination. Can you see that the S^2 value of 0.7860 is 4.8% spin contamination and that the corrected value of 0.7508 is now only 0.1%?

Please can you now go and ensure that your previous non-NBO calculation gives the same result and let me know.
Also, can you scour the NBO analysis log file and tell me how much of the doublet is at each location in your molecule. That will indicate whether it is truly delocalised or sitting squarely on that C-H bond. You need the physical numbers for this. Sorry, I can't remember the exact location in the NBO file this appears in.

I am still very concerned about that first spin density calculation you posted and you need to get to the bottom of that.

[pm133]

One other thought that you should try.
Your isoval seems a little low at 0.0004 and all those coloured lobes in your first picture might be very small numbers. Isoval is a cutoff number for showing density.

Try increasing it a little until to see if all those coloured lobes are actually just "noise". If they all rapidly fade away and you end up with just one lobe, that might be your unpaired electron. You can confirm that by seeing how many electrons your logfile thinks is at that location.

Don't forget your wavefunction stability calculation. Always re-check the S^2 value after every re-calculation.

[darkdevil]

Those are your S^2 values before and after Gaussian attempted projection methods to remove the spin contamination.
You can see the contamination in these numbers but projection seems to have removed it to an excellent degree. This is good news as you can rule out spin contamination. Can you see that the S^2 value of 0.7860 is 4.8% spin contamination and that the corrected value of 0.7508 is now only 0.1%?

Please can you now go and ensure that your previous non-NBO calculation gives the same result and let me know.
Also, can you scour the NBO analysis log file and tell me how much of the doublet is at each location in your molecule. That will indicate whether it is truly delocalised or sitting squarely on that C-H bond. You need the physical numbers for this. Sorry, I can't remember the exact location in the NBO file this appears in.

I am still very concerned about that first spin density calculation you posted and you need to get to the bottom of that.

1) The non-NBO calculations and all the charged compounds return a value of 0.750X after the correcetion by Gaussian
2) For how much doublet is at each location in the molecule, I am still figuring it out. Could it be something related to Lewis occupancies? Still looking around in the file

[pm133]

I don't know what Lewis occupancies are.
I think you might be looking for "spin densities" or words to that effect.

It should list atoms and the net spin density on each atom and it should also list bonds and the net spin densities on those.
It's a very big file with a lot of information so this might take you a while.

Whatever happens here, you have re-ignited my interest in computational stuff so thanks for that :-D

[darkdevil]

I don't know what Lewis occupancies are.
I think you might be looking for "spin densities" or words to that effect.

It should list atoms and the net spin density on each atom and it should also list bonds and the net spin densities on those.
It's a very big file with a lot of information so this might take you a while.

Whatever happens here, you have re-ignited my interest in computational stuff so thanks for that :-D

That's very nice of you! It is surprising to know that the computational tools can be so powerful that it might be able to explain something that is hard to prove by wet-lab experiments. I am still digging into the details of it. Hopefully, I can update you about this.

[darkdevil]

I don't know what Lewis occupancies are.
I think you might be looking for "spin densities" or words to that effect.

It should list atoms and the net spin density on each atom and it should also list bonds and the net spin densities on those.
It's a very big file with a lot of information so this might take you a while.

Whatever happens here, you have re-ignited my interest in computational stuff so thanks for that :-D

I found something called the Mulliken Atomic Spin densities:
So I copied the following data to an excel spreadsheet and the sum of the densities are equal to 1.00000 just as it says at the bottom in the list. I am a bit confused now as it say atomic spin density. Is it different from electron spin? but they should be related right?

If I sum up all the positive atomic spin densities values I will get +1.480265
and the sum of all negative values would be -0.480261
So the net value is exactly +1 which is the charge of the model. With this data am I able to tell where is the + charge distributed among the molecule?
For example carbon number 1 (C1) has an atomic spin density of 0.148771, the positive charge has a 0.148771/1.480265*100 = 12.6% distributed on the C1 atom?

but how can this/ is it even related to the electron spin of the radical in the molecule?

Mulliken atomic spin densities:
              1
     1  C    0.148771
     2  C   -0.043949
     3  C    0.103003
     4  C   -0.017994
     5  S   -0.015241
     6  C   -0.107084
     7  C    0.130569
     8  C    0.076180
     9  C   -0.027788
    10  C    0.051472
    11  C   -0.033282
    12  C    0.057998
    13  C   -0.034911
    14  C    0.164456
    15  C   -0.009643
    16  N   -0.010091
    17  O    0.059706
    18  C   -0.028911
    19  C    0.070158
    20  C    0.173872
    21  S   -0.031511
    22  C    0.021204
    23  C    0.140449
    24  C   -0.068433
    25  C    0.180296
    26  S   -0.025401
    27  O    0.046927
    28  O    0.034250
    29  C   -0.003824
    30  C   -0.002676
    31  H    0.001271
    32  H   -0.004552
    33  H    0.004127
    34  H   -0.003524
    35  H    0.000957
    36  H   -0.002384
    37  H    0.001252
    38  H    0.000064
    39  H    0.000578
    40  H    0.002289
    41  H   -0.008971
    42  H    0.003101
    43  H   -0.000059
    44  H    0.003101
    45  H    0.002107
    46  H    0.002107
    47  H   -0.000032
 Sum of Mulliken atomic spin densities =   1.00000