remove more of the BGJ spaghetti

H2O (RHF) runs correctly through the following:

task scf optimize
task scf freq
task dft optimize
task dft freq

this is hardly sufficient but identified a lot of issues.

the IBM XLF compiler will identify a bunch of argument mismatch issues once i use it.

Signed-off-by: Jeff Hammond <[email protected]>

src/cphf/cphf_poliz.F

@@ -45,7 +45,7 @@ logical function cphf_poliz(rtdb)
 c
 c     Get MO vectors etc.
 c
-      hf_job = .not. bgj_have_xc()
+      hf_job = .not. bgj_have_xc(rtdb)
       if(.not.rtdb_cget(rtdb, 'title', 1, title)) title = ' '
       if(.not.geom_create(geom, 'geometry'))
      $     call errquit('cphf_poliz: geom_create?', 0, GEOM_ERR)
@@ -70,7 +70,7 @@ logical function cphf_poliz(rtdb)
      $        call errquit('cphf_poliz: no DFT MO vectors',0, RTDB_ERR)
 c         write(*,*)'*** cphf_poliz: dft movecs: ',movecs
 c !!! Hack to set up scftype properly for later code
-         if (bgj_restricted()) then
+         if (bgj_restricted(rtdb)) then
             scftype = 'RHF'
          else
             scftype = 'UHF'
@@ -142,7 +142,7 @@ logical function cphf_poliz(rtdb)
             nob=noa
          endif
       endif
-      if (bgj_print() .gt. 0)
+      if (bgj_print(rtdb) .gt. 0)
      &     write(*,*)'*** cphf_poliz: noa nob',noa,nob
       nopen = 0
       nva=nmo-noa
@@ -152,7 +152,7 @@ logical function cphf_poliz(rtdb)
       else if (scftype .eq. 'RHF' .or. scftype.eq.'ROHF') then
          nvirt = nmo - nclosed - nopen
          vlen = nclosed*(nopen+nvirt) + nopen*nvirt
-         if (bgj_print() .gt. 0)
+         if (bgj_print(rtdb) .gt. 0)
      &   write(*,*)'nvirt,nmo,nclosed,nopen',nvirt,nmo,nclosed,nopen
       else
          call errquit('cphf: unknown SCF type',0,
@@ -237,7 +237,7 @@ logical function cphf_poliz(rtdb)
 c
       enddo
 c
-      if (bgj_print() .gt. 0)
+      if (bgj_print(rtdb) .gt. 0)
      &     write(*,*)'*** cphf_poliz: done with x y z loop'
       if (.not. ga_destroy(g_vecs_a)) call errquit('cphf_poliz: ga',1,
      &       GA_ERR)

src/cphf/cphf_solve.F

@@ -76,7 +76,7 @@ subroutine cphf_solve(rtdb)
      $        call errquit('cphf_solve: no DFT MO vectors',0, RTDB_ERR)
 c !!! Hack to set up scftype properly for later code !!! ? need this?
 #if 0
-         if (bgj_restricted()) then
+         if (bgj_restricted(rtdb)) then
             scftype = 'RHF'
          else
             scftype = 'UHF'
@@ -91,7 +91,7 @@ subroutine cphf_solve(rtdb)
      &       UNKNOWN_ERR)
       nbases = 1
       bases(1) = basis
-      if (bgj_have_j_fit()) then
+      if (bgj_have_j_fit(rtdb)) then
         nbases = 2
         bases(2) = bgj_CD_bas_han()
         write(*,*) 'cphf_solve: bases',bases

src/cphf/cphf_solve2.F

@@ -96,7 +96,7 @@ subroutine cphf_solve2(rtdb)
       nbases = 1
       bases(1) = basis
 #if 0
-      if (bgj_have_j_fit()) then
+      if (bgj_have_j_fit(rtdb)) then
         nbases = 2
         bases(2) = bgj_CD_bas_han()
         if (ga_nodeid().eq.0) then

src/cphf/cphf_solve3.F

@@ -149,7 +149,7 @@ subroutine cphf_solve3(rtdb, omega, lifetime, gamwidth)
      $        call errquit('cphf_solve3: no DFT MO vectors',0, RTDB_ERR)
 #if 0
 c !!! Hack to set up scftype properly for later code !!! ? need this?
