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