program mp2
c ****************** Program MP2 *******************
c This program is designed to read in the
c transformed one- and two-electron integrals and  
c the Fock orbital energies after which it will
c compute the second order Moller Plesset
c perturbation theory energy (MP2).
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         1 vector of dimension MAXORB
c         2 symmetric arrays of order MAXORB
c         1 array of dimension canonical(MAXORB**4)
c     total = 2*(maxorb*(maxorb+1))/2 + 
c         (maxorb*(maxorb+1))/2*(((maxorb*(maxorb+1))/2+1))/2 + 
c         maxorb    
c
      dimension pcore(2*mxo2 + mxo4 + maxorb)
      dimension ipntr(10)
      itotal = 2*mxo2 + mxo4 + maxorb + 10
      call modinfo('mp2         ')
      call shwmem('mp2         ',itotal)
      call goon
    1 continue
      write(*,*)
      write(*,*)' Input the number of orbitals in your system: '
      read(*,*)iorb
      if(iorb.gt.maxorb)then
        write(*,1010)maxorb
        goto 1
      endif
    5 continue
      write(*,*)' Please input the number of doubly occupied '
      write(*,*)' orbitals. '
      read(*,*)ndocc
      write(*,1000)ndocc
      if(ndocc.gt.iorb)then
        write(*,*)' The number of doubly occupied orbitals cannot '
        write(*,*)' exceed the total number of orbitals.  Please '
        write(*,*)' try again. '
        goto 5
      endif
      iorb2 = iorb*(iorb+1)/2
      iorb4 = iorb2*(iorb2+1)/2
      ipntr(1) = 1
      ipntr(2) = ipntr(1) + iorb2
      ipntr(3) = ipntr(2) + iorb4
      ipntr(4) = ipntr(3) + iorb
      ipntr(5) = ipntr(4) + iorb2
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 = 30
      call rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(.not.yesno)then
        write(*,*)' You must have the MO integrals available '
        write(*,*)' on disk to construct the MP2 energy. '
        call qmexit
      endif
c
c     Get the 1e- AO integrals.
c
      call get1e(pcore(ipntr(1)),iorb,iorb2)
c
c     Get the 2e- AO integrals.
c
      call get2e(pcore(ipntr(2)),iorb,iorb4)
c
c     Close 1e- and 2- MO integral file.
c
      close(unit=30)
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 rfile(fname,ifn,dscrpt,ids,dfname,idf,iflnum,yesno)
      if(.not.yesno)then
        write(*,*)' You must have the Fock orbital energies available '
        write(*,*)' on disk to construct the MP2 energy. '
        call qmexit
      endif
c
c     Get the orbital energies.
c
      rewind(unit=33)
      do 10 i = 1,iorb
        read(33,*)pcore(ipntr(3)+i-1)
   10 continue
c
c     Close orbital energies file.
c
      close(unit=33)
c
c     Compute MP2 energy.
c
      call domp2(pcore(ipntr(1)),pcore(ipntr(2)),pcore(ipntr(3)),
     &           pcore(ipntr(4)),ndocc,iorb,iorb2,iorb4)
      call qmexit
      stop
 1000 format(' There is(are) ',i2,' doubly occupied orbital(s). ')
 1010 format(' Maximum number of orbitals is ',I2,'. ',/,
     &       ' Please try again. ')
      end
      subroutine domp2(h,g,e,ep,nc,n,n2,n4)
      implicit real (a-h,o-z)
      include "limits.h"
      dimension h(n2), g(n4), e(n), ep(n2)
      isym2(i,j)=(max0(i,j)*(max0(i,j)-1))/2 + min0(i,j)
      isym4(i,j,k,l)=(max0(isym2(i,j),isym2(k,l))*
     &               (max0(isym2(i,j),isym2(k,l))-1))/2 + 
     &               min0(isym2(i,j),isym2(k,l))
      call ezero(ep,n2)
c
c     need to do all spin combos e.g., aaaa, bbbb, abab, etc.
c     ib, ia, is, and ir label spin orbitals
c     ibspac, iaspac, isspac, and irspac label spacial orbitals
c     ibs, ias, iss, and irs give alpha or beta labels, 0=alpha, 1=beta
c
      do 40 ib = 1, 2*nc
        ibspac = (ib+1)/2
        if(ib-2*(ib/2).eq.0)then
          ibs=1
        else
          ibs=0
        endif
        do 30 ia = 1, ib-1
          iaspac = (ia+1)/2
          if(ia-2*(ia/2).eq.0)then
            ias=1
          else
            ias=0
          endif
          ibsias = isym2(ibspac,iaspac)
          do 20 is = 2*nc+1, 2*n
            isspac = (is+1)/2
            if(is-2*(is/2).eq.0)then
              iss=1
            else
              iss=0
            endif
            do 10 ir = 2*nc+1, is-1
              irspac = (ir+1)/2
              if(ir-2*(ir/2).eq.0)then
                irs=1
              else
                irs=0
              endif
              iabrs = isym4(iaspac,irspac,ibspac,isspac)
              iabsr = isym4(iaspac,isspac,ibspac,irspac)
              if(ias.eq.irs)then
                arspin = 1.0e0
              else
                arspin = 0.0e0
              endif
              if(ibs.eq.iss)then
                bsspin = 1.0e0
              else
                bsspin = 0.0e0
              endif
              if(ias.eq.iss)then
                asspin = 1.0e0
              else
                asspin = 0.0e0
              endif
              if(ibs.eq.irs)then
                brspin = 1.0e0
              else
                brspin = 0.0e0
              endif
              giabrs = g(iabrs)*arspin*bsspin
              giabsr = g(iabsr)*asspin*brspin
              et = ((giabrs-giabsr)**2)/
     &             (e(irspac)+e(isspac)-e(iaspac)-e(ibspac))
c              write(*,*)' ia,ib,ir,is,iabrs,iabsr,giabrs,giabsr,et ',
c     &                    ia,ib,ir,is,iabrs,iabsr,giabrs,giabsr,et
              ep(ibsias) = ep(ibsias) - et
   10       continue
   20     continue
   30   continue
   40 continue
      write(*,*)' Pair correlation energies: '
      call tprntd(ep,nc)
      emp2 = 0.0e0
      do 50 i = 1, n2
        emp2 = emp2 + ep(i)
   50 continue
      write(*,*)' MP2 energy = ',emp2
      return
      end