      SUBROUTINE new4ind_N5_osu(rtdb,d_v2,
     1                                kax_v2_alpha_offset,
     1                                size_2e)
C     $Id: new_ga4ind_N5.F 25748 2014-06-08 07:53:05Z d3y133 $
C     This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
C     Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
C     t ( p1 p2 h3 h4 )_t
      IMPLICIT NONE
#include "rtdb.fh"
#include "global.fh"
#include "mafdecls.fh"
#include "stdio.fh"
#include "util.fh"
#include "bas.fh"
#include "schwarz.fh"
#include "sym.fh"
#include "sf.fh"
#include "errquit.fh"
#include "tce.fh"
#include "tce_main.fh"
c
c written by K. Kowalski
c
c
c     max. number of p2 groups =200
c
c
      integer rtdb                 ! Run-time database
      integer d_v2                 ! MO integrals
      integer kax_v2_alpha_offset  ! MO integrals offset
      integer size_2e              ! 2e file size
      integer iter
      INTEGER azone1,azone2,azone3,azone4
      INTEGER g1b,g2b,g3b,g4b
      INTEGER igi1,igi2,igi3,igi4
      INTEGER ii,i,j,k,l,N,ipos1,ipos2,ipos3,ipos4
      INTEGER size_4a,l_4a,k_4a
      INTEGER size_aaaa, size_full_length
      INTEGER d_agaa,size_agaa,l_offset_agaa,k_offset_agaa
      INTEGER d_ggaa,size_ggaa,l_offset_ggaa,k_offset_ggaa
      INTEGER d_ggga,size_ggga,l_offset_ggga,k_offset_ggga
      INTEGER l_loc_aaaa,k_loc_aaaa,size_loc_aaaa
      INTEGER l_loc_agaa,k_loc_agaa,size_loc_agaa
      INTEGER d_agaa_f,size_agaa_f,l_offset_agaa_f,k_offset_agaa_f
      INTEGER l_loc_agaa_f,k_loc_agaa_f,size_loc_agaa_f
      INTEGER l_loc_ggaa,k_loc_ggaa,size_loc_ggaa
      INTEGER l_loc_ggga,k_loc_ggga,size_loc_ggga
      INTEGER l_loc_gggg,k_loc_gggg,size_loc_gggg
c
      integer key_aaaa,key_agaa,key_ggaa,key_ggga
      integer offset_aaaa,offset_agaa,offset_ggaa,offset_ggga
      integer key_gggg,offset_gggg
c
      integer max_size_temp,sumx
c
      integer tot_azone1_sh,tot_azone2_sh
      integer tot_azone3_sh,tot_azone4_sh
      integer tot_zone(1000)  !it was d.prec.
c
      integer iha,ihb !number of corr. alpha, beta holes
      integer ipa,ipb !number of corr. alpha, beta particles
c
      integer mu,nu,rho,sigma
      integer mu_lo,mu_hi
      integer nu_lo,nu_hi
      integer rho_lo,rho_hi
      integer sigma_lo,sigma_hi
      integer mu_range
      integer nu_range
      integer rho_range
      integer sigma_range
      integer mu1,nu1,rho1,sigma1
      integer shift_mu,shift_nu
      integer shift_rho,shift_sigma
      integer work1,work2          ! Work array sizes
      integer l_work1,k_work1      ! Work array 1
      integer l_work2,k_work2      ! Work array 2
      integer imu1,inu1,irho1,isigma1
c
      integer l_movecs_orb,k_movecs_orb
c
      integer l_integral,l_coeff
      integer k_integral,k_coeff
      integer l_aux,k_aux,size_aux
      integer size_ic,size_icc,size_integral,size_coeff,max_na
c
      integer INDEX_PAIR,icol,irow
c
ccx      double precision tot_zone(1000)
c 
      integer l_4af_offset,k_4af_offset,d_4af
      integer sf_chunk,request
      integer key_4af,offset_4af,size_4af
      character*255 filename
c 
      logical parallel
c
      INTEGER length
      INTEGER next
      INTEGER nprocs
      INTEGER count
      integer nnn, ilo,ihi,jlo,jhi
      integer nxtask
      external nxtask
      logical nodezero
      logical idiskl
      logical compute_gggg
c
c
      nodezero=(ga_nodeid().eq.0)
c
c
      max_size_temp=imaxsize**4
c
      do ii=1,1000
       tot_zone(ii)=0.0d0
      enddo
      if(atpart.gt.1000)
     &  call errquit('tce_zones: atpart too big',1,MA_ERR)
      sumx=0
      do ii=1,atpart
       tot_zone(ii)=sumx
       sumx=sumx+nalength(ii)
      enddo
c
c
c this module is called only if intorb = .true.
c N is the number of correlated orbitals
        N = nmo(1) - nfc(1) - nfv(1)
        iha = nocc(1)-nfc(1)
        ihb = nocc(ipol)-nfc(ipol)
        ipa = nmo(1)-nocc(1)-nfv(1)
        ipb = nmo(ipol)-nocc(ipol)-nfv(ipol)
c
c     pre-compute size_loc_agaa & size_loc_aaaa

      size_full_length=0
      DO azone3 =1,atpart
         size_full_length = size_full_length+nalength(azone3)
      enddo
c
      DO azone2 = 1,atpart
      size_loc_agaa=0
      size_loc_aaaa=0
      size_loc_agaa_f=0
        DO azone1 = 1,atpart
            nnn=nalength(azone1)*nalength(azone2)
            DO g3b = 1,noa+nva
                  size_loc_agaa_f= max(size_loc_agaa_f,
     V                 nnn*
     1                 size_full_length*int_mb(k_range_alpha+g3b-1))
               DO azone3 = 1,atpart
                  size_loc_agaa= max(size_loc_agaa,
     V                 nnn*
     1                 nalength(azone3)*int_mb(k_range_alpha+g3b-1))
                  DO azone4=1,atpart
                     
                     size_loc_aaaa=max(size_loc_aaaa,
     V                    nnn*
     1                    nalength(azone3)*nalength(azone4))
                  enddo
               enddo
            enddo
         enddo
c      enddo
c     pre-compute size_loc_ggaa & size_loc_ggga
c
      size_loc_ggaa=0
      size_loc_ggga=0
