Subject: Spin Contamination in DFT -- How to compute?

From: Jan Labanowski <[email-masked]>
Date: Sat, 30 Sep 1995 18:43:19 -0400

Dear Netters,

I am trying to compute spin contamination in DFT. My starting point
is the A. Szabo, N.S. Ostlund -- Modern Quantum Chemistry", pp. 104-107,
McGraw-Hill, New York, 1989.

Background:
Multiplicity is
   MULT = 2*S + 1
where S is the total spin quantum number having possible values
  0.0, 0.5, 1.0, 1.5, ....
S is simply 0.5 times the number of unpaired electrons. Assuming that
NALPHA is the number of electrons with spin alpha, and NBETA is the
number of electrons with spin beta, S = 0.5 * (NALPHA - NBETA) if
we assume that NALPHA >= NBETA.
The eigenvalue of the operator of square of total spin is denoted
by S2, and is NOT EQUAL to S*S but is equal to:
   S2 = S * (S + 1)
To calculate S from S2 the following formula (solution of quadratic
equation) can be used:
     S = 0.5*(SQRT(1.0 + 4.0*S2) - 1.0)
Another quantum number:
   MS = -S, -S+1, ... S-1, S
is the quantum number (and also the eigenvalue) of the z-component
of the total spin. Now,
   S2 = 0.5*(NALPHA-NBETA)*(0.5*(NALPHA-NBETA) + 1)
i.e., the "exact value" of eigenvalue of square of spin (SEXCT) should be:
   S2EXCT = S2 = 0.25*(NALPHA - NBETA)(NALPHA - NBETA + 2)
and the "exact" multiplictity should be:
   MULTEX = NALPHA - NBETA + 1
I assume that the formula for the expectation value of S2 (S2EXPT),
given on page 107 of Szabo & Ostlund as:
   S2EXPT = S2EXCT + NBETA - sum{i=1,N} sum{j=1,N} |(SAB[i,j])|^2
for the case of NALPHA >= NBETA (the meaning of N is not defined),
can be extended to fractional occupation numbers as:
S2EXPT = S2EXCT + NBETA -
   sum{i=1,NLA} sum{j=1,NLB} DOCA[i] * DOCB[j] * |(SAB[i,j])^2|
where NLA is the number of occupied (possibly with fractional
occupation numbers) ALPHA MO's, NLB is the number of occupied (possibly
fractionally) BETA MO's, DOCA[i] is the (fractional) occupation number on i-th
ALPHA orbital, DOCB[j] is the (possibly) fractional occupation number on j-th
orbital, SAB[i,j] = <PSIA[i]|PSIB[j]> is the overlap integral between
PSIA[i] (the i-th MO for spin ALPHA) and PSIB[j] (the j-th MO for spin BETA).
Taking this into account that
    PSIA[i] = sum{k=1,NCNTRT} CALPHA[k,i]*PHI[k]
    PSIB[j] = sum{k'=1,NCNTRT} CBETA[k',j]*PHI[k']
where NCNTRT is the number of basis functions (contractions), CALPHA[k,i]
is the MO coefficient in the i-th ALPHA MO for k-th basis function,
CBETA[k',j] is the MO coefficient in the j-th BETA MO for k'-th basis function,
and PHI[k] is the k-th basis function; one gets:

SAB[i,j] =  sum{k=1,NCNTRT} sum{k'=1,NCNTRT} CALPHA[k,i]*CBETA[k',j]*S[k,k']

and hence:

   S2EXPT = S2EXCT + NBETA - 
      sum{i=1,NLA} sum{j=1,NLB} sum{k=1,NCNTRT} sum{k'=1,NCNTRT}
      DOCA[i] * CALPHA[k,i] * DOCB[j]  * CBETA[k',j] * S[k,k']
where S[k,k'] is the oridinary overlap integral between k-th and k'-th basis
function: S[k,k'] = <PHI[k]|PHI[k']>.

Now, is this formula correct? Or I messed it up...

Yes, I understand that fractional occupation numbers are suspect,
and I really should not have them when my SCF converges.
However, I need to calculate SOMETHING for the case when I do level smearing.
On the other hand, I am not that sure if the formula is even correct for
the intergal occupation numbers.

Can you help?

Jan Labanowski
[email-masked]

P.S. Sorry, to some who got a portion of this message on my first try.