Most organic chemistry deals with closed shell calculations, but every once in a while you want to calculate carbenes, free radicals or radical transition states coming from a homolytic bond break, which means your structure is now open shell.

Closed shell systems are characterized by having doubly occupied molecular orbitals, that is to say the calculation is ‘restricted’: Two electrons with opposite spin occupy the same orbital. In open shell systems, unrestricted calculations have a complete set of orbitals for the electrons with alpha spin and another set for those with beta spin. Spin contamination arises from the fact that wavefunctions obtained from unrestricted calculations are no longer eigenfunctions of the total spin operator <*S*^2>. In other words, one obtains an artificial mixture of spin states; up until now we’re dealing only with single reference methods. With each step of the SCF procedure the value of <*S*^2> is calculated and compared to *s*(*s*+1) where *s* is half the number of unpaired electrons (0.75 for a radical and 2.0 for triplets, and so on); if a large deviation between these two numbers is found, the then calculation stops.

Gaussian includes an annihilation step during SCF to reduce the amount of spin contamination but it’s not 100% reliable. Spin contaminated wavefunctions aren’t reliable and lead to errors in geometries, energies and population analyses.

One solution to overcome spin contamination is using Restricted Open Shell calculations (ROHF, ROMP2, etc.) for which singly occupied orbitals is used for the unpaired electrons and doubly occupied ones for the rest. These calculations are far more expensive than the unrestricted ones and energies for the unpaired electrons (the interesting ones) are unreliable, specially spin polarization is lost since dynamical correlation is hardly accounted for. The IOP(5/14=2) in Gaussian uses the annihilated wavefunction for the population analysis if acceptable but since Mulliken’s method is not reliable either I don’t advice it anyway.

The case of DFT is different since rho.alpha and rho.beta can be separated (similarly to the case of unrestricted ab initio calculations), but the fact that both densities are built of Kohn-Sham orbitals and not true canonical orbitals, compensates the contamination somehow. That is not to say that it never shows up in DFT calculations but it is usually less severe, of course for the case of hybrid functional the more HF exchange is included the more important spin contamination may become.

So, in short, for spin contaminated wavefunctions you want to change from restricted to unrestricted and if that doesn’t work then move to Restricted Open Shell; if using DFT you can use the same scheme and also try changing from hybrid to pure orbitals at the cost of CPU time. There is a last option which is using spin projection methods but I’ll discuss that in a following post.

Thank you for posting these helpful information about spin contamination. I enjoyed reading it and it helped through understanding the basics of it.

I’m glad the blog was useful to you.

Hello Prof. Barroso, very clear and nice post! I have a question concerning the spin-contamination: I am currently using the DFT method to calculate the S1 and T1 states of a purely organic molecule by employing the ΔSCF-DFT method, instead of making use of the TDDFT. Apart from the computational methodology employed, I am wondering how do you decide that your wavefunctions are too spin contaminated. Let’s say we are dealing with an open-shell system, thus the expected S^2 value is S^2=0.75. How far do I need to be to say I am contaminated (e.g. 5%, 10%) or 0.8, 1? is there a general rule-of-thumb in the DFT community? I am working @B3LYP level, so I would like to assess “how much” contamination I can possible get…

Thanks a lot in advance for any advice!

Hello Tommy,

That’s a great question, as far as I know the rule of thumb is around 10%, still most codes like Gaussian will abort the calculation if the contamination is too high but I ignore what the specific criterion of Gaussian is.

I hope this helps

I have done DFT (B3LYP) optimization for cobalt complexes, same level of theory and geom (optimized) used for single point analysis such as NBO, but the results of opt and sp calculation shows different energy and Mullikan spin density values. how we can rectify this problem?

Thanks in advance

Ho can we make input file on Gaussian to check spin contamination