1      subroutine dft_s2_value(geom,basis,symmetry,na,nb,
2     &     nbf,ga_densa,ga_densb)
3C     $Id$
4      implicit none
5#include "errquit.fh"
6#include "mafdecls.fh"
7#include "global.fh"
8#include "tcgmsg.fh"
9#include "stdio.fh"
10      integer geom,basis
11      logical symmetry
12      integer nbf,ga_densa,ga_densb,na,nb
13c
14      integer  ga_create_atom_blocked
15      external ga_create_atom_blocked
16c
17      integer iga_s,iga_tmp,iga_tmp2
18      double precision ana,anb
19      double precision s2,s2_exact
20      double precision one,zero,two
21      parameter(one=1.d0,zero=0.d0,two=2.d0)
22c
23c     prod S*Pb*S*Pa
24c
25c      call gacrea(geom,basis,nbf,nbf,'ovl ',iga_s,'atom')
26      iga_s = ga_create_atom_blocked(geom, basis, 'S2 temp 1')
27
28c      call gacrea(geom,basis,nbf,nbf,'ovl ',iga_tmp,'atom')
29      iga_tmp = ga_create_atom_blocked(geom, basis, 'S2 temp 2')
30
31c      call gacrea(geom,basis,nbf,nbf,'ovl ',iga_tmp2,'atom')
32      iga_tmp2 = ga_create_atom_blocked(geom, basis, 'S2 temp 3')
33
34      call ga_zero(iga_s)
35      call int_1e_ga(basis,basis,iga_s,'overlap',symmetry)
36
37      call ga_dgemm('N','N',nbf,nbf,nbf,one,iga_s,ga_densb,zero,iga_tmp)
38
39      call ga_dgemm('N','N',nbf,nbf,nbf,one,iga_tmp,iga_s,zero,iga_tmp2)
40
41      ana=dble(na)
42      anb=dble(nb)
43      s2_exact= 0.5d0*(ana-anb)*(.5d0*(ana-anb) + one)
44      s2= 0.25d0*(ana-anb)**two+0.5d0*(ana+anb) -
45     &     ga_ddot(iga_tmp2,ga_densa)
46
47c      call gadest(iga_s)
48      if (.not. ga_destroy(iga_s)) call errquit
49     &   ('dft_s2_value: could not destroy iga_s', 0, GA_ERR)
50c      call gadest(iga_tmp)
51      if (.not. ga_destroy(iga_tmp)) call errquit
52     &   ('dft_s2_value: could not destroy iga_tmp', 0, GA_ERR)
53c      call gadest(iga_tmp2)
54      if (.not. ga_destroy(iga_tmp2)) call errquit
55     &   ('dft_s2_value: could not destroy iga_tmp2', 0, GA_ERR)
56c
57      if(ga_nodeid().eq.0) then
58        call dft_header('Expectation value of S2:  ')
59        write (LuOut,12) s2,s2_exact
60   12   format('      <S2> = ',f11.4,' (Exact =',f11.4,')')
61        write(LuOut,*)
62      endif
63
64      return
65      end
66