      function MECI(EIGS,COEFF)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension EIGS(NORBS), COEFF(NORBS,NORBS)
***********************************************************************
*
*                 PROGRAM MECI
*
*   A MULTI-ELECTRON CONFIGURATION INTERACTION CALCULATION
*
*   WRITTEN BY JAMES J. P. STEWART, AT THE
*              FRANK J. SEILER RESEARCH LABORATORY
*              USAFA, COLORAdo SPRINGS, CO 80840
*
*              1985
*
***********************************************************************
C
      double precision MECI
C
C   MATRICES FOR PERMUTATION WORK
C
      dimension NFA(2*NMECI), NPERMA(NMECI,6*NMECI),
     1NPERMB(NMECI,6*NMECI), EIGA(NMECI)
C
C   MATRICES FOR ONE AND TWO ELECTRON INTEGRALS
C
      common /RJKS  /RJKAB(NMECI,NMECI), RJKAA(NMECI,NMECI)
C
C   SPIN MATRICES
C
      dimension SPIN(NMECI**2)
      logical DEBUG,  LARGE, PRNT, LSPIN, LSPIN1,
     1 FIRST1, BIGPRT, SING, DOUB, TRIP, QUAR, QUIN, SEXT,
     2 PRNT2, FORCE, GEOOK
      character KEYWRD*241, TSPIN(7)*8, LINE*80
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,
     2                NDUMMY(2), NCLOSE, NOPEN, NDUMY, FRACT
      common /LAST  / LAST
      common /SPQR/ ISPQR(NMECI**2,NMECI),IS,I,K
      common /KEYWRD/ KEYWRD
C
C   MATRICES FOR SEC.DET., VECTORS, AND EIGENVALUES.
C
      common /WORK2 / CIMAT(NMECI**4), EIG(NMECI**2), DIAG(2*NMECI**3)
      common /BASEOC/ OCCA(NMECI)
     1       /WORK3 / DIJKL(MPACK*4)
      common /CIVECT/ VECTCI(NMECI**2),CONF(NMECI**4)
      common /NALMAT/ NALPHA(NMECI**2)
      common /MICROS/ MICROA(NMECI,4*NMECI**2), MICROB(NMECI,4*NMECI**2)
      common /CIBITS/ NMOS,LAB,NELEC, NBO(3)
      common /XYIJKL/ XY(NMECI,NMECI,NMECI,NMECI)
      common /NUMCAL/ NUMCAL
      save FIRST1, TSPIN
      save J,L, DEBUG, PRNT2,MDIM, LSPIN1, FORCE
      save LARGE, LROOT, SING, DOUB, QUAR, QUIN, SEXT, SMULT, NE
      save NLEFT, GEOOK
      data ICALCN/0/
      data TSPIN/'SINGLET ','DOUBLET ','TRIPLET ','QUARTET ','QUINTET ',
     1'SEXTET  ','SEPTET  '/
      if (ICALCN.NE.NUMCAL) then
         ICALCN=NUMCAL
         FIRST1=.TRUE.
         MDIM=NMECI**2
         GEOOK=(INDEX(KEYWRD,'GEO-OK').NE.0)
         FORCE=(INDEX(KEYWRD,'FORCE').NE.0)
         LSPIN1=(INDEX(KEYWRD,'ESR').NE.0)
         DEBUG=(INDEX(KEYWRD,'DEBUG').NE.0)
         PRNT2=(INDEX(KEYWRD,'MECI').NE.0)
         DEBUG=(DEBUG.AND.PRNT2)
         LARGE=(INDEX(KEYWRD,'LARGE').NE.0)
         NDOUBL=99
         if (INDEX(KEYWRD,'C.I.=(').NE.0) then
            NDOUBL = READA(KEYWRD,INDEX(KEYWRD,'C.I.=(')+7)
            NMOS = READA(KEYWRD,INDEX(KEYWRD,'C.I.=(')+5)
         else if (INDEX(KEYWRD,'C.I.=').NE.0) then
            NMOS = READA(KEYWRD,INDEX(KEYWRD,'C.I.=')+5)
         else
            NMOS=NOPEN-NCLOSE
         end if
         LROOT=1
         if (INDEX(KEYWRD,'EXCI').NE.0) LROOT=2
         I=INDEX(KEYWRD,'ROOT')
         if (I.NE.0) LROOT = READA(KEYWRD,I)
         if (NDOUBL.EQ.99) then
            J=MAX(MIN((NCLOSE+NOPEN+1)/2-(NMOS-1)/2,NORBS-NMOS+1),1)
         else
            J=NCLOSE-NDOUBL+1
            if (FRACT.GT.1.99D0) J=J+1
         end if
         L=0
         do 10 I=J,NCLOSE
            L=L+1
   10    OCCA(L)=1
         do 20 I=NCLOSE+1,NOPEN
            L=L+1
   20    OCCA(L)=FRACT*0.5D0
         do 30 I=NOPEN+1,J+NMOS-1
            L=L+1
   30    OCCA(L)=0.D0
