program fnrgy
c **************** Program FENERGY *****************
c This program is designed to read in the LCAO-MO
c coefficient matrix, the one- and two-electron AO
c integrals (in Dirac <12|12> convention), and the
c Fock orbital energies.  Upon transformation of the
c one- and two-electron integrals from the AO to the
c MO basis, the closed shell Hartree - Fock energy
c is calculated in two ways.  First, the energy is 
c calculated with the MO integrals,
c 
c Sum(k) 2*<k|h|k> + Sum(k,l) (2*<k,l|k,l> -
c <k,l|l,k>) + ZuZv/Ruv.
c 
c Secondly, the energy is calculated with the Fock
c orbital energies and one electron energies in the
c MO basis,
c 
c Sum(k) (eps(k) + <k|h|k>) + ZuZv/Ruv.
c
c Copyright (C) 1992 Jeff Nichols and Jack Simons
c
c **************************************************
      implicit real (a-h,o-z)
      logical yesno
      character*30 dfname, fname, dscrpt
      include "limits.h"
c
c     Set core needed to:
c         2 square arrays of order maxorb
c         1 symmetric array of order maxorb
c         2 arrays of dimension:
c           maxorb*maxorb*((maxorb*(maxorb+1))/2)
c     total = 2*maxorb**2 + (maxorb*(maxorb+1))/2
c             2*maxorb*maxorb*((maxorb*(maxorb+1))/2)
c
      dimension pcore(2*mxosq + mxo2 + 2*mxosq*mxo2)
      dimension ipntr(10)
      itotal = 2*mxosq + mxo2 + 2*mxosq*mxo2 + 10
      call modinfo('fenergy     ')
      call shwmem('fenergy     ',itotal)
   10 continue
      call goon
      write(*,*)
      write(*,*)' Input the number of orbitals in your system: '
      read(*,*)iorb
      if(iorb.gt.maxorb)then
        write(*,1010)maxorb
        goto 10
      endif
      iorbsq = iorb*iorb
      iorb2  = (iorb*(iorb+1))/2
      ipntr(1) = 1
      ipntr(2) = ipntr(1) + iorbsq
      ipntr(3) = ipntr(2) + iorb2
      ipntr(4) = ipntr(3) + iorbsq
      ipntr(5) = ipntr(4) + iorbsq*iorb2
      ipntr(6) = ipntr(5) + iorbsq*iorb2
c     
c     Open the MO-AO coefficient file.
c     
      dscrpt = 'MO-AO coefficients'
      ids = 18
      dfname = 'mocoefs.dat'
      idf = 11
      fname = '                              '
      ifn = 30
      iflnum = 31
      call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
        call rdmatd(fname,ifn,dscrpt,ids,pcore(ipntr(1)),iorb,
     &              iorb,iflnum)
c
c       Close MO-AO coefficients file.
c
        close(unit=31)
      else
        call rdmatu(dscrpt,ids,pcore(ipntr(1)),iorb,iorb)
      endif
      call wtmatu(dscrpt,ids,pcore(ipntr(1)),iorb,iorb)
c     
c     Open the AO integral file.
c     
      dscrpt = '1e- and 2e- AO integrals'
      ids = 24
      dfname = 'ao_integrals.dat'
      idf = 16
      fname = '                              '
      ifn = 30
      iflnum = 30
      call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(.not.yesno)then
        write(*,*)' You must have the AO integrals available '
        write(*,*)' on disk ... this can be done by running '
        write(*,*)' program INTEGRAL. '
        call qmexit
      endif
      call driver(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)),
     &            pcore(ipntr(4)),pcore(ipntr(5)),iorb,iorbsq,iorb2)
      call qmexit
      stop
 1010 format(' Maximum number of orbitals is ',I2,'. ',/,
     &       ' Please try again. ')
      end
      subroutine driver(c,h,p,q,g,n,nsq,n2)
      implicit real (a-h,o-z)
      logical yesno
      character*30 dfname, fname, dscrpt
      dimension c(n,n), h(n2), p(n,n), q(nsq*n2), g(nsq*n2)
      n4 = n2*(n2+1)/2
      call get1e(h,n,n2)
      call get2e(g,n,n4)
c
c     Close 1e- and 2- AO integral file.
c
      close(unit=30)
c
c     Transform the one-electron integrals.
c
      write(*,*)' Begin 1e- integral transformation. '
      call t1(c,h,p,n,nsq,n2)
      write(*,*)' End 1e- integral transformation. '
c
c     Transform the two-electron integrals.
c
      write(*,*)' Begin 2e- integral transformation. '
c
c     Do the first quarter transformation.
c
      write(*,*)' Begin 1st quarter transformation. '
      call q1(c,q,g,n,nsq,n2)
      write(*,*)' End 1st quarter transformation. '