-         if (bgj_restricted()) then
+         if (bgj_restricted(rtdb)) then
             scftype = 'RHF'
          else
             scftype = 'UHF'
@@ -164,7 +164,7 @@ subroutine cphf_solve3(rtdb, omega, lifetime, gamwidth)
      &       UNKNOWN_ERR)
       nbases = 1
       bases(1) = basis
-      if (bgj_have_j_fit()) then
+      if (bgj_have_j_fit(rtdb)) then
         nbases = 2
         bases(2) = bgj_CD_bas_han()
         write(*,*) 'cphf_solve3: bases',bases

src/cphf/cphf_solve4.F

@@ -77,7 +77,7 @@ subroutine cphf_solve4(rtdb)
      $        call errquit('cphf_solve: no DFT MO vectors',0, RTDB_ERR)
 #if 0
 c !!! Hack to set up scftype properly for later code !!! ? need this?
-         if (bgj_restricted()) then
+         if (bgj_restricted(rtdb)) then
             scftype = 'RHF'
          else
             scftype = 'UHF'
@@ -92,7 +92,7 @@ subroutine cphf_solve4(rtdb)
      &       UNKNOWN_ERR)
       nbases = 1
       bases(1) = basis
-      if (bgj_have_j_fit()) then
+      if (bgj_have_j_fit(rtdb)) then
         nbases = 2
         bases(2) = bgj_CD_bas_han()
         write(*,*) 'cphf_solve4: bases',bases

src/ddscf/fock_2e.F