C#         write(MIOB,'('' INITIAL ORBITAL OCCUPANCIES'')')
C#         write(MIOB,'(6F12.6)') (OCCA(L),L=1,NMOS)
         SING=(INDEX(KEYWRD,'SING')+
     1         INDEX(KEYWRD,'EXCI')+
     2         INDEX(KEYWRD,'BIRAD').NE.0)
         DOUB=(INDEX(KEYWRD,'DOUB').NE.0)
         TRIP=(INDEX(KEYWRD,'TRIP').NE.0)
         QUAR=(INDEX(KEYWRD,'QUAR').NE.0)
         QUIN=(INDEX(KEYWRD,'QUIN').NE.0)
         SEXT=(INDEX(KEYWRD,'SEXT').NE.0)
C
C  DEFINE MAGNETIC COMPONENT OF SPIN
C
         MSDEL=INDEX(KEYWRD,' MS')
         if (MSDEL.NE.0) then
            MSDEL=1.0001D0*READA(KEYWRD,INDEX(KEYWRD,' MS'))
         else
            if (TRIP.OR.QUAR) MSDEL=1
            if (QUIN.OR.SEXT) MSDEL=2
         end if
         SMULT=-.5D0
         if (SING) SMULT=0.00D0
         if (DOUB) SMULT=0.75D0
         if (TRIP) SMULT=2.00D0
         if (QUAR) SMULT=3.75D0
         if (QUIN) SMULT=6.00D0
         if (SEXT) SMULT=8.75D0
         X=0.D0
         do 40 J=1,NMOS
   40    X=X+OCCA(J)
         XX=X+X
         NE=XX+0.5
         NELEC=(NELECS-NE+1)/2
         NLEFT=NORBS-NMOS-NELEC
      end if
      PRNT=(DEBUG.OR.LAST.EQ.3.AND.PRNT2)
      BIGPRT=(PRNT.AND.LARGE)