c      DO azone2 = 1,atpart
         DO g2b = 1,noa+nva
            DO g3b = 1,noa+nva
               DO g4b = g3b,noa+nva
                  size_loc_ggga= max(size_loc_ggga,
     N                 nalength(azone2)*int_mb(k_range_alpha+g2b-1)*
     1 int_mb(k_range_alpha+g3b-1)*int_mb(k_range_alpha+g4b-1))

                  DO azone1=1,atpart ! azone1
                     size_loc_ggaa=max(size_loc_ggaa,
     N                    nalength(azone1)*nalength(azone2)*
     1  int_mb(k_range_alpha+g3b-1)*int_mb(k_range_alpha+g4b-1))

                  enddo
               enddo
            enddo
         enddo
c      enddo
c *** debug ***
c        if(nodezero) then
c         write(6,*) 'in step1'
c         call util_flush(6)
c        endif
c *************
c
c     Offset for 4a file
c
      sf_chunk=(imaxsize)**4
      call tce_4a_offset_ns_osu(l_4af_offset,k_4af_offset,size_4af,
     &    azone2)

      call tce_gacreate_sloc_osu(d_4af, size_4af, 
     S     size_loc_aaaa, 'd_4af')
ccx      call ga_zero(d_4af)
c *** debug ***
c        if(nodezero) then
c         write(6,*) 'in step2'
c         call util_flush(6)
c        endif
c *************
c
c alpha orbitals only
c
      if (.not.ma_push_get(mt_dbl,nbf*(iha+ipa)
     1  ,"sorted MO coeffs",
     2  l_movecs_orb,k_movecs_orb))
     3  call errquit('tce_mo2e_zone: MA problem 1',0,
     2    BASIS_ERR)
c      call dcopy(nbf*(iha+ipa),0.0d0, 0,dbl_mb(k_movecs_orb), 1)
      call dcopy(nbf*iha,
     c     dbl_mb(k_movecs_sorted),1,
     C     dbl_mb(k_movecs_orb),1)
      call dcopy(nbf*ipa,
     c     dbl_mb(k_movecs_sorted+(iha+ihb)*nbf),1,
     C     dbl_mb(k_movecs_orb+iha*nbf),1)
c
c
      call int_mem_2e4c(work1,work2)
      if (.not.ma_push_get(mt_dbl,work1,'work1',l_work1,k_work1))
     1  call errquit('tce_ao2e: MA problem work1',0,MA_ERR)
      if (.not.ma_push_get(mt_dbl,work2,'work2',l_work2,k_work2))
     1  call errquit('tce_ao2e: MA problem work2',1,MA_ERR)
c
c
c *** debug ***
c        if(nodezero) then
c         write(6,*) 'in step3'
c         call util_flush(6)
c        endif
c *************
c 
c 4af file formed here
c
c
c       if(nodezero) then
c          write(6,'(A,F20.2,A)') ' starting step 0 at ',
c     C         util_wallsec(), ' secs '
c          call util_flush(6)
c       endif
      nprocs = GA_NNODES()
      count = 0
      call  ga_distribution(d_4af,ga_nodeid(),ilo,ihi,jlo,jhi)
      if(jlo.gt.0) then
cold      next = NXTASK(nprocs, 1)
      DO azone1 = 1,atpart      !nu
c      DO azone2 = 1,atpart !mu
      DO azone3 = 1,atpart      !sigma
      DO azone4 = azone3,atpart !rho
         size_4a = nalength(azone1)*nalength(azone2)*
     1            nalength(azone3)*nalength(azone4)
         key_4af=azone4 - 1 + atpart * (azone3 - 1 +
     &          atpart * (azone2 - 1 + atpart * (azone1 - 1)))
        call tce_hash(int_mb(k_4af_offset),key_4af,offset_4af)
       if((ilo.le.offset_4af+size_4a).and.
     A      (ihi.ge.offset_4af+1)) then
cold      IF (next.eq.count) THEN
c ---------------------------
        if(.not.ma_push_get(mt_dbl,size_4a,'4a',l_4a,k_4a))
     1     call errquit('tce_4af_zones1: MA problem',0,MA_ERR)
        call dfill(size_4a, 0.0d0, dbl_mb(k_4a), 1)
         shift_mu = 0
         do mu    = a2length(azone2)+1,a2length(azone2+1)
            if (.not.bas_cn2bfr(ao_bas_han,mu,mu_lo,mu_hi))
     1      call errquit('tce_ao2e: basis fn range problem 1',0,
     2      BASIS_ERR)
            mu_range = mu_hi - mu_lo + 1
         shift_nu = 0
         do nu    = a2length(azone1)+1,a2length(azone1+1)
            if (.not.bas_cn2bfr(ao_bas_han,nu,nu_lo,nu_hi))
     1      call errquit('tce_ao2e: basis fn range problem 1',0,
     2      BASIS_ERR)
            nu_range = nu_hi - nu_lo + 1
         shift_rho = 0
         do rho   = a2length(azone4)+1,a2length(azone4+1)
            if (.not.bas_cn2bfr(ao_bas_han,rho,rho_lo,rho_hi))
     1      call errquit('tce_ao2e: basis fn range problem 1',0,
     2      BASIS_ERR)
            rho_range = rho_hi - rho_lo + 1
         shift_sigma = 0
         do sigma = a2length(azone3)+1,a2length(azone3+1)
            if (.not.bas_cn2bfr(ao_bas_han,sigma,sigma_lo,sigma_hi))
     1      call errquit('tce_ao2e: basis fn range problem 1',0,
     2      BASIS_ERR)
            sigma_range = sigma_hi - sigma_lo + 1
            if (schwarz_shell(rho,sigma)*schwarz_shell(mu,nu)
     1          .ge. tol2e) then
            call int_2e4c(ao_bas_han,mu,nu,ao_bas_han,rho,sigma,
     1           work2,dbl_mb(k_work2),work1,dbl_mb(k_work1))
c
            i=0
            do mu1     = 1,mu_range
               imu1=mu1+shift_mu
               do nu1     = 1,nu_range
                  inu1=nu1+shift_nu
                  do rho1    = 1,rho_range
                     irho1=rho1+shift_rho
                     do sigma1  = 1,sigma_range
                        i=i+1
                        isigma1=sigma1+shift_sigma
