spin contrained UHF calculations.

[topherson]

Dear Turbomole Users,

In the webpage http://www.uni-koeln.de/themen/chemie/software/turbomole/DOK/html/node107.html (Version 5.6)  there is a suggestion that one can constrain the spin of a UHF wavefunction calculation by writing in the control file:
"$uhf {spin constraint={real}}"

if I understood correctly, which I’m not sure I have, {real} should be the value one wishes the wavefunction to give for the S² operator.
If that is right, I am not sure I understood the format It requires for one to define the {real} variable.

I wished to optimise a doublet (s=1/2) state, so that S² =  s(s+1) = 0.75.
I tried "$uhf 0.75",  "$uhf spin constraint=0.75", and a few other combinations.
On all occasions, I got the warning:
"*** WARNING: Could not read properly <logical> value for $uhf ***"

On the web, I’ve seen output files that print the following statement:
"UHF modus switched on ! SPIN CONSTRAINT is : tau=    .00000"

Which suggests this can be done. On the other hand, the pdf version of the manual does not mention this option, which made me wonder whether the option exists.

I have run calculations where I start with a wavefunction that gives S² = 0.75.
but by the end of the RIDFT-MARIJ calculation I often get values above 2, presumably due to spin contamination.

Is there a way to constrain the spin during the optimization?

Thank you for your time!

[uwe]

Hello,

yes, the keyword does still work (I checked it for energy calculations using HF or DFT). The approach is based on:

http://www.sciencedirect.com/science/article/pii/000926149190405X

Spin contamination in single-determinant wavefunctions,  Jamie S. Andrews, Dylan Jayatilaka, Richard G.A. Bone, Nicholas C. Handy, Roger D. Amos, Chemical Physics Letters, Volume 183, Issue 5, September 1991, Pages 423–431

The real value you have to enter as option for the $uhf keyword is the Lagrange multiplier lambda given in the paper. The higher the value, the more ROHF-like (for HF) the calculation will be, zero is equivalent to UHF. Small values of about 0.01 seem to half the <S^2> error by almost keeping the UHF energy. Larger values (the authors use ranges between 0.01 and 500) bring down <S^2> to the exact expected value by moving the total energy towards the ROHF energy.

The original paper deals with HF and MP2 only, not DFT. To me it looks like it could be done in a similar way for DFT since both use a single-determinant wavefunction and the Lagrange method should work equally well. ridft and dscf both have this method implemented.

Just try it. Besided the <S^2> value which is printed to the output and to the control file ($ssquare), a good check are the occupation numbers of natural orbitals. If they are closer to 0, 1 and 2, the spin contamination did decrease.

Good luck, and please tell us about your results.

Regards,

Uwe

[topherson]

Dear Uwe,

Thank you for the detailed reply, it was very helpful. The method indeed works for RI-DFT.
For future users, I will just add a few observations with regards the spin constrained UHF method.

I tested it on metallic clusters with organic molecules.
Here is a sample of results:
S²                       tau
0.75070898     0.1
0.75650968     0.01
0.77702148     0.001
0.78549727     0.0001
0.79130847     0.00001

This is not the best example, since the contamination is not dramatic. However, I observed similar trends  for systems with more "contaminated" spins.
With tau > 0.5 it always gave ~ 0.75.
These deceptively small changes can have a dramatic effect of the PES, however. As they suggest in the paper, as tau -> infinity the wavefunction becomes closer to the ROHF wavefunction, and the energy can change significantly. For example, optimising the NO radical with different tau values can lead to staggeringly different values in the NO stretch bond.
It can also affect binding energies to the point where there is no binding.
Fortunately, I found that values around 0.01-0.02 typically keep S²  under ~5% of its correct value while keeping the surfaces close to the pure UHF.

Finally, the keyword in the control file must read:
$uhf on spin contraint={tau}

Hope this helps!