program fock
c ****************** Program FOCK ******************
c This program is designed to read in the LCAO-MO
c coefficient matrix, the one- and two-electron AO
c integrals and to form a CLOSED SHELL Fock matrix
c (i.e., a Fock matrix for species with all doubly
c occupied orbitals).
c
c Copyright (C) 1992 Jeff Nichols and Jack Simons
c
c **************************************************
      implicit real (a-h,o-z)
      include "limits.h"
      logical yesno
      character*30 dfname, fname, dscrpt
c
c     Set core needed to:
c         3 square arrays of order maxorb
c         1 array of dimension canonical(maxorb**4)
c     total = 3*mxosq + mxo4
c
      dimension pcore(3*mxosq + mxo4)
      dimension ipntr(10)
      itotal = 3*mxosq + mxo4 + 10
      call modinfo('fock        ')
      call shwmem('fock        ',itotal)
      call goon
   10 continue
      write(*,*)
      write(*,*)' Input the number of orbitals in your system: '
      read(*,*)iorb
      if(iorb.gt.maxorb)then
        write(*,1010)maxorb
        goto 10
      endif
      ipntr(1) = 1
      ipntr(2) = ipntr(1) + iorb*iorb
      ipntr(3) = ipntr(2) + iorb*iorb
      ipntr(4) = ipntr(3) + iorb*iorb
      iorb2 = iorb*(iorb+1)/2
      iorb4 = iorb2*(iorb2+1)/2
      ipntr(5) = ipntr(4) + iorb4
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 integral file.
c     
      dscrpt = '1e- and 2e- integrals'
      ids = 21
      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)),iorb,iorb2,iorb4)
      call qmexit
      stop
 1010 format(' Maximum number of orbitals is ',I2,'. ',/,
     &       ' Please try again. ')
      end
      subroutine driver(c,p,f,g,n,n2,n4)
      implicit real (a-h,o-z)
      logical yesno
      character*30 dfname, fname, dscrpt
      dimension c(n,n),p(n,n),f(n,n),g(n4)
c
c     Get the 1e- integrals.
c
      call get1e(g,n,n2)
c
c     Put the 1e- Hamiltonian directly into the Fock matrix.
c
      ij = 0
      do 2 i = 1,n
        do 1 j = 1,i
          ij = ij + 1
          f(i,j) = g(ij)
          f(j,i) = f(i,j)
    1   continue
    2 continue
c
c     Get the 2e- integrals.
c
      call get2e(g,n,n4)
c
c     Close 1e- and 2e- integrals file.
c
      close(unit=30)
    5 continue
c
c     Determine the number of closed shell orbitals.
c
      write(*,*)' Please input the number of doubly occupied '
      write(*,*)' 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     Calculate the charge bond order matrix P(MR,MS).
c
        do 30 mr = 1,n
          do 20 ms = 1,n
            p(mr,ms) = 0.0e0
            do 10 k = 1,ndocc
              p(mr,ms) = p(mr,ms) + c(mr,k)*c(ms,k)
   10       continue
   20     continue
   30   continue
c
c     Write out the charge bond order matrix.
c
      write(*,*)' The charge bond order matrix: '
      call matout(p,n,n,n)
c
c     Form the CLOSED SHELL Fock matrix.
c     
      do 70 mu = 1,n
        do 60 nu = 1,n
          do 50 mr = 1,n
            do 40 ms = 1,n
              icx1 = icanon(mu,mr,nu,ms)
              icx2 = icanon(mu,mr,ms,nu)
              f(mu,nu) = f(mu,nu) + 
     &                   p(mr,ms)*(2.0e0*g(icx1) - g(icx2))
   40       continue
   50     continue
   60   continue
   70 continue
      dscrpt = 'Fock matrix'
      ids = 11
      call  wtmatu(dscrpt,ids,f,n,n)
      dfname = 'fock.out'
      idf = 8
      fname = '                              '
      ifn = 30
      iflnum = 32
      call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
        call wtmatd(fname,ifn,dscrpt,ids,f,n,n,
     &              iflnum)
c
c       Close Fock matrix file.
c
        close(unit=32)
      endif
      return
 1000 format(' There is ',I2,' doubly occupied orbital(s). ')
      end