c     (isigma1,irho1|inu1, imu1)
                        ipos1=(((imu1-1)*nalength(azone1)+inu1-1)*
     1                       nalength(azone4)+irho1-1)*nalength(azone3)
     2                       +isigma1
                        dbl_mb(k_4a+ipos1-1)=dbl_mb(k_work1+i-1)
                     enddo
                  enddo
               enddo
            enddo
            end if !schwarz  screening
         shift_sigma = shift_sigma + sigma_range
         enddo !sigma
         shift_rho   = shift_rho + rho_range
         enddo !rho
         shift_nu    = shift_nu + nu_range
         enddo !nu
         shift_mu    = shift_mu + mu_range
         enddo !mu
c
c fixing offsets and sf_writing
        call ga_put(d_4af,offset_4af+1,offset_4af+size_4a,1,1,
     1    dbl_mb(k_4a),1)
c       write(6,'(A,I5,A,I5,A,I5,A,I5)') "azone4 : ",azone4," azone3 : ",
c     C  azone3,"azone2 : ",azone2," azone1 : ",azone1
c       write(6,'(A,I5,A,I5,A,I5,A,I5)') "azone4 : ",azone4," azone3 : ",
c     C  azone3,"azone2 : ",azone2," azone1 : ",azone1
c        do iter = 1,20
c          write(6,'(F5.2,A,$)') dbl_mb(k_4a+iter),' , '
c          call util_flush(6)
c        enddo
c closing l_4a file
        if (.not.ma_pop_stack(l_4a))
     1   call errquit('tcc_mo2e_4af2: l_4a',15,MA_ERR)
c ---------------------------
cold      next = NXTASK(nprocs, 1)
      END IF
      count = count + 1
      ENDDO !azone4
      ENDDO !azone3
c      ENDDO !azone2
      ENDDO !azone1
      endif
c
cold      next = NXTASK(-nprocs, 1)
      call ga_sync()
c
      call tce_offset_ggaa_osu(l_offset_ggaa,k_offset_ggaa,
     &     size_ggaa,azone2)
c     open offset l_ggga
c      call tce_offset_ggga(l_offset_ggga,k_offset_ggga,size_ggga,azone2)
c
      max_na=0
      do i=1,atpart
       if(nalength(i).gt.max_na) max_na=nalength(i)
      enddo
      size_icc=tile_dim*max_na
c
       if (.not.ma_push_get(mt_dbl,size_icc,'l_coeff',
     1  l_coeff,k_coeff))
     1  call errquit('tce_4s: MA problem l_coeff',0,MA_ERR)
c
c     step1
c
c *** debug ***
c        if(nodezero) then
c         write(6,*) 'in step4'
c         call util_flush(6)
c        endif
c *************
c
c     open ga d_agaa (size_agaa)
c       if(nodezero) then
c          write(6,'(A,F20.2,A)') ' starting step 1+2 at ',
c     C         util_wallsec(), ' secs '
c          call util_flush(6)
c       endif
      call tce_gacreate_sloc_osu(d_ggaa, size_ggaa,
     S     size_loc_ggaa, 'd_ggaa')
c
      nprocs = GA_NNODES()
      count = 0
cold      next = NXTASK(nprocs, 1)
      call  ga_distribution(d_ggaa,ga_nodeid(),ilo,ihi,jlo,jhi)
c       
      if(jlo.gt.0) then
       if (.not.ma_push_get(mt_dbl,size_loc_agaa_f,'loc_agaa_f',
     1  l_loc_agaa_f,k_loc_agaa_f))
     1  call errquit('step0:0',0,MA_ERR)
       if (.not.ma_push_get(mt_dbl,size_loc_agaa,'loc_agaa',
     1  l_loc_agaa,k_loc_agaa))
     1  call errquit('step1:1',0,MA_ERR)
        if (.not.ma_push_get(mt_dbl,size_loc_aaaa,'loc_aaaa',
     1  l_loc_aaaa,k_loc_aaaa))
     1  call errquit('step1:2',0,MA_ERR)
       if (.not.ma_push_get(mt_dbl,
     M   max(size_loc_aaaa,size_loc_agaa), 'auxaaaa',
     1   l_aux,k_aux))
     1   call errquit('step1:3',0,MA_ERR)
c
c *** debug ***
c        if(nodezero) then
c         write(6,*) 'in step5'
c         call util_flush(6)
c        endif
c *************
      nprocs = GA_NNODES()
      count = 0
cold      next = NXTASK(nprocs, 1)
c do parallel 
c      DO azone2 = 1,atpart
      DO azone1 = 1,atpart
      DO g3b = 1,noa+nva
c the above loops are loops fused for local memory enhancement
       g4b=g3b
       key_ggaa=g4b-1+(noa+nva)*(g3b-1+
     &  (noa+nva)*(azone1-1+atpart*(azone2-1)))
c       write(6,'(A,I5,A,I5,A,I5,A,I5)') "Before g3b : ",g3b," g4b : ",
c     C  g4b,"azone2 : ",azone2," azone1 : ",azone1       
c *** debug ***
c        if(nodezero) then
c         write(6,*) 'in step6'
c         call util_flush(6)
c        endif
c *************
       call tce_hash_n(int_mb(k_offset_ggaa),key_ggaa,offset_ggaa)
c *** debug ***
c        if(nodezero) then
c         write(6,*) 'in step6a'
c         call util_flush(6)
c        endif
c *************
       size_loc_ggaa= nalength(azone1)*nalength(azone2)*
     1 int_mb(k_range_alpha+g3b-1)*int_mb(k_range_alpha+g4b-1)
c       write(6,'(A,I5,A,I5,A,I5,A,I5)') "After g3b : ",g3b," g4b : ",
c     C  g4b,"azone2 : ",azone2," azone1 : ",azone1       
c the condition identifies, if this processor is responsible
c for computing ggaa[azone2,azone1,g3b,*] tiles. If it is responsible for 
c combuting g4b=1, then it is responsible for computing all g4b
       if((ilo.le.offset_ggaa+size_loc_ggaa).and.
     A      (ihi.ge.offset_ggaa+1)) then         
