Subject: Spin contamination in DFT - FORTRAN code

From: Jan Labanowski <[email-masked]>
Date: Sun, 1 Oct 1995 12:16:11 -0400

Once more about spin contamination in DFT. I hope that my FORTRAN is
better than my broken English. I wrote a routine to compute it.
Can you look at it and see if it is a nonsense or a pearl {:-)}?
Maybe it is easier for you to look at it than at my previous message
(I actually messed up in the previous message, since I did not square
the SAB[i,j] at the very end). Sorry for my oldfashioned FORTRAN,
but I learned it before optimizing compilers on a 64k machine and old
habits die slow (but I still have a few years till retirement, I simply
started late... {:-)}).

Jan Labanowski
[email-masked]

---------------- Variable names as in my previous message -------
      SUBROUTINE SPNCNT(S2EXCT, S2EXPT)

      IMPLICIT REAL*8(A-H,O-Z)
      INCLUDE "param.h"

      COMMON/BASIS/HMAT(NFOCK),TVECTR(NAUXFN),NDTT,NCNTRT,NNNS,NNNP,NNND

      COMMON/ORBEN/ELGALP(NCONTS),ELGBET(NCONTS),DOCA(NCONTS),
     &             DOCB(NCONTS),NLA,NLB

      COMMON/MATRX/SINV(NCONTS,NCONTS),FOCK1(NFOCK),FOCK2(NFOCK),
     &             FOCK3(NFOCK)
      COMMON/MODAT/AMAT(NFOCK),BMAT(NFOCK),
     &             CALPHA(NCONTS,NCONTS),CBETA(NCONTS,NCONTS)
      COMMON/INFOA/CENT(3,NATOMS),CHARGE(NATOMS),NATOMC(NATOMS),NCENTR,
     &             NALPHA,NBETA,NSTATE,NHOLEA(NCONTS),NHOLEB(NCONTS)

      DIMENSION WORK(NCONTS, NCONTS), WVEC(NCONTS)

C  --  Read in overlap matrix
      REWIND 43
      DO 100 I=1,NCNTRT
  100   CALL MYREAD(43,SINV(1,I),NCNTRT)

C  --  "exact spin square"
      S2EXCT = 0.25D0*(NALPHA - NBETA)*(NALPHA - NBETA + 2.0D0)

C  --  Calculate MO overlap matrix SAB = CALPHA' * SINV * CBETA
C  --  CALPHA' is a transpose of CALPHA

C  --  Make a matrix WORK(1:NCNTRT,1:NLB) = SINV * CBETA
      DO 400 K1 = 1, NCNTRT
        DO 150 K2 = 1, NCNTRT
C --        since SINV(K1,K2) = SINV(K2,K1) and FORTRAN orders matrices along
C --        columns, it is better for vectorization to do SINV(K2,K1)
          WVEC(K2) = SINV(K1,K2)
  150   CONTINUE
        DO 300 J =  1, NLB
          SUM = 0.0
          DO 200 K2 = 1, NCNTRT
            SUM = SUM + CBETA(K2,J)*WVEC(K2)
  200     CONTINUE
          WORK(K1,J) = SUM
  300   CONTINUE
  400 CONTINUE

C  --  Calculate SAB = CALPHA' * WORK (save it in SINV)
      DO 700 I = 1, NLA
        DO 600 J = 1, NLB
          SUM = 0.0
          DO 500 K1 = 1, NCNTRT
             SUM = SUM + CALPHA(K1,I)*WORK(K1,J)
  500     CONTINUE
          SINV(I,J) = SUM
  600   CONTINUE
  700 CONTINUE

C  --  Calculate expectation value of square of total spin
      S2EXPT = S2EXCT + NBETA
      DO 900 I = 1, NLA
        DO 800 J = 1, NLB
          X = SINV(I,J)
          S2EXPT = S2EXPT - DOCA(I)* DOCB(J) * X*X
  800   CONTINUE
  900 CONTINUE

C  --  Restore overlap matrix
      REWIND 43
      DO 1000 I=1,NCNTRT
 1000   CALL MYREAD(43,SINV(1,I),NCNTRT)

      RETURN
      END