@@ -1,4 +1,4 @@
-      subroutine fock_2e( geom1, ao_basis, mfock, jfac, kfac,
+      subroutine fock_2e( rtdb, geom1, ao_basis, mfock, jfac, kfac,
      $                    tol2e1, oskel1, vg_dens, vg_fock, asym )
 c$Id$
 c
@@ -19,6 +19,7 @@ subroutine fock_2e( geom1, ao_basis, mfock, jfac, kfac,
 #include "util.fh"
 #include "stdio.fh"
 
+      integer rtdb
       integer geom1, ao_basis        ! [input] parameter handles
       integer nfock,mfock            ! [input] number of Fock matrices
       double precision tol2e1        ! [input] integral selection threshold
@@ -41,7 +42,7 @@ subroutine fock_2e( geom1, ao_basis, mfock, jfac, kfac,
       logical xc_active, jfit,oprint,odft,cphf_uhf
       oprint= util_print('fock_2e',print_debug)
       cphf_uhf = .false.
-      if (.not. rtdb_get(bgj_get_rtdb_handle(), 
+      if (.not. rtdb_get(rtdb, 
      &     'cphf_solve:cphf_uhf', mt_log, 1, cphf_uhf)) then
          cphf_uhf = .false.
       endif
@@ -59,23 +60,23 @@ subroutine fock_2e( geom1, ao_basis, mfock, jfac, kfac,
 c
 c     See if XC is active for this call
 c
-      if (.not.rtdb_get(bgj_get_rtdb_handle(),'bgj:xc_active',MT_LOG,1,
+      if (.not.rtdb_get(rtdb,'bgj:xc_active',MT_LOG,1,
      &     xc_active)) xc_active = .false.
       if (oprint)
      &   write(luout,*)'*** fock_2e: xc_active ',xc_active
 c
 c     Set K contribution to be correct even for pure DFT and hybrids
 c
-      if (xc_active .and. bgj_have_xc()) then
-         call dscal(nfock, bgj_kfac(), kfac, 1)
+      if (xc_active .and. bgj_have_xc(rtdb)) then
+         call dscal(nfock, bgj_kfac(rtdb), kfac, 1)
          if (oprint)
      &      write(luout,*)'*** fock_2e: set kfac to ',
      &        (kfac(ifock),ifock=1,nfock)
       endif
 c
 c     Determine whether J fitting is involved
 c
-      jfit = bgj_have_j_fit().and.(jfac(1).ne.0d0)
+      jfit = bgj_have_j_fit(rtdb).and.(jfac(1).ne.0d0)
 c !!! Edo commented this out before since it caused bugs - ask him about
 c !!! it when it is necessary to turn back on for J fit CPKS
       jfit = .false.
@@ -90,8 +91,8 @@ subroutine fock_2e( geom1, ao_basis, mfock, jfac, kfac,
      &      write(luout,*)'*** fock_2e: jfit case - original jfac = ',
      &        jfac(1)
 cold         jfac(1)=0d0
-cold         kfac(1)=-0.5d0*bgj_kfac()
-         call ao_fock_2e( geom1, ao_basis, 1, 0d0, -0.5d0*bgj_kfac(),
+cold         kfac(1)=-0.5d0*bgj_kfac(rtdb)
+         call ao_fock_2e(geom1, ao_basis, 1, 0d0, -0.5d0*bgj_kfac(rtdb),
      $        tol2e1, oskel1, vg_dens, vg_fock, asym )
 
          call fock_j_fit(nfock, vg_dens, vg_fockc)
@@ -137,15 +138,15 @@ subroutine fock_2e( geom1, ao_basis, mfock, jfac, kfac,
 c !!! Exit if no XC part to do - this can be changed to just do XC
 c !!! part and print it out for comparison by removing xc_active
 c !!! from the test below
-      if (.not. (xc_active .and. bgj_have_xc())) then
+      if (.not. (xc_active .and. bgj_have_xc(rtdb))) then
 c         write(*,*)'*** fock_2e: no xc, returning'
          return
       endif
       if(nfock.ne.1.and.nfock.ne.4) then
          if(ga_nodeid().eq.0) then
             write(0,*) ' WARNING: likely restart: nfock = ',nfock
          endif
-         if (.not.rtdb_put(bgj_get_rtdb_handle(),'bgj:xc_active',
+         if (.not.rtdb_put(rtdb,'bgj:xc_active',
      M        MT_LOG,1,.false.)) 
      C        call errquit(' fock2e: bgjrtdbput failed',0,0)
          return
@@ -183,7 +184,7 @@ subroutine fock_2e( geom1, ao_basis, mfock, jfac, kfac,
       Exc(1) = 0.0d0
       Exc(2) = 0.0d0
       nExc = 1
-      call fock_xc(geom1, nbf,ao_basis,
+      call fock_xc(rtdb,geom1, nbf,ao_basis,
      &             nfock, vg_dens, g_xc,Exc,nExc,.false.)
       ffac=1d0
       if (xc_active) then

src/ddscf/fock_2e_cam.F

@@ -1,4 +1,4 @@
-      subroutine fock_2e_cam( geom1, ao_basis, mfock, jfac, kfac,
+      subroutine fock_2e_cam( rtdb, geom1, ao_basis, mfock, jfac, kfac,
      $           tol2e1, oskel1, vg_dens, vg_fock, asym, doxc )
 c$Id$
 c
@@ -42,13 +42,10 @@ subroutine fock_2e_cam( geom1, ao_basis, mfock, jfac, kfac,
 c
       external ga_create_atom_blocked
       logical xc_active, jfit,oprint,odft,cphf_uhf
-c
-c     Get a handle to the rtdb
-      rtdb = bgj_get_rtdb_handle()
 c
       oprint= util_print('fock_2e',print_debug)
       cphf_uhf = .false.
-      if (.not. rtdb_get(bgj_get_rtdb_handle(), 
+      if (.not. rtdb_get(rtdb, 
      &     'cphf_solve:cphf_uhf', mt_log, 1, cphf_uhf)) then
          cphf_uhf = .false.
       endif
@@ -66,23 +63,23 @@ subroutine fock_2e_cam( geom1, ao_basis, mfock, jfac, kfac,
 c
 c     See if XC is active for this call
 c
-      if (.not.rtdb_get(bgj_get_rtdb_handle(),'bgj:xc_active',MT_LOG,1,
+      if (.not.rtdb_get(rtdb,'bgj:xc_active',MT_LOG,1,
      &     xc_active)) xc_active = .false.
       if (oprint)
      &   write(luout,*)'*** fock_2e: xc_active ',xc_active
 c
 c     Set K contribution to be correct even for pure DFT and hybrids
 c
-      if (xc_active .and. bgj_have_xc()) then
-         call dscal(nfock, bgj_kfac(), kfac, 1)
+      if (xc_active .and. bgj_have_xc(rtdb)) then
+         call dscal(nfock, bgj_kfac(rtdb), kfac, 1)
          if (oprint)
      &      write(luout,*)'*** fock_2e: set kfac to ',
      &        (kfac(ifock),ifock=1,nfock)
       endif
 c
 c     Determine whether J fitting is involved
 c
-      jfit = bgj_have_j_fit().and.(jfac(1).ne.0d0)
+      jfit = bgj_have_j_fit(rtdb).and.(jfac(1).ne.0d0)
 c !!! Edo commented this out before since it caused bugs - ask him about
 c !!! it when it is necessary to turn back on for J fit CPKS
       jfit = .false.
@@ -97,8 +94,8 @@ subroutine fock_2e_cam( geom1, ao_basis, mfock, jfac, kfac,
      &      write(luout,*)'*** fock_2e: jfit case - original jfac = ',
      &        jfac(1)
 cold         jfac(1)=0d0
-cold         kfac(1)=-0.5d0*bgj_kfac()
-         call ao_fock_2e( geom1, ao_basis, 1, 0d0, -0.5d0*bgj_kfac(),
+cold         kfac(1)=-0.5d0*bgj_kfac(rtdb)
+         call ao_fock_2e(geom1, ao_basis, 1, 0d0, -0.5d0*bgj_kfac(rtdb),
      $        tol2e1, oskel1, vg_dens, vg_fock, asym )
 
          call fock_j_fit(nfock, vg_dens, vg_fockc)
@@ -135,7 +132,7 @@ subroutine fock_2e_cam( geom1, ao_basis, mfock, jfac, kfac,
      $          tol2e1, oskel1, vg_dens, vg_fock,asym )
         endif
       else
-         if (bgj_have_j_fit()) call fock_force_direct(rtdb)
+         if (bgj_have_j_fit(rtdb)) call fock_force_direct(rtdb)
 
          call ao_fock_2e( geom1, ao_basis, nfock, jfac, kfac,
      $                    tol2e1, oskel1, vg_dens, vg_fock, asym )
@@ -146,15 +143,15 @@ subroutine fock_2e_cam( geom1, ao_basis, mfock, jfac, kfac,
 c !!! Exit if no XC part to do - this can be changed to just do XC
 c !!! part and print it out for comparison by removing xc_active
 c !!! from the test below
-      if (.not. (xc_active .and. bgj_have_xc())) then
+      if (.not. (xc_active .and. bgj_have_xc(rtdb))) then
 c         write(*,*)'*** fock_2e: no xc, returning'
          return
       endif
       if(nfock.ne.1.and.nfock.ne.4) then
          if(ga_nodeid().eq.0) then
             write(0,*) ' WARNING: likely restart: nfock = ',nfock
          endif
-         if (.not.rtdb_put(bgj_get_rtdb_handle(),'bgj:xc_active',
+         if (.not.rtdb_put(rtdb,'bgj:xc_active',
      M        MT_LOG,1,.false.)) 
      C        call errquit(' fock2e: bgjrtdbput failed',0,0)
          return
@@ -193,7 +190,7 @@ subroutine fock_2e_cam( geom1, ao_basis, mfock, jfac, kfac,
       Exc(2) = 0.0d0
       nExc = 1
       if (doxc) then
-      call fock_xc(geom1, nbf,ao_basis,
+      call fock_xc(rtdb,geom1, nbf,ao_basis,
      &             nfock, vg_dens, g_xc,Exc,nExc,.false.)
       ffac=1d0
       if (xc_active) then

src/ddscf/fock_j_fit.F

@@ -62,20 +62,19 @@ subroutine fock_j_fit_gen(nmat, g_dens, g_j, Tvec, have_Tvec,
       double precision ddum
       logical IOLGC, old_incore, old_direct
 c
-      if (bgj_print() .gt. 0)
+      if (bgj_print(rtdb) .gt. 0)
      &     write(*,*)'--------Entered fock_j_fit-------------'
-      rtdb = bgj_get_rtdb_handle()
 c      write(*,*)'nmat rtdb geom',nmat,rtdb,geom
 c
 c     Get fock_j_fit variables
 c
       if (.not. rtdb_get(rtdb, 'fock_j:derfit', mt_log, 1,
      &     derfit)) then
-         if (bgj_print() .gt. 0)
+         if (bgj_print(rtdb) .gt. 0)
      &      write(*,*)'*** fock_j_fit: derfit not set: setting to false'
          derfit = .false.       !!!
       endif
-      if (bgj_print() .gt. 0)
+      if (bgj_print(rtdb) .gt. 0)
      &     write(*,*)'fock_j_fit: derfit =',derfit
 c
 c     Initialize J matrices
@@ -143,7 +142,7 @@ subroutine fock_j_fit_gen(nmat, g_dens, g_j, Tvec, have_Tvec,
       direct = old_direct
       incore = old_incore
 c
-      if (bgj_print() .gt. 0) then
+      if (bgj_print(rtdb) .gt. 0) then
          do imat = 1, nmat
             write(*,*)'Fitted J matrix',imat
             call ga_print(g_j(imat))

src/ddscf/fock_xc.F

@@ -1,7 +1,7 @@
 C>
 C> \brief Wrapper for AO-basis XC matrices without fitting
 C>
-      subroutine fock_xc(geom, nbf_ao, ao_bas_han,
+      subroutine fock_xc(rtdb, geom, nbf_ao, ao_bas_han,
      &   nfock, g_dens, g_xc, Exc, nExc, l3d)
 c
 c     $Id$
@@ -87,7 +87,6 @@ subroutine fock_xc(geom, nbf_ao, ao_bas_han,
 c
       ndens=1
       if(xc_getipol().eq.2) ndens=2
-      rtdb = bgj_get_rtdb_handle()
       alo(2) = 1
       ahi(2) = nbf_ao
       alo(3) = 1
@@ -234,7 +233,7 @@ subroutine fock_xc(geom, nbf_ao, ao_bas_han,
          if(ndens.eq.2) int_mb(igxcd+1)=g_dens(3)
       endif
       if (calc_type .eq. 2) then
-         if (.not. bgj_get_scf_dens(g_dens_scf))
+         if (.not. bgj_get_scf_dens(rtdb,g_dens_scf))
      &        call errquit('fock_xc: cant get scf density handles',0,
      &       UNKNOWN_ERR)
          alo(1) = 1

src/ddscf/rhf_dens_mo.F

@@ -85,7 +85,7 @@ subroutine rhf_dens_to_mo(rtdb, geom, basis, nelec, nbf, nmo,
       call int_1e_ga(basis, basis, g_hcore,'potential', oskel)
 c     
       if (nelec.gt.1 .and. abs(dens_norm) .gt. 1d-2)
-     $     call rhf_fock_2e(geom, basis, g_dens, g_fock, tol2e,
+     $     call rhf_fock_2e(rtdb, geom, basis, g_dens, g_fock, tol2e,
      &     .true., .true., oskel)
 c     
 c     Compute contributions to the energy and symmetrize

src/ddscf/rhf_fock_2e.F

@@ -1,4 +1,4 @@
-      subroutine rhf_fock_2e( geom, basis, g_dens, g_fock, tol2e,
+      subroutine rhf_fock_2e( rtdb, geom, basis, g_dens, g_fock, tol2e,
      $     ocoul, oexch, oskel )
 *
 * $Id$
@@ -12,7 +12,7 @@ subroutine rhf_fock_2e( geom, basis, g_dens, g_fock, tol2e,
 c     Now also called from moints_full to generate an effective
 c     frozen-core hamiltonian
 c     
-      integer geom, basis
+      integer rtdb, geom, basis
       integer g_dens, g_fock    ! [input] handles
       double precision tol2e    ! [input] integral selection threshold
       logical ocoul, oexch      ! [input] compute coulomb/exchange 
@@ -38,7 +38,7 @@ subroutine rhf_fock_2e( geom, basis, g_dens, g_fock, tol2e,
       vg_dens(1) = g_dens
       vg_fock(1) = g_fock
 
-      call fock_2e( geom, basis, nfock, jfac, kfac, tol2e, oskel,
+      call fock_2e( rtdb, geom, basis, nfock, jfac, kfac, tol2e, oskel,
      $     vg_dens, vg_fock, .false.)
 
       return

src/ddscf/rohf_fock.F

@@ -229,7 +229,7 @@ subroutine rohf_fock(rtdb, geom, basis, nclosed, nopen, nmo,
       if ( nopen .eq. 0 ) call riscf_dovecs(.true., g_vecs)
 c
       if (.not.cam_exch) then
-        call fock_2e(geom, basis, nfock, jfactor, kfactor,
+        call fock_2e(rtdb, geom, basis, nfock, jfactor, kfactor,
      $     tol2e, oskel, iv_dens, iv_fock, .false.)
       else 
 c       for attenuated calculations