c
      call dfill(size_loc_agaa_f,0.0d0, dbl_mb(k_loc_agaa_f), 1)
      k_offset_agaa_f = 0;
      DO azone3 = 1,atpart
       size_loc_agaa= nalength(azone1)*nalength(azone2)*
     1 nalength(azone3)*int_mb(k_range_alpha+g3b-1)
       call dfill(size_loc_agaa,0.0d0, dbl_mb(k_loc_agaa), 1)
       DO azone4=1,atpart
        tot_azone4_sh=tot_zone(azone4)
        size_loc_aaaa=nalength(azone1)*nalength(azone2)*
     1  nalength(azone3)*nalength(azone4)
c --- C_[g3][azone4]
        i = 0
        do irho1 =  1,nalength(azone4)
        do igi3    =  1,int_mb(k_range_alpha+g3b-1)
        i = i+1
        ipos1=(int_mb(k_offset_alpha+g3b-1)+igi3-1)*nbf+tot_azone4_sh
     &       +irho1
        dbl_mb(k_coeff+i-1)=dbl_mb(k_movecs_orb+ipos1-1)
        enddo
        enddo
c ---
        if(azone4.le.azone3) then ! azone4 <= azone3
         key_aaaa=azone3-1+atpart*(azone4-1+
     &          atpart*(azone2-1+atpart*(azone1-1)))
         call tce_hash(int_mb(k_4af_offset),key_aaaa,offset_aaaa)
         call ga_get(d_4af,offset_aaaa+1,offset_aaaa+size_loc_aaaa,1,1,
     1    dbl_mb(k_loc_aaaa),1)
c now C([g3][azone4]) ([azone4]<=[azone3]|[azone1]<=[azone2]) 
        call  dgemm('N','N',
     1   int_mb(k_range_alpha+g3b-1),  !m
     1   nalength(azone1)*nalength(azone2)*nalength(azone3), !n
     1   nalength(azone4), !k
     1   1.0d0,
     1   dbl_mb(k_coeff),                    !A
     1   int_mb(k_range_alpha+g3b-1), !lda
     1   dbl_mb(k_loc_aaaa),nalength(azone4),          ! B,ldb
     1   1.0d0,dbl_mb(k_loc_agaa),
     1   int_mb(k_range_alpha+g3b-1))
        else    ! azone4 > azone3
         key_aaaa=azone4-1+atpart*(azone3-1+
     &          atpart*(azone2-1+atpart*(azone1-1)))
         call tce_hash(int_mb(k_4af_offset),key_aaaa,offset_aaaa)
         call ga_get(d_4af,offset_aaaa+1,offset_aaaa+size_loc_aaaa,1,1,
     1    dbl_mb(k_loc_aaaa),1)
         CALL TCE_SORT_4KG_(dbl_mb(k_loc_aaaa),dbl_mb(k_aux),
     &   nalength(azone3),nalength(azone4),
     &   nalength(azone1),nalength(azone2),
     &   2,1,3,4,1.0d0)
c C([g3][azone4]) ([azone4]>[azone3]|[azone1]<=[azone2])
        call  dgemm('N','N',
     1   int_mb(k_range_alpha+g3b-1),  !m
     1   nalength(azone1)*nalength(azone2)*nalength(azone3), !n
     1   nalength(azone4), !k
     1   1.0d0,
     1   dbl_mb(k_coeff),                    !A
     1   int_mb(k_range_alpha+g3b-1), !lda
     1   dbl_mb(k_aux),nalength(azone4),          ! B,ldb
     1   1.0d0,dbl_mb(k_loc_agaa),
     1   int_mb(k_range_alpha+g3b-1))
        end if  ! azone4 <= azone3
       ENDDO  !azone4
C TRANSPOSITION HERE (g3 az3|az1 <=az2) => (az3 g3|az1 <=az2)
       CALL TCE_SORT_4KG_(dbl_mb(k_loc_agaa),
     &  dbl_mb(k_loc_agaa_f+k_offset_agaa_f),
     &  int_mb(k_range_alpha+g3b-1),nalength(azone3),
     &  nalength(azone1),nalength(azone2),
     &  2,1,3,4,1.0d0)
c       write(6,'(A,I5,A,I5,A,I5,A,I5)') "g3b : ",g3b," azone3 : ",
c     C  azone3,"azone2 : ",azone2," azone1 : ",azone1
c       write(6,'(A,I5,A,I5,A,I5,A,I5)') "g3b : ",g3b," azone3 : ",
c     C  azone3,"azone2 : ",azone2," azone1 : ",azone1
c        do iter = 1,20
c          write(6,'(F5.2,A,$)') dbl_mb(k_aux+iter),' , '
c          call util_flush(6)
c        enddo
c
cold      next = NXTASK(nprocs, 1)
      k_offset_agaa_f = k_offset_agaa_f + size_loc_agaa
      count = count + 1
      ENDDO !azone3
c       if(nodezero) then
c          write(6,'(A,F20.2,A)') ' starting step 2 at ',
c     C         util_wallsec(), ' secs '
c          call util_flush(6)
c       endif
      
      DO g4b = g3b,noa+nva
cold      IF (next.eq.count) THEN
       key_ggaa=g4b-1+(noa+nva)*(g3b-1+
     &  (noa+nva)*(azone1-1+atpart*(azone2-1)))
c       write(6,'(A,I5,A,I5,A,I5,A,I5)') "g3b : ",g3b," g4b : ",
c     C  g4b,"azone2 : ",azone2," azone1 : ",azone1
c *** debug ***
c        if(nodezero) then
c         write(6,*) 'in step7'
c         call util_flush(6)
c        endif
c *************
       call tce_hash_n(int_mb(k_offset_ggaa),key_ggaa,offset_ggaa)
c *** debug ***
c        if(nodezero) then
c         write(6,*) 'in step7a'
c         call util_flush(6)
c        endif
c *************
c       write(6,'(A,F20.2,A)') ' Done ',
c     C      util_wallsec(), ' secs '
c       call util_flush(6)
       size_loc_ggaa= nalength(azone1)*nalength(azone2)*
     1 int_mb(k_range_alpha+g3b-1)*int_mb(k_range_alpha+g4b-1)
       if (.not.ma_push_get(mt_dbl,size_loc_ggaa,'loc_ggaa',
     1 l_loc_ggaa,k_loc_ggaa))
     1 call errquit('step2:1',0,MA_ERR)
      call dfill(size_loc_ggaa,0.0d0, dbl_mb(k_loc_ggaa), 1)
      k_offset_agaa_f = 0