C
C    TEST TO SEE IF THE SET OF ENERGY LEVELS USED IN MECI IS COMPLETE,
C    I.E., ALL COMPONENTS OF DEGENERATE IRREDUCIBLE REPRESENTATIONS
C    ARE USED.  IF NOT, THEn RESULTS WILL BE NONSENSE.  GIVE USERS A
C    CHANCE TO REALLY FOUL THINGS UP BY ALLOWING JOB TO continue IF
C    'GEO-OK' IS SPECifIED.
C
      do 50 I=1,NMOS
         IN=I+NELEC
   50 EIGA(I)=EIGS(IN)
      LSPIN=(LSPIN1.AND. LAST.EQ.3)
      if (BIGPRT) then
         write(MIOB,'(''  INITIAL EIGENVALUES'')')
         write(MIOB,'(5F12.6)') (EIGA(I),I=1,NMOS)
         write(MIOB,'(//10X,''NUMBER OF ELECTRONS IN C.I. ='',F5.1)') XX
      end if
      if (.NOT.GEOOK.AND.NELEC.GT.0) then
         if (ABS(EIGS(NELEC+1)-EIGS(NELEC)).LT.1.D-1.OR.
     1ABS(EIGS(NELEC+1+NMOS)-EIGS(NELEC+NMOS)).LT.1.D-1) then
            write(MIOB,'(///10X,A)')
     %        'DEGENERATE ENERGY LEVELS DETECTED IN MECI'
            write(MIOB,'(10X,A)')
     %        'SOME OF THESE LEVELS WOULD BE TREATED BY MECI,'
            write(MIOB,'(10X,A)')
     %        'WHILE OTHERS WOULD NOT.  THIS WOULD RESULT IN'
            write(MIOB,'(10X,A)')
     %        'NON-REPRODUCIBLE ELECTRONIC ENERGIES.'
            write(MIOB,'(10X,A)')
     %        '  JOB STOPPED.  TO continue, SPECifY "GEO-OK"'
            call filclo ! J/M
            stop ' in MECI - degenerate energy levels found'
         end if
      end if
      if ( BIGPRT ) then
         write(MIOB,'(//10X,''EIGENVECTORS'',/)')
         do 60 I=1,NORBS
   60    write(MIOB,'(6F12.6)') (COEFF(I,J+NELEC),J=1,NMOS)
      end if
      NFA(2)=1
      NFA(1)=1
      do 70 I=3,NMECI+1
   70 NFA(I)=NFA(I-1)*(I-1)
      call IJKL(COEFF(1,NELEC+1),COEFF,NELEC,NMOS,DIJKL)
      do 80 I=1,NMOS
         do 80 J=1,NMOS
            RJKAA(I,J)=XY(I,I,J,J)-XY(I,J,I,J)
   80 RJKAB(I,J)=XY(I,I,J,J)
      do 100 I=1,NMOS
         X=0.0
         do 90 J=1,NMOS
            X=X+(RJKAA(I,J)+RJKAB(I,J))*OCCA(J)
   90    continue
         EIGA(I)=EIGA(I)-X
