program trans
c ****************** Program TRANS *****************
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 to
c transform the integrals from the AO to the MO
c basis, and write these MO integrals to a file.
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
      character*1 ans
      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)
c
c     Calculate and show memory use.
c
      itotal = 2*mxosq + mxo2 + 2*mxosq*mxo2 + 10
      call modinfo('trans       ')
      call shwmem('trans       ',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
      iorbsq = iorb*iorb
      iorb2  = iorb*((iorb+1))/2
      iorb4  = iorb2*((iorb2+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 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
c     
c     Open the MO integral file.
c     
      dscrpt = '1e- and 2e- MO integrals'
      ids = 24
      dfname = 'mo_integrals.dat'
      idf = 16
      fname = '                              '
      ifn = 30
      iflnum = 32
      call wfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      call driver(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)),
     &            pcore(ipntr(4)),pcore(ipntr(5)),iorb,iorbsq,
     &            iorb2,iorb4)
      if(yesno)then
        call ptintd(pcore(ipntr(2)),pcore(ipntr(5)),iorb,iorb2,iorb4)
        write(*,*)' The MO integrals have been written to the ',
     &            ' file: ',FNAME
c
c       Close 1e- and 2- MO integral file.
c
        close(unit=32)
      endif
      write(*,*)' Do you want to see the 1e- or 2e- MO '
      write(*,*)' integrals?  <y or n> '
      read(*,'(a)')ans
      if(ans.eq.'Y'.or.ans.eq.'y')then
        call ptintu(pcore(ipntr(2)),pcore(ipntr(5)),iorb,iorb2,iorb4)
      endif
      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,n4)
      implicit real (a-h,o-z)
      dimension c(n,n), h(n2), p(n,n), q(nsq*n2), g(nsq*n2)
      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. '
      return
      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
      subroutine ptintd(h,g,n,n2,n4)
      implicit real (a-h,o-z)
      dimension h(n2), g(n4)
c     
c     Write a canonical list of 1 and 2 electron integrals
c     to unit 32.
c
      rewind(unit=32)
c     
c     Write the one-electron integrals.
c     
      il = n2
      do 10 i = 1,il
        write(32,1000)h(i)
   10 continue
c
c     Write the two-electron integrals.
c
      ind = 0
      do 50 i = 1,n
        do 40 k = 1,i
          do 30 j = 1,i
            lend = j
            if(i.eq.j)lend = k
            do 20 l = 1,lend
              ind = ind + 1
              write(32,1000)g(ind)
   20       continue
   30     continue
   40   continue
   50 continue
      return
 1000 format(1x,f15.8)
      end
      subroutine ptintu(h,g,n,n2,n4)
      implicit real (a-h,o-z)
      character*1 ans
      dimension h(n2), g(n4)
c     
c     Write a canonical list of 1 and 2 electron integrals
c     to the user.
c     
c     Write the one-electron integrals.
c     
      write(*,*)' Do you want to see the 1e- MO '
      write(*,*)' integrals?  <y or n> '
      read(*,'(a)')ans
      if(ans.eq.'Y'.or.ans.eq.'y')then
        write(*,*)' The one-electron MO integrals: '
        call tprntd(h,n)
      endif
c
C     Write the two-electron integrals.
C
      write(*,*)' Do you want to see the 2e- MO '
      write(*,*)' integrals?  <y or n> '
      read(*,'(a)')ans
      if(ans.eq.'Y'.or.ans.eq.'y')then
        write(*,*)' The two-electron MO integrals: '
        write(*,*)
        ind = 0
        do 50 i = 1,n
          do 40 k = 1,i
            do 30 j = 1,i
              lend = j
              if(i.eq.j)lend = k
              do 20 l = 1,lend
                ind = ind + 1
                write(*,1010)i,j,k,l,g(ind)
   20         continue
   30       continue
   40     continue
   50   continue
      endif
      return
 1000 format(1x,f15.8)
 1010 format(' i = ',i2,' j = ',i2,' k = ',i2,' l = ',i2,
     &       ' g(i,j,k,l) = ',f15.8)
      end