program scf
c ****************** Program SCF *******************
c This program is designed to read in the LCAO-MO
c coefficient matrix (or generate one), the one-
c and two-electron AO integrals and form a CLOSED
c SHELL Fock matrix (i.e., a Fock matrix for species
c with all doubly occupied orbitals). It then solves
c the Fock equations; iterating until convergence to
c six significant figures in the energy expression.
c A modified damping algorithm is used to insure
c convergence.
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         4 square arrays of order MAXORB
c         1 symmetric array of order MAXORB
c         3 vectors of order MAXORB
c         1 array of dimension canonical(MAXORB**4)
c     total = 4*mxosq + mxo2 + 4*maxorb + mxo4
c
      dimension pcore(4*mxosq + mxo2 + 4*maxorb + mxo4)
      dimension ipntr(10)
      itotal = 4*mxosq + mxo2 + 4*maxorb + mxo4 + 10
      call modinfo('scf         ')
      call shwmem('scf         ',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
      iorb2 = iorb*(iorb+1)/2
      iorb4 = iorb2*(iorb2+1)/2
      ipntr(1) = 1
      ipntr(2) = ipntr(1) + iorb*iorb
      ipntr(3) = ipntr(2) + iorb*iorb
      ipntr(4) = ipntr(3) + iorb2
      ipntr(5) = ipntr(4) + iorb
      ipntr(6) = ipntr(5) + iorb
      ipntr(7) = ipntr(6) + iorb
      write(*,*)' Using ',ipntr(7),' words of available memory. '
      call dos(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)),
     &         pcore(ipntr(4)),pcore(ipntr(5)),pcore(ipntr(6)),
     &         iorb,iorb2)
   20 continue
      write(*,*)
      write(*,*)' You must choose a method of generating some '
      write(*,*)' initial MO-AO guess vectors.  Your choices are: '
      write(*,*)
      write(*,*)'   1.) Use MO-AO guess vectors resulting from '
      write(*,*)'       diagonalization of the one-electron '
      write(*,*)'       Hamiltonian. '
      write(*,*)'   2.) Read the MO-AO guess vectors from disk. '
      write(*,*)'   3.) Input the MO-AO coefficients by hand. '
      write(*,*)
      write(*,*)' Input choice: < 1, 2, or 3 > '
      write(*,*)
      read(*,*)ichoos
      if(ichoos.eq.1)then
        call ezero(pcore(ipntr(2)),iorb*iorb)
      elseif(ichoos.eq.2)then
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(2)),iorb,
     &         iorb,iflnum)
c
c         Close MO coefficients file.
c
          close(unit=31)
        endif
      elseif(ichoos.eq.3)then
        call rdmatu(dscrpt,ids,pcore(ipntr(2)),iorb,iorb)
      else
        write(*,*)' Invalid choice, try again. '
        goto 20
      endif
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
      ipntr(1) = 1
      ipntr(2) = ipntr(1) + iorb*iorb
      ipntr(3) = ipntr(2) + iorb*iorb
      ipntr(4) = ipntr(3) + iorb*iorb
      ipntr(5) = ipntr(4) + iorb*iorb
      ipntr(6) = ipntr(5) + iorb4
      ipntr(7) = ipntr(6) + iorb2
      ipntr(8) = ipntr(7) + iorb
      ipntr(9) = ipntr(8) + iorb
      ipntr(10) = ipntr(9) + iorb
      write(*,*)' Using ',ipntr(10),' words of available memory. '
      call dof(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)),
     &         pcore(ipntr(4)),pcore(ipntr(5)),pcore(ipntr(6)),
     &         pcore(ipntr(7)),pcore(ipntr(8)),pcore(ipntr(9)),
     &         iorb,iorb2,iorb4)
c
c     Close 1e- and 2e- integral file.
c
      close(unit=30)
      call qmexit
      stop
 1010 format(' Maximum number of orbitals is ',I2,'. ',/,
     &       ' Please try again. ')
      end
      subroutine dos(s12,svecs,ssym,svals,scr1,scr2,n,n2)
      implicit real (a-h,o-z)
      logical yesno
      character*30 dfname, fname, dscrpt
      character*1 ans
      dimension s12(n,n),svecs(n,n),ssym(n2)
      dimension svals(n),scr1(n),scr2(n)
   10 continue