c
c     Do the second quarter transformation.
c     
      write(*,*)' Begin 2nd quarter transformation. '
      call q2(c,q,g,n,nsq,n2)
      write(*,*)' End 2nd quarter transformation. '
c
c     Do the third quarter transformation.
c     
      write(*,*)' Begin 3rd quarter transformation. '
      call q3(c,q,g,n,nsq,n2)
      write(*,*)' End 3rd quarter transformation. '
c
c     Do the fourth quarter transformation.
c     
      write(*,*)' Begin 4th quarter transformation. '
      call q4(c,q,g,n,nsq,n2)
      write(*,*)' End 4th quarter transformation. '
      write(*,*)' End 2e- integral transformation. '
c
c     Open the Fock orbital energy file.
c     
      dscrpt = 'Fock orbital energies'
      ids = 21
      dfname = 'orbital.evals'
      idf = 13
      fname = '                              '
      ifn = 30
      iflnum = 32
      call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
        call rdmatd(fname,ifn,dscrpt,ids,p,n,1,iflnum)
c
c       Close energies file.
c
        close(unit=32)
      else
        call rdmatu(dscrpt,ids,p,n,1)
      endif
      call wtmatu(dscrpt,ids,p,n,1)
    5 continue
c
c     Determine the number of closed shell orbitals.
c
      write(*,*)' Input the number of doubly occupied orbitals. '
      read(*,*)ndocc
      write(*,1000)ndocc
      if(ndocc.gt.n)then
        write(*,*)' The number of doubly occupied orbitals cannot '
        write(*,*)' exceed the total number of orbitals.  Please '
        write(*,*)' try again. '
        goto 5
      endif
c
c     Get the nuclear repusion.
c
      dscrpt = 'nuclear repulsion'
      ids = 17
      dfname = 'NucRep'
      idf = 6
      fname = '                              '
      ifn = 30
      iflnum = 31
      call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
c
c     Read in the nuclear repulsion.
c
      if(yesno)then
        call rdmatd(fname,ifn,dscrpt,ids,repnuc,1,1,iflnum)
c
c       Close nuclear repulsion file.
c
        close(unit=31)
      else
        write(*,*)' Input the nuclear repulsion energy. '
        read(*,*)repnuc
      endif
c
c     Calculate energy using MO integrals and 
c     using Fock energies.
c
      fe1 = 0.0e0 
      fe2 = 0.0e0 
      do 101 k = 1, ndocc
        kk = (k*(k+1))/2
        fe1 = fe1 + h(kk)
        fe2 = fe2 + p(k,1) + h(kk)
  101 continue
      fe1 = 2.0e0*fe1
      do 103 k = 1, ndocc
        do 102 l = 1, ndocc
          iklkl = icanon(k,l,k,l)
          ikllk = icanon(k,l,l,k)
          fe1 = fe1 + 2.0e0*g(iklkl) - g(ikllk)
  102   continue
  103 continue
      fe1 = fe1 + repnuc
      fe2 = fe2 + repnuc
      ediff = fe1 - fe2
      write(*,1010)fe1,fe2,ediff
      return
 1000 format(' There is ',I2,' doubly occupied orbital(s). ')
 1010 format(' The closed shell Fock energy from: ',/,
     &       ' formula: ',//,
     &       ' Sum(k) 2*<k|h|k> + Sum(k,l) (2*<k,l|k,l> - ',/,
     &       ' <k,l|l,k>) + ZuZv/Ruv ',//,
     &       ' is: ',F15.8,//,
     &       ' from formula: ',//,
     &       ' Sum(k) (eps(k) + <k|h|k>) + ZuZv/Ruv ',//,
     &       ' is: ',F15.8,//,
     &       ' the difference is: ',F15.8)
      end
      subroutine q1(c,q,g,n,nsq,n2)
      implicit real (a-h,o-z)
      dimension c(n,n), q(nsq*n2), g(nsq*n2)