c ----- Inner Loop for Computing ggaa ---------- c
       DO azone3=1,atpart
        tot_azone3_sh=tot_zone(azone3)
        size_loc_agaa=nalength(azone1)*nalength(azone2)*
     1  int_mb(k_range_alpha+g3b-1)*nalength(azone3)
c        if (.not.ma_push_get(mt_dbl,size_loc_agaa,'loc_agaa',
c     1  l_loc_agaa,k_loc_agaa))
c     1 call errquit('step2:2',0,MA_ERR)
c --- C_[g4][azone3]
        i = 0
        do isigma1 =  1,nalength(azone3)
        do igi4    =  1,int_mb(k_range_alpha+g4b-1)
        i = i+1
        ipos1=(int_mb(k_offset_alpha+g4b-1)+igi4-1)*nbf+tot_azone3_sh
     &       +isigma1
        dbl_mb(k_coeff+i-1)=dbl_mb(k_movecs_orb+ipos1-1)
        enddo
        enddo
c ---
c  C([g4b][az3]) * ([az3][g3b]|[az1][az2])
        call  dgemm('N','N',
     1   int_mb(k_range_alpha+g4b-1),  !m
     1   nalength(azone1)*nalength(azone2)*int_mb(k_range_alpha+g3b-1), !n
     1   nalength(azone3), !k
     1   1.0d0,
     1   dbl_mb(k_coeff),                    !A
     1   int_mb(k_range_alpha+g4b-1), !lda
     1   dbl_mb(k_loc_agaa_f+k_offset_agaa_f),nalength(azone3),          ! B,ldb
     1   1.0d0,dbl_mb(k_loc_ggaa),
     1   int_mb(k_range_alpha+g4b-1))
c
       k_offset_agaa_f = k_offset_agaa_f+size_loc_agaa 
       ENDDO
c     
       call ga_put(d_ggaa,offset_ggaa+1,offset_ggaa+size_loc_ggaa,1,1,
     1  dbl_mb(k_loc_ggaa),1)

       if (.not.ma_pop_stack(l_loc_ggaa))
     1 call errquit('g4ind:MA1',15,MA_ERR)

       ENDDO
c     
c       if((.not.nodezero).and.(g3b.eq.1)) then 
c          write(6,'(A)') "a"
c          write(6,'(A,I2,A,I2,A,I2,A,I2)') "g3b : ",g3b," g4b : ",
c     C         g4b,"azone2 : ",azone2," azone1 : ",azone1
cc     
c          do iter = 1,10
c             write(6,'(F5.2,A,$)') dbl_mb(k_loc_ggaa+iter),' , '
c             call util_flush(6)
c          enddo
c       end if
cold      next = NXTASK(nprocs, 1)
      END IF
      ENDDO !g3b
      ENDDO !azone1


c      ENDDO !azone2
        if (.not.ma_pop_stack(l_aux))
     1  call errquit('l_aux:MA1',15,MA_ERR)
        if (.not.ma_pop_stack(l_loc_aaaa))
     1  call errquit('l_loc_aaaa:MA1',15,MA_ERR)
       if (.not.ma_pop_stack(l_loc_agaa))
     1 call errquit('l_loc_agaa:MA1',15,MA_ERR)
       if (.not.ma_pop_stack(l_loc_agaa_f))
     1 call errquit('l_loc_agaa_f:MA1',15,MA_ERR)
cold      next = NXTASK(-nprocs, 1)

      endif
      call ga_sync()      
c     delete d_4af
      call deletefile(d_4af)

c
c     step3
c
c     open ga d_ggga (size_ggga)
c       if(nodezero) then
c          write(6,'(A,F20.2,A)') ' starting step 3 at ',
c     C         util_wallsec(), ' secs '
c          call util_flush(6)
c       endif
      call  ga_distribution(d_v2,ga_nodeid(),ilo,ihi,jlo,jhi)
      nprocs = GA_NNODES()
      count = 0
c      if(nodezero) then
c         write(6,'(A,I5,A,I5,A,I5,A,I5)')  
c     C        "High Low :",ilo,", ",
c     C        ihi,", ",jlo,", ",jhi
c      endif
      if(jlo.gt.0) then
      DO g3b = 1,noa+nva
      DO g4b = g3b,noa+nva
      compute_gggg=.FALSE.
c     Checking if these g3b and g4b are computed on this node
      DO g1b = 1,noa+nva
      DO g2b = g1b,noa+nva
       IF (int_mb(k_spin_alpha+g3b-1)+int_mb(k_spin_alpha+g4b-1).eq.
     &        int_mb(k_spin_alpha+g1b-1)
     &        +int_mb(k_spin_alpha+g2b-1)) THEN
          IF (ieor(int_mb(k_sym_alpha+g3b-1),
     &         ieor(int_mb(k_sym_alpha+g4b-1)
     &         ,ieor(int_mb(k_sym_alpha+g1b-1),
     &         int_mb(k_sym_alpha+g2b-1)))) .eq.
     &         irrep_v) THEN
             IROW=INDEX_PAIR(g4b,g3b)
             ICOL=INDEX_PAIR(g2b,g1b)
             IF(IROW.GE.ICOL) THEN
                size_loc_gggg= int_mb(k_range_alpha+g1b-1)*
     1          int_mb(k_range_alpha+g2b-1)*
     1          int_mb(k_range_alpha+g3b-1)*int_mb(k_range_alpha+g4b-1)
                key_gggg=g2b - 1 + (noa+nva) *
     &               (g1b - 1 + (noa+nva) * (g4b-
     &               1 + (noa+nva) * (g3b - 1)))
                call tce_hash_v2(int_mb(k_v2_alpha_offset),
     &               key_gggg,offset_gggg)
                if((ilo.le.offset_gggg+size_loc_gggg).and.
     A               (ihi.ge.offset_gggg+1)) then
                   compute_gggg = .TRUE.
c                   if(nodezero) then 
c                      write(6,'(I5,A,I5,A,I5,A,I5,A,I5,A,I5)')  
c     C                     ,g3b,", ",
c     C                     g4b,", ",g1b,", ",g2b,") : ",offset_gggg,
c     C                     ", ",offset_gggg+size_loc_gggg
c                      endif
c
                endif
                go to 10
             endif
          endif
       endif
      ENDDO
      ENDDO
 10   continue