c     
c     Open the overlap file.
c     
      dscrpt = 'overlap integrals'
      ids = 17
      dfname = 'overlap'
      idf = 7
      fname = '                              '
      ifn = 30
      iflnum = 32
      call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
c
c     Read in the overlap matrix.
c
      if(yesno)then
        call rdmatd(fname,ifn,dscrpt,ids,s12,n,n,
     &              iflnum)
c
c       Close overlap integral file.
c
        close(unit=32)
      else
        write(*,*)' Would you like to input the overlap '
        write(*,*)' matrix by hand? '
        read(*,'(a)')ans
        if(ans.eq.'Y'.or.ans.eq.'y')then
          call rdmatu(dscrpt,ids,s12,n,n)
        else
          call qmexit
        endif
      endif
c
c     Put the input array into a symmetric array.
c      
      ij = 0
      do 40 i = 1,n
        do 30 j = 1,i
          ij = ij + 1
          ssym(ij) = s12(i,j)
   30   continue
   40 continue
c
c     Diagonalize the array.
c
      call diaghe(ssym,n,svals,n,svecs,n,scr1)
c      nerr = 0
c      call jacobi(n,n,ssym,svals,svecs,nerr)
c      if(nerr.ne.0)then
c        write(*,*)' Error in diagonalization routine (JACOBI). '
c        call qmexit
c      endif
c
c     Construct S exp(-1/2).
c
      do 70 i = 1,n
        svals(i) = sqrt(svals(i))
        svals(i) = 1.0e0/svals(i)
   70 continue
      do 90 i = 1,n
        do 80 j = 1,n
          s12(i,j) = svals(i)*svecs(j,i)
   80   continue
   90 continue
      ij = 0
      do 120 i = 1,n
        do 110 j = 1,i
          ij = ij + 1
          ssym(ij) = 0.0e0
          do 100 k = 1,n
            ssym(ij) = ssym(ij) + svecs(i,k)*s12(k,j)
  100     continue
  110   continue
  120 continue
c
c     Move the symmetric output array into a full size 
c     square array (SCR3).
c      
      ij = 0
      do 140 i = 1,n
        do 130 j = 1,i
          ij = ij + 1
          s12(i,j) = ssym(ij)
          s12(j,i) = ssym(ij)
  130   continue
  140 continue
      return
      end
      subroutine dof(s12,c,p,f,g,s1symn,s1n,s2n,s3n,n,n2,n4)
      implicit real (a-h,o-z)
      logical yesno
      character*30 dfname, fname, dscrpt
      character*1 ans
      dimension s12(n,n)
      dimension c(n,n),p(n,n),f(n,n)
      dimension s1symn(n2),s1n(n),s2n(n),s3n(n)
      dimension g(n4)
      fesi = 0.0e0
      ifcnt = 0
c
c     Determine the number of closed shell orbitals.
c
    5 continue
      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     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
      open(unit=99,file='scratch',access='sequential',
     &     status='unknown',form = 'unformatted')
  999 continue
c
c     Get the 1e- AO 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- AO integrals.
c
      call get2e(g,n,n4)
c
      if(ifcnt.gt.1.and.ediff.gt.0)then
        rewind(unit=99)
        read(99)p
        damp = 0.6e0
        do 4 i = 1,n
          do 3 j = 1,n
            p(i,j) = damp*p(i,j)
    3     continue
    4   continue
        adamp=1.0e0-damp
      else
        call ezero(p,n*n)
        adamp = 1.0e0
      endif
c
c     Calculate the charge bond order matrix P(MR,MS).
c
      do 30 mr = 1,n
        do 20 ms = 1,n
c         p(mr,ms) = 0.0e0
          do 10 k = 1,ndocc
            p(mr,ms) = p(mr,ms) + adamp*c(mr,k)*c(ms,k)
   10     continue
   20   continue
   30 continue