c
c     Do the first quarter transformation.
c
      nl = nsq*n2
      call ezero(q,nl)
      nsq = n*n
      iabsym = 0
      icx = 0
      do 50 ia = 1,n
        iam1 = ia - 1
        iasq = iam1*n
        do 40 ib = 1,ia
          ibm1 = ib - 1
          ibsq = ibm1*n
          iabsym = iabsym + 1
          iabsnn = (iabsym - 1)*nsq
          do 30 ic = 1,ia
            icm1 = ic - 1
            icsq = icm1*n
            icsym = ic*icm1/2
            idend = ic
            if(ic.eq.ia)idend = ib
            do 20 id = 1,idend
              idm1 = id - 1
              idsq = idm1*n
              icdsym = (icsym + id)
              icdsnn = (icdsym - 1)*nsq
              icx = icx + 1
              xg = g(icx)
              do 10 l = 1,n
                issx = iabsnn + icsq + l
                q(issx) = q(issx) + xg*c(id,l)
                if(ic.ne.id)then
                  issx = iabsnn + idsq + l 
                  q(issx) = q(issx) + xg*c(ic,l)
                endif
                if(iabsym.ne.icdsym)then
                  issx = icdsnn + iasq + l
                  q(issx) = q(issx) + xg*c(ib,l)
                  if(ia.ne.ib)then
                    issx = icdsnn + ibsq + l
                    q(issx) = q(issx) + xg*c(ia,l)
                  endif
                endif
   10         continue
   20       continue
   30     continue
   40   continue
   50 continue
      return
      end
      subroutine q2(c,q,g,n,nsq,n2)
      implicit real (a-h,o-z)
      dimension c(n,n), q(nsq*n2), g(nsq*n2)
      nl = nsq*n2
      call ezero(g,nl)
c
c     Do the second quarter transformation.
c     
      nsq = n*n
      iabsym = 0
      do 150 ia = 1,n
        do 140 ib = 1,ia
          iabsym = iabsym + 1
          iabsnn = (iabsym - 1)*nsq
          do 130 ic = 1,n
            icm1 = ic - 1
            icsq = icm1*n
            do 120 l = 1,n
              issx = iabsnn + icsq + l
              xq = q(issx)
              do 110 k = l,n
                 km1 = k - 1
                 ksq = km1*n
                issx = iabsnn + ksq + l
                g(issx) = g(issx) + xq*c(ic,k)
  110         continue
  120       continue
  130     continue
  140   continue
  150 continue
      return
      end
      subroutine q3(c,q,g,n,nsq,n2)
      implicit real (a-h,o-z)
      dimension c(n,n), q(nsq*n2), g(nsq*n2)
      nl = nsq*n2
      call ezero(q,nl)
c
c     Do the third quarter transformation.
c     
      nsq = n*n
      klsym = 0
      do 250 k = 1,n
        km1 = k - 1
        ksq = km1*n
        do 240 l = 1,k
          klsym = klsym + 1
          klsnn = (klsym - 1)*nsq
          iabsym = 0
          do 230 ia = 1,n
            iam1 = ia - 1
            iasq = iam1*n
            do 220 ib = 1,ia
              ibm1 = ib - 1
              ibsq = ibm1*n
              iabsym = iabsym + 1
              iabsnn = (iabsym - 1)*nsq
              issx = iabsnn + ksq + l
              xg = g(issx)
              do 210 j = 1,n
                issx = klsnn + iasq + j
                q(issx) = q(issx) + xg*c(ib,j)
                if(ia.ne.ib)then
                  issx = klsnn + ibsq + j
                  q(issx) = q(issx) + xg*c(ia,j)
                endif
  210         continue
  220       continue
  230     continue
  240   continue
  250 continue
      return
      end
      subroutine q4(c,q,g,n,nsq,n2)
      implicit real (a-h,o-z)
      dimension c(n,n), q(nsq*n2), g(nsq*n2)
      nl = nsq*n2
      call ezero(g,nl)
c
c     Do the fourth quarter transformation.
c     
      nsq = n*n
      klsym = 0
      do 350 k = 1,n
        do 340 l = 1,k
          klsym = klsym + 1
          klsnn = (klsym - 1)*nsq
          do 330 ia = 1,n
            iam1 = ia - 1
            iasq = iam1*n
            do 320 j = 1,n
              issx = klsnn + iasq + j
              xq = q(issx)
              istart = max(j,k)
              if(l.gt.j)istart = k + 1
              do 310 i = istart,n
                icx = icanon(i,k,j,l)
                g(icx) = g(icx) + xq*c(ia,i)
  310         continue
  320       continue
  330     continue
  340   continue
  350 continue
      return
      end
      subroutine t1(c,h,p,n,nsq,n2)
      implicit real (a-h,o-z)
      dimension c(n,n),h(n2),p(n,n)
c
c     Transform the one-electron integrals.
c
      nl = nsq
      call ezero(p,nl)
      iaib = 0
      do 430 ia = 1,n
        do 420 ib = 1,ia
          iaib = iaib + 1
          do 410 j = 1,n
            p(ia,j) = p(ia,j) + h(iaib)*c(ib,j)
            if(ia.ne.ib)then
              p(ib,j) = p(ib,j) + h(iaib)*c(ia,j)
            endif
  410     continue
  420   continue
  430 continue
      nl = n2
      call ezero(h,nl)
      do 530 ia = 1,n
        ij = 0
        do 520 j = 1,n
          DO 510 I = 1,J
            ij = ij + 1
            h(ij) = h(ij) + p(ia,j)*c(ia,i)
  510     continue
  520   continue
  530 continue
      return
      end