c The tile must be computed at this node
      IF(compute_gggg) THEN
      DO g2b = 1,noa+nva
cold      IF (next.eq.count) THEN
       size_loc_ggga= nalength(azone2)*int_mb(k_range_alpha+g2b-1)*
     1 int_mb(k_range_alpha+g3b-1)*int_mb(k_range_alpha+g4b-1)
       if (.not.ma_push_get(mt_dbl,size_loc_ggga,'loc_ggga',
     1 l_loc_ggga,k_loc_ggga))
     1 call errquit('step3:1',0,MA_ERR)
       call dfill(size_loc_ggga,0.0d0, dbl_mb(k_loc_ggga), 1)
       DO azone1=1,atpart ! azone1
        tot_azone1_sh=tot_zone(azone1)
        size_loc_ggaa=nalength(azone1)*nalength(azone2)*
     1  int_mb(k_range_alpha+g3b-1)*int_mb(k_range_alpha+g4b-1) 
        if (.not.ma_push_get(mt_dbl,size_loc_ggaa,'loc_ggaa',
     1  l_loc_ggaa,k_loc_ggaa))
     1  call errquit('step3:2',0,MA_ERR)
c --- C_[azone1][g2]
        ii=0
        do igi2 =  1,int_mb(k_range_alpha+g2b-1)
        do inu1 =  1,nalength(azone1)
         ii = ii+1
         ipos1=(int_mb(k_offset_alpha+g2b-1)+igi2-1)*nbf+tot_azone1_sh
     &           +inu1
         dbl_mb(k_coeff+ii-1)=dbl_mb(k_movecs_orb+ipos1-1)
        enddo
        enddo
c ---
         size_loc_ggaa= nalength(azone1)*nalength(azone2)*
     1   int_mb(k_range_alpha+g3b-1)*int_mb(k_range_alpha+g4b-1)
         key_ggaa=g4b-1+(noa+nva)*(g3b-1+
     &   (noa+nva)*(azone1-1+atpart*(azone2-1)))
         call tce_hash(int_mb(k_offset_ggaa),key_ggaa,offset_ggaa)
         call ga_get(d_ggaa,offset_ggaa+1,offset_ggaa+size_loc_ggaa,1,1,
     1   dbl_mb(k_loc_ggaa),1)
         if (.not.ma_push_get(mt_dbl,size_loc_ggaa,'auxaaaa',
     1   l_aux,k_aux))
     1   call errquit('step3:3',0,MA_ERR)
         CALL TCE_SORT_4KG_(dbl_mb(k_loc_ggaa),dbl_mb(k_aux),
     &   int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1),
     &   nalength(azone1),nalength(azone2),
     &   1,2,4,3,1.0d0)
c ([g4]=>[g3]|[azone2]>[azone1])*C_([azone1][g2])
         call dgemm('N','N',
     1   int_mb(k_range_alpha+g4b-1)*int_mb(k_range_alpha+g3b-1)*
     1   nalength(azone2),                                        !m
     1   int_mb(k_range_alpha+g2b-1),                             !n
     3   nalength(azone1),                                        !k
     4   1.0d0,dbl_mb(k_aux),
     5   int_mb(k_range_alpha+g4b-1)*int_mb(k_range_alpha+g3b-1)*
     5   nalength(azone2),
     6   dbl_mb(k_coeff),nalength(azone1),
     7   1.0d0,dbl_mb(k_loc_ggga),
     8   int_mb(k_range_alpha+g4b-1)*int_mb(k_range_alpha+g3b-1)*
     8   nalength(azone2))
         if (.not.ma_pop_stack(l_aux))
     1   call errquit('g4ind:MA1',15,MA_ERR)
         if (.not.ma_pop_stack(l_loc_ggaa))
     1   call errquit('g4ind:MA1',15,MA_ERR)
       ENDDO ! azone1
C TRANSPOSITION HERE (g4=>g3|az2 g2) => (g4=>g3|g2 az2)
       if (.not.ma_push_get(mt_dbl,size_loc_ggga,'auxggga',
     1  l_aux,k_aux))
     1 call errquit('step3:4',0,MA_ERR)
       CALL TCE_SORT_4KG_(dbl_mb(k_loc_ggga),dbl_mb(k_aux),
     & int_mb(k_range_alpha+g4b-1),int_mb(k_range_alpha+g3b-1),
     & nalength(azone2),int_mb(k_range_alpha+g2b-1),
     & 1,2,4,3,1.0d0)
c       if(.NOT.nodezero) then
c          write(6,'(I5,A,I5,A,I5,A,I5)')  
c     C     ga_nodeid(), ": (",g3b,", ",
c     C         g4b,", ",g2b
c          write(6,'(A,F20.2,A)') ' starting step 4 at ',
c     C         util_wallsec(), ' secs '
c          call util_flush(6)
c       endif
      DO g1b = 1,g2b
       IF (int_mb(k_spin_alpha+g3b-1)+int_mb(k_spin_alpha+g4b-1).eq.
     & int_mb(k_spin_alpha+g1b-1)+int_mb(k_spin_alpha+g2b-1)) THEN
       IF (ieor(int_mb(k_sym_alpha+g3b-1),ieor(int_mb(k_sym_alpha+g4b-1)
     & ,ieor(int_mb(k_sym_alpha+g1b-1),int_mb(k_sym_alpha+g2b-1)))) .eq.
     & irrep_v) THEN
       IROW=INDEX_PAIR(g4b,g3b)
       ICOL=INDEX_PAIR(g2b,g1b)
       IF(IROW.GE.ICOL) THEN
       size_loc_gggg= int_mb(k_range_alpha+g1b-1)*
     1 int_mb(k_range_alpha+g2b-1)*
     1 int_mb(k_range_alpha+g3b-1)*int_mb(k_range_alpha+g4b-1)
       key_gggg=g2b - 1 + (noa+nva) *
     &(g1b - 1 + (noa+nva) * (g4b-
     & 1 + (noa+nva) * (g3b - 1)))
       call tce_hash_v2(int_mb(k_v2_alpha_offset),key_gggg,offset_gggg)
       if (.not.ma_push_get(mt_dbl,size_loc_gggg,'loc_gggg',
     1 l_loc_gggg,k_loc_gggg))
     1 call errquit('step4:1',0,MA_ERR)
      call dfill(size_loc_gggg,0.0d0, dbl_mb(k_loc_gggg), 1)
        tot_azone2_sh=tot_zone(azone2)