c     
c     Put copy of current iteration charge bond order matrix 
c     on disk.
c
      rewind(unit=99)
      write(99)p
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
c
c     Transform this matrix by S^(-1/2).
c
      do 100 i = 1,n
        do 90 j = 1,n
          p(i,j) = 0.0e0
          do 80 k = 1,n
            p(i,j) = p(i,j) + f(i,k)*s12(k,j)
   80     continue
   90   continue
  100 continue
      do 130 i = 1,n
        do 120 j = 1,n
          c(i,j) = 0.0e0
          do 110 k = 1,n
            c(i,j) = c(i,j) + s12(k,i)*p(k,j)
  110     continue
  120   continue
  130 continue
c
c     Pack this array (symmetric) in G.
c      
      ij = 0
      do 150 i = 1,n
        do 140 j = 1,i
          ij = ij + 1
          s1symn(ij) = c(i,j)
  140   continue
  150 continue
c
c     Diagonalize this matrix.
c
      call diaghe(s1symn,n,s1n,n,p,n,s2n)
c      nerr = 0
c      call jacobi(n,n,s1symn,s1n,p,nerr)
c      if(nerr.ne.0)then
c        write(*,*)' Error in diagonalization routine (JACOBI). '
c        call qmexit
c      endif
c
c     Back-transform to obtain the new MO coefficients.
c
      do 180 i = 1,n
        do 170 j = 1,n
          c(i,j) = 0.0e0
          do 160 k = 1,n
            c(i,j) = c(i,j) + s12(i,k)*p(k,j)
  160     continue
  170   continue
  180 continue
c
c     Transform the 1e- integrals.
c
c     Get the 1e- AO integrals.
c
      nsq = n*n
      call get1e(s1symn,n,n2)
      call t1(c,s1symn,p,n,nsq,n2)
c
c     Calculate total energy.
c
      fe2 = 0.0e0 
      do 190 k = 1, ndocc
        kk = k*(k+1)/2
        fe2 = fe2 + s1n(k) + s1symn(kk)
  190 continue
      fe2 = fe2 + repnuc
      ifcnt = ifcnt + 1
c
c     Single precision ... 7 sig figs max
c
      iconve = abs(nint(fe2))
      iconve = max0(iconve,1)
      conve = float(iconve)*5.0e-7
      ediff = fe2 - fesi 
      if(ediff.gt.0.0e0)then
        write(*,1919)ifcnt,fe2
      else
        write(*,1921)ifcnt,fe2
      endif
      if(abs(ediff).lt.conve)then
        write(*,*)' Convergence is met. '
      elseif(ifcnt.lt.30)then
        fesi  = fe2
        goto 999
      endif
 1919 format(1X,' Iteration: ',I2,' Total energy is: ',F15.7,
     &          ' up --> damp next.')
 1921 format(1X,' Iteration: ',I2,' Total energy is: ',F15.7,
     &          ' down --> OK. ')
c
c     Write the converged Fock orbital energies to user.
c
      write(*,*)' Would you like to see the orbital '
      write(*,*)' energies? '
      read(*,'(a)')ans
      if(ans.eq.'Y'.or.ans.eq.'y')then
        write(*,*)' The eigenvalues of the converged '
        write(*,*)' Fock matrix are: '
        call matout(s1n,n,1,n)
      endif
c
c     Open the Fock orbital energy file.
c     
      dscrpt = 'Fock orbital energies'
      ids = 21
      dfname = 'orbital.evals'
      idf = 13
      fname = '                              '
      ifn = 30
      iflnum = 33
      call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
        call wtmatd(fname,ifn,dscrpt,ids,s1n,n,1,iflnum)
c
c       Close orbital energies file.
c
        close(unit=33)
      endif
c
c     Write the MO-AO coefficients to user.
c
      write(*,*)' Would you like to see the converged '
      write(*,*)' MO-AO coefficients? '
      read(*,'(a)')ans
      if(ans.eq.'Y'.or.ans.eq.'y')then
        write(*,*)' The converged MO-AO coefficients are: '
        call matout(c,n,n,n)
      endif
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 wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(yesno)then
        call wtmatd(fname,ifn,dscrpt,ids,c,n,n,iflnum)
c
c       Close MO coefficients file.
c
        close(unit=31)
      endif
      close(unit=99,status='delete')
      return
 1000 format(' There is(are) ',I2,' doubly occupied orbital(s). ')
      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