1c===============================================================c 2c c 3c Socrates - an NWChem module for teaching c 4c Copyright © 2009 Jeff Hammond c 5c c 6c Developed by: c 7c c 8c Jeff R. Hammond c 9c Leadership Computing Facility c 10c Argonne National Laboratory c 11c jhammond@mcs.anl.gov c 12c c 13c===============================================================c 14 logical function socrates_scf_ga(rtdb,geom,num_alfa,num_beta) 15c 16c $Id$ 17c 18 implicit none 19#include "mafdecls.fh" 20#include "global.fh" 21#include "rtdb.fh" 22#include "errquit.fh" 23#include "stdio.fh" 24#include "bas.fh" 25#include "geom.fh" 26#include "sym.fh" 27#include "schwarz.fh" 28c 29c object handles 30c 31 integer rtdb ! RTDB handle 32 integer geom ! GEOM handle 33 integer basis ! BASIS handle 34c 35c GA variables 36c 37 integer bytes ! Number of bytes in a double 38c 39c INT variables 40c 41 integer i,j,k,l 42 integer nbas ! number of basis sets used 43 parameter (nbas=1) 44 integer nbf,nshells 45 integer s_int2e,k_int2e,l_int2e 46 integer s_scr2e,k_scr2e,l_scr2e 47 integer ish,ilo,ihi,irng 48 integer jsh,jlo,jhi,jrng 49 integer ksh,klo,khi,krng 50 integer lsh,llo,lhi,lrng 51 integer ao_offset 52c 53c SCF variables 54c 55 integer g_smat 56 integer g_hmat 57 integer g_dmata 58 integer g_dmatb 59 integer g_fmata 60 integer g_fmatb 61 integer num_elec,num_alfa, num_beta 62c 63c timing variables 64c 65 double precision cpu ! CPU sec counter 66 double precision wall ! WALL sec counter 67c 68c INT variables 69c 70 double precision int_thresh 71c 72c SCF variables 73c 74 double precision tol2e 75 double precision schwarz_ij,schwarz_kl 76c 77c primitive variables 78c 79 logical nodezero ! am I node 0? 80 logical debug ! debug mode? 81 logical stat ! status 82 logical use_symm ! use symmetry in fock building 83c 84c SYM variables 85c 86 character*8 group 87c 88c external routines 89c 90 logical int_normalize 91 external int_normalize 92 integer ga_create_atom_blocked 93 external ga_create_atom_blocked 94 double precision dnrm2 95 external dnrm2 96c 97 parameter (use_symm = .false.) 98c 99#ifdef DEBUG_PRINT 100 debug = (GA_NodeId().eq.0) ! debug print on nodezero only 101c debug = .true. ! debug print everywhere 102#else 103 parameter (debug = .false.) 104#endif 105c 106 if (debug) then 107 write(LuOut,*) 'top of socrates_scf_ga' 108 call util_flush(LuOut) 109 endif 110c 111 socrates_scf_ga = .false. 112c 113 nodezero=(ga_nodeid().eq.0) 114c 115c bytes = ma_sizeof(mt_dbl,1,mt_byte) 116c 117c===============================================================c 118c c 119c Initialize BASIS and INTEGRAL objects c 120c c 121c===============================================================c 122c 123c --------- 124c Basis set 125c --------- 126c 127 if (.not.bas_create(basis,'ao basis')) then 128 call errquit(__FILE__,__LINE__,BASIS_ERR) 129 endif 130 if (.not.bas_rtdb_load(rtdb,geom,basis,'ao basis')) then 131 call errquit(__FILE__,__LINE__,BASIS_ERR) 132 endif 133 if (.not.bas_numbf(basis,nbf)) then 134 call errquit(__FILE__,__LINE__,BASIS_ERR) 135 endif 136 if (.not.bas_numcont(basis,nshells)) then 137 call errquit(__FILE__,__LINE__,BASIS_ERR) 138 endif 139c 140 if (nodezero) then 141 call bas_print_labels(basis) 142 write(LuOut,*) 143 call util_flush(LuOut) 144 endif 145c 146c --------- 147c Integrals 148c --------- 149c 150 if (.not.rtdb_get(rtdb,'socrates:tol2e',mt_dbl,1,tol2e)) then 151 tol2e = 1.0d-10 152 endif 153 if (.not.rtdb_get(rtdb,'socrates:int_thresh',mt_dbl,1, 154 1 int_thresh)) then 155 int_thresh = 1.0d-20 156 endif 157 call int_acc_set(int_thresh) 158c 159 call int_init(rtdb,nbas,basis) 160 call schwarz_init(geom,basis) 161 if (.not.int_normalize(rtdb,basis)) then 162 if (nodezero) write(LuOut,*) 'int_normalize failed' 163 call errquit(__FILE__,__LINE__,INT_ERR) 164 endif 165c 166 call int_mem_2e4c(s_int2e,s_scr2e) 167 if (.not.ma_push_get(mt_dbl,s_int2e,'int2e',l_int2e,k_int2e)) then 168 call errquit(__FILE__,__LINE__,MA_ERR) 169 endif 170 if (.not.ma_push_get(mt_dbl,s_scr2e,'scr2e',l_scr2e,k_scr2e)) then 171 call errquit(__FILE__,__LINE__,MA_ERR) 172 endif 173c 174c===============================================================c 175c c 176c Begin SCF code c 177c c 178c===============================================================c 179c 180c initialize overlap matrix 181c 182 g_smat = ga_create_atom_blocked(geom, basis,'smat') 183 call int_1e_ga(basis,basis,g_smat,'overlap',use_symm) 184c 185c initialize Hcore matrix 186c 187 g_hmat = ga_create_atom_blocked(geom, basis,'hmata') 188 call int_1e_ga(basis,basis,g_hmat,'kinetic',use_symm) 189 call int_1e_ga(basis,basis,g_hmat,'potential',use_symm) 190c 191c initialize density matrices 192c 193 g_dmata = ga_create_atom_blocked(geom, basis,'dmata') 194 g_dmatb = ga_create_atom_blocked(geom, basis,'dmatb') 195! 196! initialize density matrix somehow here 197! 198c 199c create Fock matrices 200c 201 g_fmata = ga_create_atom_blocked(geom, basis,'fmata') 202 g_fmatb = ga_create_atom_blocked(geom, basis,'fmatb') 203 204c 205c 206c 207 208 209c 210#ifdef DEBUG_PRINT 211 if (nodezero) then 212 write(LuOut,*) 'Printing all non-negligible AOs' 213 call util_flush(LuOut) 214 endif 215#endif 216 do ish=1,nshells 217 if (.not.bas_cn2bfr(basis,ish,ilo,ihi)) then 218 call errquit(__FILE__,__LINE__,BASIS_ERR) 219 endif 220 irng = ihi - ilo + 1 221 do jsh=1,nshells 222 if (.not.bas_cn2bfr(basis,jsh,jlo,jhi)) then 223 call errquit(__FILE__,__LINE__,BASIS_ERR) 224 endif 225 jrng = jhi - jlo + 1 226 schwarz_ij = schwarz_shell(ish,jsh) 227 do ksh=1,nshells 228 if (.not.bas_cn2bfr(basis,ksh,klo,khi)) then 229 call errquit(__FILE__,__LINE__,BASIS_ERR) 230 endif 231 krng = khi - klo + 1 232 do lsh=1,nshells 233 if (.not.bas_cn2bfr(basis,lsh,llo,lhi)) then 234 call errquit(__FILE__,__LINE__,BASIS_ERR) 235 endif 236 lrng = lhi - llo + 1 237 schwarz_kl = schwarz_shell(ksh,lsh) 238 if ((schwarz_ij*schwarz_kl).ge.tol2e) then 239 call int_2e4c(basis,ish,jsh,basis,ksh,lsh, 240 1 s_scr2e,dbl_mb(k_scr2e), 241 2 s_int2e,dbl_mb(k_int2e)) 242#ifdef DEBUG_PRINT 243 call socrates_print_ao2e(ilo,ihi,jlo,jhi, 244 1 klo,khi,llo,lhi, 245 2 dbl_mb(k_int2e),tol2e) 246#endif 247 endif 248 enddo 249 enddo 250 enddo 251 enddo 252#ifdef DEBUG_PRINT 253 if (nodezero) then 254 write(LuOut,*) 'End of integral printing' 255 write(LuOut,*) 256 call util_flush(LuOut) 257 endif 258#endif 259c 260c ========== 261c Deallocate 262c ========== 263c 264 stat = ga_destroy(g_dmata) 265 if (.not.stat) then 266 call errquit(__FILE__,__LINE__,GEOM_ERR) 267 endif 268 stat = ga_destroy(g_dmatb) 269 if (.not.stat) then 270 call errquit(__FILE__,__LINE__,GEOM_ERR) 271 endif 272 stat = ga_destroy(g_fmata) 273 if (.not.stat) then 274 call errquit(__FILE__,__LINE__,GEOM_ERR) 275 endif 276 stat = ga_destroy(g_fmatb) 277 if (.not.stat) then 278 call errquit(__FILE__,__LINE__,GEOM_ERR) 279 endif 280 stat = ga_destroy(g_hmat) 281 if (.not.stat) then 282 call errquit(__FILE__,__LINE__,GEOM_ERR) 283 endif 284 stat = ga_destroy(g_smat) 285 if (.not.stat) then 286 call errquit(__FILE__,__LINE__,GEOM_ERR) 287 endif 288c 289c===============================================================c 290c c 291c Close any open object references and stack blocks c 292c c 293c===============================================================c 294c 295#ifdef DETAILED_FREE 296 if (.not.ma_pop_stack(l_scr2e)) then 297 call errquit('ncc_driver: MA problem scr2e ',0,MA_ERR) 298 endif 299 if (.not.ma_pop_stack(l_int2e)) then 300 call errquit('ncc_driver: MA problem int2e ',0,MA_ERR) 301 endif 302#else 303 if (.not. ma_chop_stack(l_int2e)) then 304 call errquit('ncc_driver: stack corrupt ',0, MA_ERR) 305 endif 306#endif 307c 308 call schwarz_tidy() 309 call int_terminate() 310c 311 if (.not.bas_destroy(basis)) then 312 call errquit(__FILE__,__LINE__,BASIS_ERR) 313 endif 314c 315c===============================================================c 316c c 317c THE END c 318c c 319c===============================================================c 320c 321 socrates_scf_ga = .true. 322c 323 if (debug) then 324 write(LuOut,*) 'end of socrates_scf_ga' 325 call util_flush(LuOut) 326 endif 327c 328 return 329c 330 2001 format(/,3x,'for ',a5,' MO vectors:',/, 331 1 3x,'number of basis functions = ',i8,/, 332 2 3x,'number of molecular orbitals = ',i8,/) 333 334c 335 end 336