c --- C_[azone2][g1]
        ii=0
        do igi1 =  1,int_mb(k_range_alpha+g1b-1)
        do imu1 =  1,nalength(azone2)
         ii = ii+1
         ipos1=(int_mb(k_offset_alpha+g1b-1)+igi1-1)*nbf+tot_azone2_sh
     &           +imu1
         dbl_mb(k_coeff+ii-1)=dbl_mb(k_movecs_orb+ipos1-1)
        enddo
        enddo
c ---
        call dgemm('N','N',
     1   int_mb(k_range_alpha+g4b-1)*int_mb(k_range_alpha+g3b-1)*
     1   int_mb(k_range_alpha+g2b-1),                             !m
     1   int_mb(k_range_alpha+g1b-1),                             !n
     3   nalength(azone2),                                        !k
     4   1.0d0,dbl_mb(k_aux),
     5   int_mb(k_range_alpha+g4b-1)*int_mb(k_range_alpha+g3b-1)*
     5   int_mb(k_range_alpha+g2b-1),
     6   dbl_mb(k_coeff),nalength(azone2),
     7   1.0d0,dbl_mb(k_loc_gggg),
     8   int_mb(k_range_alpha+g4b-1)*int_mb(k_range_alpha+g3b-1)*
     8   int_mb(k_range_alpha+g2b-1))
c
       call ga_acc(d_v2,offset_gggg+1,offset_gggg+size_loc_gggg,1,1,
     1  dbl_mb(k_loc_gggg),1,1.0D0)
c       write(6,'(A,I5,A,I5,A,I5,A,I5)') "g3b : ",g3b," g4b : ",
c     C  g4b,"g2 : ",g2b," g1 : ",g1b
c       if(.NOT.nodezero) then
c          write(6,'(A,I5,A,I5,A,I5,A,I5,A,I5)') "Node ID : ", 
c     C     ga_nodeid(), "g3b : ",g3b," g4b : ",
c     C         g4b,"g2 : ",g2b," g1 : ",g1b
c          do iter = 1,20
c             write(6,'(F5.2,A,$)') dbl_mb(k_loc_gggg+iter),' , '
c             call util_flush(6)
c          enddo
c       endif
       if (.not.ma_pop_stack(l_loc_gggg))
     1 call errquit('g4ind:MA1',15,MA_ERR)
cold      next = NXTASK(nprocs, 1)
      END IF
      END IF
      END IF
      ENDDO !g1b
       if (.not.ma_pop_stack(l_aux))
     1  call errquit('g4ind:MA1',15,MA_ERR)
       if (.not.ma_pop_stack(l_loc_ggga))
     1 call errquit('g4ind:MA1',15,MA_ERR)
      ENDDO !g2b
      ENDIF
      ENDDO                     !g3b
      ENDDO !g4b
      endif
      call ga_sync()
c     delete l_coeff 
      if (.not.ma_pop_stack(l_coeff))
     1  call errquit('tcc_off_4a: MA problem',15,MA_ERR)
c     delete l_ggaa
      if (.not.ma_pop_stack(l_offset_ggaa))
     1  call errquit('ga4ind:ggaa',15,MA_ERR)
c     delete l_agaa    
c      if (.not.ma_pop_stack(l_offset_agaa))
c     1  call errquit('ga4ind:agaa',15,MA_ERR)
c
      if (.not.ma_pop_stack(l_work2))
     1  call errquit('tcc_ao2e: MA problem',14,MA_ERR)
c
      if (.not.ma_pop_stack(l_work1))
     1  call errquit('tcc_ao2e: MA problem',15,MA_ERR)
c
      if (.not.ma_pop_stack(l_movecs_orb))
     1  call errquit('tcc_ao2e: MA problem',15,MA_ERR)
c
      if (.not.ma_pop_stack(l_4af_offset))
     1  call errquit('ga4ind:4a',15,MA_ERR)
      
      ENDDO
c      if(nodezero) then
c          write(6,'(A,F20.2,A)') ' done step 4 at ',
c     C         util_wallsec(), ' secs '
c          call util_flush(6)
c       endif
c
c
c
      RETURN
      END
c
c
c
c
c
c
c
c
      SUBROUTINE tce_offset_agaa_osu(l_a_offset,k_a_offset,size,
     & azone2_para)
      IMPLICIT NONE
#include "global.fh"
#include "mafdecls.fh"
#include "sym.fh"
#include "errquit.fh"
#include "tce.fh"
#include "tce_main.fh"
      INTEGER l_a_offset
      INTEGER k_a_offset
      INTEGER size
      INTEGER length
      INTEGER addr
      INTEGER g3b,azone1,azone2,azone3,azone2_para
c ([az3],[g3]|[az1]<=[az2])
      length = 0
      azone2 = azone2_para
c      DO azone2 = 1,atpart      
      DO azone1 = 1,atpart
      DO g3b = 1,noa+nva
      DO azone3 = 1,atpart
      length = length + 1
      END DO
      END DO
c      END DO
      END DO
      IF (.not.MA_PUSH_GET(mt_int,2*length+1,'noname',l_a_offset,k_a_off
     &set)) CALL ERRQUIT('tce_t2_offset',0,MA_ERR)
      int_mb(k_a_offset) = length
      addr = 0
      size = 0
c      DO azone2 = 1,atpart      
      DO azone1 = 1,atpart 
      DO g3b = 1,noa+nva
      DO azone3 = 1,atpart
      addr = addr + 1
      int_mb(k_a_offset+addr) = azone3 - 1 + atpart * (g3b - 1 +
     &  (noa+nva) * (azone1 - 1 + atpart * (azone2 - 1)))
      int_mb(k_a_offset+length+addr) = size
      size = size + nalength(azone2) * nalength(azone1) *
     &  int_mb(k_range_alpha+g3b-1) * nalength(azone3)
      END DO
c      END DO
      END DO
      END DO
      RETURN
      END
c
c
c
      SUBROUTINE tce_offset_ggaa_osu(l_a_offset,k_a_offset,size,
     & azone2_para)
      IMPLICIT NONE
