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