C#      if (ABS(OCCA(I)-0.5).LT.1.D-4) EIGA(I)=EIGA(I)+XY(I,I,I,I)*0.25D0
  100 continue
      if (BIGPRT) then
         write(MIOB,110)
  110    format(/,5X,'EIGENVALUES AFTER REMOVAL OF INTER-ELECTRONIC INTE
     1RACTIONS',/)
         write(MIOB,'(6F12.6)') (EIGA(I),I=1,NMOS)
         write(MIOB,'(///10X,''TWO-ELECTRON J-INTEGRALS'',/)')
         do 120 I1=1,NMOS
  120    write(MIOB,'(10F10.4)') (RJKAB(I1,J1),J1=1,NMOS)
         write(MIOB,'(///10X,''TWO-ELECTRON K-INTEGRALS'',/)')
         do 130 I1=1,NMOS
  130    write(MIOB,'(10F10.4)') (RJKAB(I1,J1)-RJKAA(I1,J1),J1=1,NMOS)
      end if
      NATOMS=NUMAT
      do 140 I=1,NMOS
         do 140 J=1,NMOS
            RJKAA(I,J)=RJKAA(I,J)*0.5D0
  140 continue
      if (FIRST1) then
         I=INDEX(KEYWRD,'MICROS')
         if (I.NE.0) then
            K = READA(KEYWRD,I)
            LAB=K
            if (PRNT) write(MIOB,'(''    MICROSTATES READ IN'')')
            NTOT=XX+0.5
            rewind MIOA
            do 150 I=1,1000
               read(MIOA,'(A)') LINE
  150       if (INDEX(LINE,'MICRO').NE.0) go to 160
  160       do 170 I=1,1000
               read(MIOA,'(A)') LINE
  170       if (INDEX(LINE,'MICRO').NE.0) go to 180
  180       do 210 I=1,LAB
               read(MIOA,'(A)') LINE
               IZERO=MAX(0,MIN(INDEX(LINE,'0'),INDEX(LINE,'1'))-1)
               do 190 J=1,NMOS
                  if (LINE(J+IZERO:J+IZERO).NE.'1')
     1               LINE(J+IZERO:J+IZERO)='0'
                  if (LINE(J+NMOS+IZERO:J+NMOS+IZERO).NE.'1')
     1               LINE(J+NMOS+IZERO:J+NMOS+IZERO)='0'
                  MICROA(J,I)=ICHAR(LINE(J+IZERO:J+IZERO))-
     1          ICHAR('0')
                  MICROB(J,I)=ICHAR(LINE(J+NMOS+IZERO:J+NMOS+IZERO))-
     1          ICHAR('0')
  190          continue
               if (PRNT) write(MIOB,'(20I6)') (MICROA(J,I),J=1,NMOS),
     1                                       (MICROB(J,I),J=1,NMOS)
               K=0
               do 200 J=1,NMOS
  200          K=K+MICROA(J,I)+MICROB(J,I)
               if (K.NE.NTOT) then
                  NTOT=K
                  XX=K
                  write(MIOB,
     %              '(/,''NUMBER OF ELECTRONS IN C.I. REDEFINED TO'',
     %                  '':'',I4,/)') K
               end if
  210       continue
            FIRST1=.FALSE.
            go to 260
         end if
         NUPP=(NE+1)/2 +MSDEL
         NDOWN=NE-NUPP
         AMS=(NUPP-NDOWN)*0.5D0
         if (PRNT) write(MIOB,220) AMS
  220    format(10X,'COMPONENT OF SPIN  = ',F4.1)
         if (NUPP*NDOWN.LT.0) then
            write(MIOB,'(/10X,''IMPOSSIBLE VALUE OF DELTA S'')')
            call filclo !  J/M
            stop ' in MECI - IMPOSSIBLE VALUE OF DELTA''S'
         end if
         LIMA=NFA(NMOS+1)/(NFA(NUPP+1)*NFA(NMOS-NUPP+1))
         LIMB=NFA(NMOS+1)/(NFA(NDOWN+1)*NFA(NMOS-NDOWN+1))
         LAB=LIMA*LIMB
         if (PRNT) write(MIOB,230) LAB
  230    format(//10X,35H NO OF CONFIGURATIONS CONSIDERED = ,I4)
C#      if (LAB.LT.101) go to 240
C#      write(MIOB,230)
C#  230 format(10X,24H TOO MANY CONFIGURATIONS/)
C#      go to 160
C#  240 continue
         call PERM(NPERMA, NUPP, NMOS, NMECI, LIMA)
         call PERM(NPERMB, NDOWN, NMOS, NMECI, LIMB)
         K=0
         do 240 I=1,LIMA
            do 240 J=1,LIMB
               K=K+1
               do 240 L=1,NMOS
                  MICROA(L,K)=NPERMA(L,I)
  240    MICROB(L,K)=NPERMB(L,J)
  250    format(10I1)
  260    continue
         LIMA=LAB
         LIMB=LAB
      end if
      GSE=0.0D0
      do 270 I=1,NMOS
         GSE=GSE+EIGA(I)*OCCA(I)*2.D0
         GSE=GSE+XY(I,I,I,I)*OCCA(I)*OCCA(I)
         do 270 J=I+1,NMOS
  270 GSE=GSE+2.D0*(2.D0*XY(I,I,J,J) - XY(I,J,I,J))*OCCA(I)*OCCA(J)
      J=0
      do 280 I=1,LAB
         DIAG(I)=DIAGI(MICROA(1,I),MICROB(1,I),EIGA,XY,NMOS)-GSE
  280 continue
  290 continue
      if (LAB.LE.MDIM) go to 330
      X=-100.D0
      do 300 I=1,LAB
         if (DIAG(I).GT.X) then
            X=DIAG(I)
            J=I
         end if
  300 continue
      if (J.NE.LAB) then
         do 320 I=J,LAB
            I1=I+1
            do 310 K=1,NMOS
               MICROA(K,I)=MICROA(K,I1)
  310       MICROB(K,I)=MICROB(K,I1)
  320    DIAG(I)=DIAG(I1)
      end if
      LAB=LAB-1
      go to 290
  330 continue
C
C     BUILD SPIN AND NUMBER OF ALPHA SPIN TABLES.
C     -------------------------------------------
      do 350 I=1,LAB
         K=0
         X=0.D0
         do 340 J=1,NMOS
            X=X+MICROA(J,I)*MICROB(J,I)
  340    K=K+MICROA(J,I)
         NALPHA(I)=K
  350 SPIN(I)=4.D0*X-(XX-2*NALPHA(I))**2
C
C   BEFORE STARTING, CHECK THAT THE ROOT WANTED CAN EXIST
C
      if (LAB.LT.LROOT) then
         write(MIOB,
     %     '(//10X,''C.I. IS OF SIZE LESS THAN ROOT SPECIFIED'')')
         write(MIOB,'(10X,''MODifY SIZE OF C.I. OR ROOT NUMBER'')')
         call filclo !  J/M
         stop ' in MECI - check C.I. size or root #'
      end if
      if (PRNT) then
         write(MIOB,'(/,'' CONFIGURATIONS CONSIDERED IN C.I.      '',/
     %                  '' M.O. NUMBER :      '',10I4)')
     %                                        (I,I=NELEC+1,NELEC+NMOS)
         write(MIOB,'(''          ENERGY'')')
         do 360 I=1,LAB
            write(MIOB,'(/10X,I4,6X,10I4)') I,(MICROA(K,I),K=1,NMOS)
  360    write(MIOB,'(6X,F10.4,4X,10I4)')
     %                                DIAG(I),(MICROB(K,I),K=1,NMOS)
      end if
      call MECIH(DIAG,CIMAT,NMOS,LAB)
      if (BIGPRT) then
         write(MIOB,'(//,'' C.I. MATRIX'')')
         I=MIN(LAB,MAXORB)
         if (I.NE.LAB)
     1       write(MIOB,'(''   (OUTPUT HAS BEEN TRUNCATED)'')')
         call VECPRT(CIMAT,-I)
      else
         if (PRNT) write(MIOB,'(//,'' DIAGONAL OF C.I. MATRIX'')')
         if (PRNT) write(MIOB,'(5F13.6)') (CIMAT((I*(I+1))/2),I=1,LAB)
      end if
C#       call TIMER('SEC. DET. CONSTRUCTED')
      LABSIZ=MIN(LAB,LROOT+10)
      call HQRII(CIMAT,LAB,LABSIZ,EIG,CONF)
C#       call TIMER('DIAG. DONE')
C
C   DECIDE WHICH ROOT TO EXTRACT
C
      KROOT=0
      if (SMULT.LT.-0.1D0) then
         MECI=EIG(LROOT)
         do 370 J=1,LAB
  370    VECTCI(J)=CONF(J+LAB*(LROOT-1))
         KROOT=LROOT
      end if
      if (BIGPRT)  then
         write(MIOB,'(//20X,''STATE VECTORS'',//)')
         I=MIN(LAB,NORBS)
         J=MIN(LABSIZ,NORBS)
         call MATOUT(CONF,EIG,J,I,LAB)
      end if
      if (PRNT) then
         write(MIOB,380)
  380    format(///,' STATE ENERGIES '
     1,' EXPECTATION VALUE OF S**2  S FROM S**2=S(S+1)',//)
      end if
      IROOT=0
      do 390 I=1,9
  390 CIMAT(I)=0.1D0
      do 440 I=1,LABSIZ
         X=0.5D0*XX
         II=(I-1)*LAB
         do 420 J=1,LAB
            JI=J+II
            X=X-CONF(JI)*CONF(JI)*SPIN(J)*0.25D0
            K=ISPQR(J,1)
            if (K.EQ.1)  go to  410
            do 400 L=2,K
               LI=ISPQR(J,L)+II
  400       X=X+CONF(JI)*CONF(LI)*2.D0
  410       continue
  420    continue
         Y=(-1.D0+SQRT(1.D0+4.D0*X))*0.5D0
         if (ABS(SMULT-X).LT.0.01) then
            IROOT=IROOT+1
            if (IROOT.EQ.LROOT) then
               KROOT=I
               MECI=EIG(I)
               do 430 J=1,LAB
  430          VECTCI(J)=CONF(J+LAB*(I-1))
            end if
         end if
         J=Y*2.D0+1.5D0
         CIMAT(J)=CIMAT(J)+1
  440 if (PRNT) write(MIOB,460) I,EIG(I),TSPIN(J),X,Y
      if (KROOT.EQ.0) then
         write(MIOB,
     %     '(//10X,''THE STATE REQUIRED IS NOT PRESENT IN THE'')')
         write(MIOB,'(10X,  ''    SET OF CONFIGURATIONS AVAILABLE'')')
         write(MIOB,
     %    '(/ 4X,''NUMBER OF STATES ACCESSIBLE USING CURRENT KEY-'',
     %           ''WORDS'',/)')
         do 450 I=1,7
  450    if (CIMAT(I).GT.0.5D0)
     1     write(MIOB,'((24X,A8,I4))')   TSPIN(I),NINT(CIMAT(I))
         call filclo !  J/M
         stop ' in MECI - state required is not present in conf. avail.'
      end if
  460 format(I5,F12.6,3X,A8,F15.5,F10.5)
  470 continue
      MAXVEC=0
      if (LSPIN) MAXVEC=MIN(4,LAB)
      if (LSPIN.AND.(NE/2)*2.EQ.NE) then
         write(MIOB,
     %     '(''   ESR SPECIFIED FOR AN EVEN-ELECTRON SYSTEM'')')
      end if
C#      do 570 I=1,NMOS
C#         do 570 J=1,NORBS
C#  570 COEFF(J,I+NELEC)=COEFF(J,I+NELEC)**2
      do 540 IUJ=1,MAXVEC
         IOFSET=(IUJ-1)*LAB
         write(MIOB,'(//,''      MICROSTATE CONTRIBUTIONS TO '',
     1                   ''STATE EIGENFUNCTION'',I3)')  IUJ
         write(MIOB,'(5F13.6)') (CONF(I+IOFSET),I=1,LAB)
         do 480 I=1,LAB
  480    CONF(I)=CONF(I+IOFSET)**2
C                                             SECOND VECTOR!
         do 500 I=1,NMOS
            SUM=0.D0
            do 490 J=1,LAB
  490       SUM=SUM+(MICROA(I,J)-MICROB(I,J))*CONF(J)
  500    EIGA(I)=SUM
         write(MIOB,'(/,''    SPIN DENSITIES FROM EACH M.O., ENERGY:''
     %         ,F7.3)') EIG(IUJ)
         write(MIOB,'(5F12.6)') (EIGA(I),I=1,NMOS)
         write(MIOB,*)
         write(MIOB,*) '     SPIN DENSITIES FROM EACH ATOMIC ORBITAL'
         write(MIOB,*)
     %     '                              S        PX        '//
     %     'PY        PZ        TOTAL'
         do 530 I=1,NATOMS
            IL=NFIRST(I)
            IU=NLAST(I)
            L=0
            SUMM=0.D0
            do 520 K=IL,IU
               L=L+1
               SUM=0.D0
               do 510 J=1,NMOS
  510          SUM=SUM+COEFF(K,J+NELEC)*EIGA(J)
               SUMM=SUMM+SUM
  520       EIGS(L)=SUM
            if (L.EQ.4) then
               write(MIOB,
     %           '(''  ATOM'',I4,''    SPIN DENSITY  '',5F10.7)')
     %                                           I,(EIGS(K),K=1,L),SUMM
            else
               write(MIOB,
     %           '(''  ATOM'',I4,''    SPIN DENSITY  '',F10.7,30X,
     %             F10.7)') I,EIGS(1),SUMM
            end if
  530    continue
  540 continue
      return
      end