#include "global.fh"
#include "mafdecls.fh"
#include "sym.fh"
#include "errquit.fh"
#include "tce.fh"
#include "tce_main.fh"
      INTEGER l_a_offset
      INTEGER k_a_offset
      INTEGER size
      INTEGER length
      INTEGER addr
      INTEGER g3b,g4b,azone1,azone2,azone2_para
      azone2 = azone2_para
c ([g4]=>[g3]|[az1]<=[az2])
      length = 0
c      DO azone2 = 1,atpart
      DO azone1 = 1,atpart
      DO g3b = 1,noa+nva
      DO g4b = g3b,noa+nva
      length = length + 1
      END DO
      END DO
c      END DO
      END DO
      IF (.not.MA_PUSH_GET(mt_int,2*length+1,'noname',l_a_offset,k_a_off
     &set)) CALL ERRQUIT('tce_t2_offset',0,MA_ERR)
      int_mb(k_a_offset) = length
      addr = 0
      size = 0
c      DO azone2 = 1,atpart
      DO azone1 = 1,atpart
      DO g3b = 1,noa+nva
      DO g4b = g3b,noa+nva
      addr = addr + 1
      int_mb(k_a_offset+addr) = g4b - 1 + (noa+nva) * (g3b - 1 +
     &  (noa+nva) * (azone1 - 1 + atpart * (azone2 - 1)))
      int_mb(k_a_offset+length+addr) = size
      size = size + nalength(azone2) * nalength(azone1) *
     &  int_mb(k_range_alpha+g3b-1) * int_mb(k_range_alpha+g4b-1)
      END DO
c     END DO
      END DO
      END DO
      RETURN
      END
c
c
      SUBROUTINE tce_4a_offset_ns_osu(l_a_offset,k_a_offset,size,
     & azone2_pa)
C     $Id: tce_mo2e_zones_4a_disk_2s_new.F 19706 2010-10-29 17:52:31Z d3y133 $
C     This is a Fortran77 program generated by Tensor Contraction Engine v.1.0
C     Copyright (c) Battelle & Pacific Northwest National Laboratory (2002)
      IMPLICIT NONE
#include "global.fh"
#include "mafdecls.fh"
#include "sym.fh"
#include "errquit.fh"
#include "tce.fh"
#include "tce_main.fh"
      INTEGER l_a_offset
      INTEGER k_a_offset
      INTEGER size
      INTEGER length
      INTEGER addr
      INTEGER mu
      INTEGER nu
      INTEGER rho
      INTEGER sigma
      INTEGER azone1,azone2,azone3,azone4,azone2_pa
      length = 0
      azone2 = azone2_pa
      DO azone1 = 1,atpart      !nu
c      DO azone2 = 1,atpart !mu
      DO azone3 = 1,atpart      !sigma
      DO azone4 = azone3,atpart !rho
      length = length + 1
      END DO
      END DO
c     END DO
      END DO
      IF (.not.MA_PUSH_GET(mt_int,2*length+1,'noname',l_a_offset,k_a_off
     &set)) CALL ERRQUIT('tce_t2_offset',0,MA_ERR)
      int_mb(k_a_offset) = length
      addr = 0
      size = 0
      DO azone1 = 1,atpart      !nu
c      DO azone2 = 1,atpart !mu
      DO azone3 = 1,atpart      !sigma
      DO azone4 = azone3,atpart !rho
      addr = addr + 1
      int_mb(k_a_offset+addr) = azone4 - 1 + atpart * (azone3 - 1 + 
     &  atpart * (azone2 - 1 + atpart * (azone1 - 1)))
      int_mb(k_a_offset+length+addr) = size
      size = size + nalength(azone1) * nalength(azone2) * 
     &  nalength(azone3) * nalength(azone4)
      END DO
      END DO
c      END DO
      END DO
      RETURN
      END



c
      SUBROUTINE tce_offset_ggga_osu(l_a_offset,k_a_offset,size,
     & azone2_para)
      IMPLICIT NONE
#include "global.fh"
#include "mafdecls.fh"
#include "sym.fh"
#include "errquit.fh"
#include "tce.fh"
#include "tce_main.fh"
      INTEGER l_a_offset
      INTEGER k_a_offset
      INTEGER size
      INTEGER length
      INTEGER addr
      INTEGER g3b,g4b,g2b,azone2, azone2_para
c ([g4]=>[g3]|[g2][az2])
      length = 0
      azone2 = azone2_para
c      DO azone2 = 1,atpart
      DO g2b = 1,noa+nva
      DO g3b = 1,noa+nva
      DO g4b = g3b,noa+nva
      length = length + 1
      END DO
      END DO
c      END DO
      END DO
      IF (.not.MA_PUSH_GET(mt_int,2*length+1,'noname',l_a_offset,k_a_off
     &set)) CALL ERRQUIT('tce_t2_offset',0,MA_ERR)
      int_mb(k_a_offset) = length
      addr = 0
      size = 0
c      DO azone2 = 1,atpart
      DO g2b = 1,noa+nva
      DO g3b = 1,noa+nva
      DO g4b = g3b,noa+nva
      addr = addr + 1
      int_mb(k_a_offset+addr) = g4b - 1 + (noa+nva) * (g3b - 1 +
     &  (noa+nva) * (g2b - 1 + (noa+nva) * (azone2 - 1)))
      int_mb(k_a_offset+length+addr) = size
      size = size + nalength(azone2) * int_mb(k_range_alpha+g2b-1) *
     &  int_mb(k_range_alpha+g3b-1) * int_mb(k_range_alpha+g4b-1)
      END DO
      END DO
c      END DO
      END DO
      RETURN
      END
c
      subroutine tce_gacreate_sloc_osu(g_a, size, size_loc, gname)
      implicit none
#include "global.fh"
#include "mafdecls.fh"
#include "errquit.fh"
#include "stdio.fh"
      integer g_a
      character*(*) gname
      integer size
      integer size_loc
c
      if (.not.ga_create(mt_dbl,size,1,gname,
     1     size_loc,1,g_a)) then
         write(LuOut,*) ' available GA memory ',
     1        ga_memory_avail(),' bytes'
         call errquit ('tce_gacreate: failed ga_create size/nproc bytes'
     S      ,   (size*ma_sizeof(mt_dbl,1,mt_byte))/ga_nnodes(),
     1        GA_ERR)
      endif
      return
      end
