C 14 NOV 05 - DGF - PAD COMMON BLOCK SCFOPT
C 19 SEP 05 - MWS - ADJUSTMENT FOR H OR I FUNCTION RUNS
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 30 APR 05 - DGF - REPAD COMMON BLOCK DMPING
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - PAD COMMON BLOCK NSHEL, FRGINF
C  5 FEB 05 - MWS - PAD COMMON BLOCK DMPING
C  7 SEP 04 - MWS - PAD COMMON BLOCK INTFIL
C 23 JUL 04 - MWS - RENAME -PRINTM- TO -PRINTSCFMI-
C  7 APR 04 - MS  - INCLUDE MULTIFRAGMENT IMPLEMENTATION OF SCF-MI
C
C      ****************************************************
C      *                     SCF-MI                       *
C      * SELF CONSISTENT FIELD FOR MOLECULAR INTERACTIONS *
C      ****************************************************
C               (DIRECT AND CONVENTIONAL)
C     ----- CLOSED SHELL SCF-MI CALCULATION -----
C      E. GIANINETTI, M. RAIMONDI, E. TORNAGHI
C      INT. J. QUANTUM CHEM. 60, 157-166(1996)
C
C THIS MULTI-FRAGMENT IMPLEMENTATION FROM MARIUZIO SIRONI
C REPLACES A 2-FRAGMENT CODE THAT WAS INCLUDED WITH GAMESS
C BY NINO FAMULARI IN OCTOBER 1997.
C
C*MODULE SCFMI   *DECK MIINP
      SUBROUTINE MIINP
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,DIISON
C
      PARAMETER (MXATM=2000, MAXNFRAGS=10)
C
      COMMON /BASSPH/ QMTTOL,ISPHER
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NAINF,NBINF,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJKAO,IJKMO,IDAF,NAV,IODA(950)
      COMMON /MIOPT / DTOL,CNVLOC,DIISTL,ALPHA,
     *                NF(MAXNFRAGS),MF(MAXNFRAGS),NFRAGS,MSHIFT,ITER,
     *                ITLOC,IGUESS,MXDIIS,
     *                IOPT,MIFLAG,DIISON
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
C     ----- SET UP NAMELIST SIMULATION FOR $SCFMI-----
C
      PARAMETER (NNAM=14)
      DIMENSION QNAM(NNAM), KQNAM(NNAM)
      DATA SCFMINP/8HSCFMI   /
      DATA QNAM/8HNFRAGS  ,8HNF      ,8HMF      ,8HMSHIFT  ,
     *          8HITER    ,8HDTOL    ,8HITLOC   ,8HCNVLOC  ,
     *          8HDIISON  ,8HMXDIIS  ,8HDIISTL  ,8HALPHA   ,
     *          8HIGUESS  ,8HIOPT     /
      DATA KQNAM/1,-1,-1,1,1,3,1,3,0,1,3,3,1,1/
C
      DATA RNONE,RHF/8HNONE    ,8HRHF     /
C
C        ---- READ POSSIBLE $SCFMI INPUT GROUP -----
C
      NFRAGS = 0
      KQNAM(2)=MAXNFRAGS*10 + 1
      DO 5 I=1,MAXNFRAGS
         NF(I) = 0
    5 CONTINUE
      KQNAM(3)=MAXNFRAGS*10 + 1
      DO 6 I=1,MAXNFRAGS
         MF(I) = 0
    6 CONTINUE
C
C        MSHIFT IS AN UNDOCUMENTED "MEMORY SHIFT", WHICH MIGHT
C        LET A DEBUGGER APPLY A POSITIVE SHIFT TO ALL MEMORY
C        POINTERS.  NOT REALLY USEFUL TO ANYONE ELSE.
C
      MSHIFT = 0
      ITER = 50
      DTOL = 1.0D-6
      ITLOC = 50
      CNVLOC = 0.02D+00
      IGUESS = 0
      DIISON = .TRUE.
      MXDIIS = 10
      DIISTL = 0.01D+00
      ALPHA = 0.0D+00
      IOPT = 0
C
      JRET=0
      CALL NAMEIO(IR,JRET,SCFMINP,NNAM,QNAM,KQNAM,
     *            NFRAGS,NF,MF,MSHIFT,ITER,DTOL,ITLOC,CNVLOC,DIISON,
     *     MXDIIS,DIISTL,ALPHA,IGUESS,IOPT,  0,0,0,0,0,  0,0,0,0,0,
     *     0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *     0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
C
C         IF NOT FOUND IN INPUT, WE ARE NOT DOING SCF-MI CALCULATION
C
      MIFLAG=0
      IF(JRET.EQ.1) RETURN
C
      MIFLAG=1
C
C         PRINT $SCFMI PARAMETERS
C
      IF(MASWRK) WRITE(IW,*)
      IF(MASWRK) WRITE(IW,*)'        **********************'
      IF(MASWRK) WRITE(IW,*)'        SCFMI INPUT PARAMETERS'
      IF(MASWRK) WRITE(IW,*)'        **********************'
      IF(MASWRK) WRITE(IW,*)
C
  846 FORMAT(1X,'FRAGMENTS: ',10(2X,I3))
      IF(MASWRK) WRITE(IW,846) (III,III=1,NFRAGS)
      IF(MASWRK) WRITE(IW,847) (NF(III),III=1,NFRAGS)
      IF(MASWRK) WRITE(IW,848) (MF(III),III=1,NFRAGS)
      IF(MASWRK) WRITE(IW,*)
      IF(MASWRK) WRITE(IW,*) ' NFRAGS =  ',NFRAGS
      IF(MASWRK) WRITE(IW,*) ' ITER =    ',ITER
      IF(MASWRK) WRITE(IW,*) ' DTOL =    ',DTOL
      IF(MASWRK) WRITE(IW,*) ' DIISON =  ',DIISON
      IF(MASWRK) WRITE(IW,*) ' MXDIIS = ',MXDIIS
      IF(MASWRK) WRITE(IW,*) ' DIISTL = ',DIISTL
      IF(MASWRK) WRITE(IW,*) ' ITLOC =   ',ITLOC
      IF(MASWRK) WRITE(IW,*) ' CNVLOC = ',CNVLOC
      IF(MASWRK) WRITE(IW,*) ' IOPT =    ',IOPT
      IF(MASWRK) WRITE(IW,*) ' ALPHA =    ',ALPHA
      IF(MASWRK) WRITE(IW,*)
C
C     IF(IGUESS.EQ.-99)THEN
C       IF(MASWRK) WRITE(IW,*) 'IGUESS =  ',IGUESS
C       IF(MASWRK) WRITE(IW,*)'SCF-MI ORBITALS WILL BE',
C    +             ' WRITTEN ON UNIT 30'
C     ENDIF
C
C     IF(IGUESS.EQ.-999)THEN
C       IF(MASWRK) WRITE(IW,*) 'IGUESS =  ',IGUESS
C       IF(MASWRK) WRITE(IW,*)'SCF-MI GUESS ORBITALS WILL BE',
C    +            ' READED FROM UNIT 31'
C       IF(MASWRK) WRITE(IW,*)'SCF-MI ORBITALS WILL BE WRITTEN',
C    +            ' ON UNIT 30'
C     ENDIF
C
      IF(MASWRK) WRITE(IW,*)
C
      NERR=0
      IF(JRET.EQ.2) NERR=NERR+1
C
      IF(ISPHER.EQ.1) THEN
         IF(MASWRK) WRITE(IW,*)
     *       'ERROR - SCF-MI CANNOT USE SPHERICAL HARMONICS'
         NERR=NERR+1
      END IF
C
      IF(NFRAGS.GT.MAXNFRAGS)THEN
         IF(MASWRK) WRITE(IW,*) 'ERROR - MAXIMUM NFRAGS IS',MAXNFRAGS
         NERR=NERR+1
      END IF
C
      IF(NFRAGS.LT.2)THEN
         IF(MASWRK) WRITE(IW,*) 'ERROR - NFRAGS MUST BE AT LEAST 2'
         NERR=NERR+1
      END IF
C
      NTOTALI=0
      MTOTALI=0
C
      DO III=1,NFRAGS
      IF(NF(III).LE.0 ) THEN
         IF(MASWRK) WRITE(IW,*) 'ERROR - NF(',III,') MUST BE',
     +             ' POSITIVE VALUES'
         NERR=NERR+1
      END IF
C
      IF(MF(III).LT.NF(III) ) THEN
         IF(MASWRK) WRITE(IW,*) 'ERROR - MF(',III,') MUST',
     +             ' EXCEED NF(',III,')'
         NERR=NERR+1
      END IF
C
      NTOTALI=NTOTALI+NF(III)
      MTOTALI=MTOTALI+MF(III)
C
      ENDDO
C
      IF(NTOTALI  .NE.  NAINF) THEN
        IF(MASWRK) WRITE(IW,*) 'ERROR - SUMMATION OF NF MUST EQUAL NO.',
     *               ' OCC. MO-S =',NAINF
         NERR=NERR+1
      END IF
C
      IF(MTOTALI.NE.  NUM) THEN
        IF(MASWRK) WRITE(IW,*) 'ERROR - SUMMATION OF MF MUST EQUAL',
     +             ' TOTAL NO. AO-S =',NUM
         NERR=NERR+1
      END IF
C
      IF(SCFTYP.NE.RHF) THEN
         IF(MASWRK) WRITE(IW,*) 'ERROR - SCFMI WORKS ONLY WITH',
     +             ' -RHF- WAVEFUNCTIONS'
         NERR=NERR+1
      END IF
C
      IF(MPLEVL.NE.0  .OR.  CITYP.NE.RNONE  .OR.
     *   DFTYPE.NE.RNONE  .OR.  CCTYP.NE.RNONE) THEN
         IF(MASWRK) WRITE(IW,*) 'ERROR - SCFMI DOES NOT PERMIT',
     +              ' ANY ELECTRON CORRELATION METHOD.'
         NERR=NERR+1
      END IF
C
C     IF(NPROC.GT.1) THEN
C        WRITE(IW,*) 'ERROR - SCFMI DOES NOT RUN IN PARALLEL'
C        NERR=NERR+1
C     END IF
C
  847 FORMAT(1X,'NR. MO',5X,10(2X,I3))
  848 FORMAT(1X,'NR. AO',5X,10(2X,I3))
C
      IF(NERR.GT.0) THEN
         IF(MASWRK) WRITE(IW,*) 'ERROR IN $SCFMI INPUT, JOB',
     +              ' IS TERMINATING.'
         CALL ABRT
      END IF
C
      RETURN
      END
C*MODULE SCFMI   *DECK SCFMI
      SUBROUTINE SCFMI
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DIISON
      LOGICAL PACK2E,DIRSCF,FDIFF
      LOGICAL DODIIS,DEM
      LOGICAL GOPARR,DSKWRK,MASWRK
      LOGICAL ABINIT,MINMEM,IZRFLS
      LOGICAL SOSCF
C
C...  MAXNFRAGS -> MAX NUMERO DI FRAMMENTI
      PARAMETER (MAXNFRAGS=10)
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30, MXG2=MXGSH*MXGSH,
     *           MXATM=2000, MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /DIISSO/ SOGTOL,ETHRSH,MAXDII,IRAF
      COMMON /FMCOM / X(1)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /INDDIP/ TMUX,TMUY,TMUZ,TMUXD,TMUYD,TMUZD,MINMEM
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /IOFILE/ IR,IW,IP,IJKAO,IJKMO,IDAF,NAV,IODA(950)
      COMMON /N2ELCT/ N2EL
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
      COMMON /MIOPT / DTOL,CNVLOC,DIISTL,ALPHA,
     *                NF(MAXNFRAGS),MF(MAXNFRAGS),NFRAGS,MSHIFT,ITER,
     *                ITLOC,IGUESS,MXDIIS,
     *                IOPT,MIFLAG,DIISON
C
C...  MF -> VETTORE CONTENENTE IL PARTIZIONAMENTO DELLA BASE
C...  NF -> VETTORE CONTENENTE IL PARTIZIONAMENTO DEGLI ELETTRONI
C
C--   COMMON /MIADDI/ I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12,I13,
C--  + I14,I15,I16,I17,I18,I19,I20,I21,I22,I23,I24,I25,I26,I27,I28,I29,
C--  + I30,I31,I32,I33,I34,I35,I36,I37,I38,I39,I40,I41,I42,I43,I44,I45,
C--  + I46,I47,I48,I49,I50,I51,I52,I53,I54,I55,I56,I57,I58,I59,I60,I61,
C--  + I62,I63,I64,I65,I66,I67,I68,I69,I70,I71,I72,I73,I74,I75,I76,I77,
C--  + I78,I79,I80,I81,I82,I83,I84,I85,I86,I87,I88,I89,I90
C
      DATA NONE/4HNONE/
C
      MMFTOT=0
      NNFTOT=0
C
      DO III=1,NFRAGS
        MMFTOT=MMFTOT+MF(III)
        NNFTOT=NNFTOT+NF(III)
      ENDDO
C
C DEFAULTS FOR NAMELIST INPUT
C
      MSHIFT=0
      IBF=5
      SOGLIA=CNVLOC
      NITM=ITLOC
      MAXDIIS=MXDIIS
      DIISTOL=DIISTL
C...
C     DEFINIZIONE DI M ED N (TOTALI)
C...
      M=0
      N=0
      DO I=1,NFRAGS
         M=M+MF(I)
         N=N+NF(I)
      END DO
C
C
      IBF=MAX0(1,IBF)*1000
C
      MTR=M*(M+1)/2
C
C...
C     DEFINIZIONE DI MMAS
C...
      MMAS=0
      DO I=1,NFRAGS
         IDUM=MAX0(MMAS,MF(I))
         IF (IDUM.GT.MMAS) THEN
            MMAS=IDUM
         END IF
      END DO
C
      IF(MASWRK) WRITE(IW,6666)
 6666 FORMAT(/5X,60('*')/
     *        5X,'SCF FOR MOLECULAR INTERACTIONS (SCF-MI)',7X,
     *           'N-BODY VERSION'/
     *       16X,'E.GIANINETTI, M.RAIMONDI, E.TORNAGHI'/
     *       15X,'INT. J. QUANTUM CHEM. 60, 157-166(1996)'/
     *        5X,60('*')/)
C
C     ----- SET POINTERS FOR PARTITIONING OF THE MEMORY -----
C
C READ THE LAST OCCUPIED POSITION IN X(1)
      CALL VALFM(LOADFM)
      IF(IOPT.GT.1) THEN
         IF(MASWRK) WRITE(6,*)'START FIRST BLOCK --> NEED1'
         IF(MASWRK) WRITE(6,*)'LOADFM -->',LOADFM
      END IF
C
      IDBL=1
      I1 = LOADFM + 1 + MSHIFT
C     IV(I1) ==> T
      I2 =I1  + N*M*IDBL + MSHIFT
C     IV(I2) ==> SAO
      I3 =I2  + M*M*IDBL + MSHIFT
C     IV(I3) ==> HAO
      I4 =I3  + MTR*IDBL + MSHIFT
C     IV(I4) ==> P
      I5 =I4  + MTR*IDBL + MSHIFT
      I6 =I5
      I7 =I6
      I8 =I7
      I9 =I8
C     IV(I9) ==> SOCC
      I10=I9  + N*N*IDBL + MSHIFT
C     IV(I10) ==> SWORK
      I11=I10 + MAX0(M*M,N*N,N*100)*IDBL + MSHIFT
C     IV(I11) ==> TTOT
      I12=I11 + M*M*IDBL + MSHIFT
      I13=I12
C     WRITE(6,*)'I13', I13
C
C
      I14=I13
      I15=I14
      I16=I15
      I17=I16
      I18=I17
C--      I19=I18
C--C     IV(I19) ==> MAP
C--      I20=I19 + M + MSHIFT
C--C     IV(I20) ==> INDXMAP
C--      I21=I20 + MTR + MSHIFT
      I21=I18
C     IV(I21) ==> SSQR
      I22=I21 + M*M*IDBL + MSHIFT
C     IV(I22) ==> PREV
      I23=I22 + MTR*IDBL + MSHIFT
      I24=I23 + MSHIFT
      I25=I24
      I26=I25
C     IV(I26) ==> TOTG
      I27=I26 + M*M*IDBL + MSHIFT
C     IV(I27) ==> ROOT
C--      I28=I27 + MAX0(MMAS,N)*IDBL + MSHIFT
C--C     IV(I28) ==> RMAT
C--      I29=I28 + MMAS*M*IDBL + MSHIFT
      I29=I27 + MAX0(MMAS,N)*IDBL + MSHIFT
C     IV(I29) ==> ZMAT
      I30=I29 + MMAS*MMAS*IDBL + MSHIFT
      I31=I30 + MSHIFT
      I32=I31 + MSHIFT
C     IV(I32) ==> U
      I33=I32 + MAX0(MTR,MMAS*(MMAS+2))*IDBL + MSHIFT
      I34=I33 + MSHIFT
      I35=I34 + MSHIFT
      I36=I35 + MSHIFT
      I37=I36 + MSHIFT
      IPUNT=0
      DO I=1,NFRAGS
         IPUNT=IPUNT+MF(I)*MF(I)
      END DO
C     IV(I37) ==> FFI(GMAT1, GMAT2,....,GMATN)
      I38=I37 + IPUNT*IDBL + MSHIFT
C     IV(I38) ==> SSI(SG1, SG2,....,SGN)
      I39=I38 + IPUNT*IDBL + MSHIFT
C     IV(I39) ==> TVEC(T1, T2,....,TN)
      I40=I39 + IPUNT*IDBL + MSHIFT
C     IV(I40) ==> DMAT(M,M)
      I41=I40 +  M*M*IDBL + MSHIFT
C     IV(I41) ==> DS(M,M)
C--      I42=I41 +  M*M*IDBL + MSHIFT
C--C     IV(I42) ==> TV(M,N)
C--      I43=I42 +  M*N*IDBL + MSHIFT
      I43=I41 +  M*M*IDBL + MSHIFT
C     IV(I43) ==>
      I44=I43
C     IV(I44) ==>WORK2
      I45=I44 + M*(M+1)/2*IDBL +MSHIFT
C     IV(I45) ==> VTTS(N,M)
      I46=I45 +  N*M*IDBL + MSHIFT
C     IV(I46) ==> TVDBVTTS(M,M)
      I47=I46 +  M*M*IDBL + MSHIFT
C     IV(I47) ==> PFOCK(M,M)
      I48=I47 +  M*M*IDBL + MSHIFT
C     IV(I48) ==> FP(M,M)
      I49=I48 +  M*M*IDBL + MSHIFT
C     IV(I49) ==>
      I50=I49
C     IV(I50) ==> FEFF(M,M)
      I51=I50 +  M*M*IDBL + MSHIFT
C     IV(I51) ==> DBS(M,M)
      I52=I51 +  M*M*IDBL + MSHIFT
C     IV(I52) ==>
      I53=I52
C     IV(I53) ==> STVDOTVTTS(M,M)
      I54=I53 +  M*M*IDBL + MSHIFT
C     IV(I54) ==> SDS(M,M)
      I55=I54 +  M*M*IDBL + MSHIFT
C     IV(I55) ==> SEFF(M,M)
      I56=I55 +  M*M*IDBL + MSHIFT
C     IV(I56) ==> OLDFFI
      I57=I56 + IPUNT*IDBL*MAXDIIS+MSHIFT
C     IV(I57) ==> OLDSSI
      I58=I57 + IPUNT*IDBL*MAXDIIS+MSHIFT
C     IV(I58) ==> OLDERR
      I59=I58  + IPUNT*IDBL*MAXDIIS+MSHIFT
C     IV(I59) ==>  BMAT
      I60=I59 + (MAXDIIS+1)**2 *IDBL +MSHIFT
C
C     IV(I60) ==> EIGEN
      I61=I60 + M*IDBL + MSHIFT
C     IV(I61) ==> EIGENWORK
      I62=I61 + M*IDBL + MSHIFT
C     IV(I62) ==> INDX
      I63=I62 + M*IDBL + MSHIFT
C     IV(I63) ==> TOTWORK
      I64=I63 + M*M*IDBL + MSHIFT
C     IV(I64) ==> FDUM
      I65=I64 + M*M*IDBL + MSHIFT
C     IV(I65) ==> FMO
      I66=I65 + N*N*IDBL + MSHIFT
C     IV(I66) ==> SMO
      I67=I66 + N*N*IDBL + MSHIFT
C     IV(I67) ==> EMO
      I68=I67 + N*IDBL + MSHIFT
C     IV(I68) ==> TMO
      I69=I68 + N*(N+2)*IDBL + MSHIFT
C     IV(I69) ==> TNEW
      I70=I69 + M*M*IDBL + MSHIFT
C     IV(I70) ==> DUM
      I71=I70 + MTR*IDBL + MSHIFT
C     IV(I71) ==> WORK
C     I72=I71 + N*N*IDBL + MSHIFT
      I72=I71 + MAX0(N*N,9*N)*IDBL + MSHIFT
C     IV(I72) ==> SINV
C--      I73=I72 + N*N*IDBL + MSHIFT
C--C     IV(I73) ==> OVER
C--      I74=I73 + M*M*IDBL + MSHIFT
      I74=I72 + N*N*IDBL + MSHIFT
C     IV(I74) ==> RIJ
      I75=I74 + N*(N+1)*NFRAGS*IDBL/2  + MSHIFT
C     IV(I75) ==> ICONT
      I76=I75 + NFRAGS + MSHIFT
C     IV(I76) ==> IORD
      I77=I76 + M + MSHIFT
C     IV(I77) ==> IIR
      I78=I77 + M   + MSHIFT
C     IV(I78) ==> MAPP
C--      I79=I78 + M + MSHIFT
C--C     IV(I79) ==> CERR
C--      I80=I79 +(MAXDIIS+1)*IDBL +MSHIFT
C--C     IV(I80) ==> A
C--      I81=I80 + M*IDBL + MSHIFT
C--C     IV(I81) ==> TOTAO
C--      I82=I81 + M*IDBL   + MSHIFT
      I82=I78 + M   + MSHIFT
C     IV(I82) ==> TRAN
      I83=I82 + N * N * IDBL + MSHIFT
C     IV(I83) ==>  AOMULL
      I84=I83 + M*M*IDBL + MSHIFT
C     IV(I84) ==> NOROT
      I85=I84 + N * N + MSHIFT
C     IV(I85) ==> SOCC1
      I86=I85 + N*(N+1)/2*IDBL + MSHIFT
C     IV(I86) ==> QPIX
      I87=I86 + M
C     IV(I87) ==> QPJX
      I88=I87 + M
C     IV(I88) ==> SAOTR
      I89=I88 + M*(M+1)*IDBL/2
C     IV(I89) ==> GROSS
      I90=I89 + NFRAGS * N * IDBL
C
      NEED2=I90-LOADFM-1
C
C ***** FROM GAMESS
C
C
      ABINIT = MPCTYP.EQ.NONE
C
C         DIFFERENTIAL FOCK FORMATION IS NOT SUPPORTED
C
      IF(DIRSCF) FDIFF=.FALSE.
C
      IZRFLS = IZRF.NE.0
C
C     ----- SELECT CONVERGENCE METHOD(S) -----
C
      DODIIS = MOD(MCONV,32) .GE.16
      DEM    = MOD(MCONV,64) .GE.32
      SOSCF  = MOD(MCONV,128).GE.64
C
C     ----- SET POINTERS FOR PARTITIONING OF MEMORY -----
C
      LL1 = NUM
      LL2 = (LL1*LL1+LL1)/2
      LL3 = LL1*LL1
C
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
C
      MAXIT2=(MAXIT*MAXIT+MAXIT)/2
      MAXIT3=MAXIT*MAXIT
C
      NSH2 = (NSHELL*NSHELL+NSHELL)/2
C
C
      LSCR  = I90
      LIWRK = LSCR   + L1*8
      LWRK  = LIWRK  + L1
      LWRK1 = LWRK   + L1
      LWRK2 = LWRK1  + LL3
      LWRK3 = LWRK2  + LL3
      LEIG  = LWRK3  + LL3
      LOCC  = LEIG   + LL1
      LFAO  = LOCC   + L1
      LD    = LFAO   + LL2
      LVEC  = LD     + LL2
      LAST  = LVEC   + LL3
C
      CALL BASCHK(LMAX)
                    NANGM =  4
      IF(LMAX.EQ.2) NANGM =  6
      IF(LMAX.EQ.3) NANGM = 10
      IF(LMAX.EQ.4) NANGM = 15
      IF(LMAX.EQ.5) NANGM = 21
      IF(LMAX.EQ.6) NANGM = 28
      MAXG = NANGM**4
C
      IF(ABINIT) THEN
         IF(DIRSCF) THEN
            LBUF  = LAST
            LIBUF = LAST
            LXINTS= LAST
            LGHOND= LXINTS + NSH2
            LDSH  = LGHOND + MAXG
            LDDIJ = LDSH   + NSH2
            LAST  = LDDIJ  + 49*MXG2
         ELSE
            LXINTS= LAST
            LGHOND= LAST
            LDSH  = LAST
            LDDIJ = LAST
            LBUF  = LAST
            LIBUF = LBUF   + NINTMX
            LAST  = LIBUF  + NINTMX
         END IF
         LTWOEI= LAST
         LPTOT2= LAST
      ELSE
         LXINTS= LAST
         LGHOND= LAST
         LDSH  = LAST
         LDDIJ = LAST
         LBUF  = LAST
         LIBUF = LAST
         LTWOEI= LAST
         LPTOT2= LTWOEI + N2EL
         LAST  = LPTOT2 + NAT*16
      END IF
C
      IF(DODIIS) THEN
         LADIIS = LAST
         LBDIIS = LADIIS + MAXIT3
         LPDIIS = LBDIIS + MAXIT2
         LXDIIS = LPDIIS + MAXIT
         LIODII = LXDIIS + MAXIT
         LAST   = LIODII + 2*MAXDII
      ELSE
         LADIIS = LAST
         LBDIIS = LAST
         LPDIIS = LAST
         LXDIIS = LAST
         LIODII = LAST
      END IF
C
      IF(DEM) THEN
         LWRK4 = LAST
         LWRK5 = LWRK4 + L3
         LAST  = LWRK5 + L3
      ELSE
         LWRK4 = LAST
         LWRK5 = LAST
      END IF
C
      IF (IZRFLS .OR. IEFP.EQ.1) THEN
         LXDINT = LAST
         LYDINT = LXDINT  + L2
         LZDINT = LYDINT  + L2
         LAST   = LZDINT  + L2
      ELSE
         LXDINT = LAST
         LYDINT = LAST
         LZDINT = LAST
      END IF
      IF((IZRFLS .OR. IEFP.EQ.1) .AND. (DIRSCF.AND.FDIFF)) THEN
         NDIPER = LAST
         LDIPER = NDIPER  + L2
         LAST   = LDIPER  + L2
      ELSE
         NDIPER = LAST
         LDIPER = LAST
      END IF
      IF(IEFP.EQ.1)THEN
         LEFLD  = LAST
         MADD   = LEFLD  + 3*NPTTPT
         LABFLD = MADD   + 3*NPTTPT
         LAST   = LABFLD + 3*NPTTPT
      ELSE
         LEFLD  = LAST
         MADD   = LAST
         LABFLD = LAST
      END IF
C
      NPR=(NUM-NA)*NA
      IF (SOSCF) THEN
         LGRAD  = LAST
         LHESS  = LGRAD  + NPR
         LPGRAD = LHESS  + NPR
         LDISPL = LPGRAD + NPR
         LAST   = LDISPL + NPR
      ELSE
         LGRAD  = LAST
         LHESS  = LAST
         LPGRAD = LAST
         LDISPL = LAST
      END IF
C
C        WE CAN REDUCE MEMORY NEED BY 2N**2 BY DOING MORE DAREADS,
C        TO OVERLAP -H1- AND -WRK1-, -Q- AND -WRK2-, -S- AND -WRK3-.
C
      NEED = LAST - LOADFM - 1
      MINMEM = .FALSE.
      LH1   = LAST
      LQ    = LH1    + LL2
      LS    = LQ     + LL3
      LAST  = LS     + LL2
      NEED  = LAST -LOADFM -1
C
C CHECK FOR TOTAL MEMORY AVAIBILITY
C
      CALL GETFM(NEED)
      CALL GOTFM(NGOTMX)
      CALL RETFM(NEED)
C
C GET THE BLOCK OF MEMORY NEED2
C
      CALL GETFM(NEED2)
      CALL GOTFM(NGOTMX)
      IF(IOPT.GT.1) THEN
        IF(MASWRK)  WRITE(6,*)'NEED2 --> ', NEED2
        IF(MASWRK)  WRITE(6,*)'NGOTMX --> ', NGOTMX
      END IF
C
      ILELE3=IPUNT
      CALL MIITER(X(I1),X(I21),X(I3),X(I4),X(I9),X(I10),X(I26),
     *            X(I27),X(I29),X(I32),X(I22),X(I11),ERN,DTOL,MMAS,
     *            ITER,N,M,MTR,IOPT,X(I37),X(I38),X(I39),MF,NF,NFRAGS,
     *            X(I60),X(I61),X(I62),X(I63),X(I64),X(I65),X(I66),
     *            X(I67),X(I68),X(I69),X(I70),X(I71),X(I72),ALPHA,
     *            X(I40),X(I41),X(I45),X(I46),X(I47),X(I48),X(I50),
     *            X(I51),X(I53),X(I54),X(I55),X(I74),X(I75),X(I76),
     *            X(I77),X(I78),X(I82),X(I83),X(I84),NITM,
     *            SOGLIA,IGUESS,X(I85),ILELE3,MAXDIIS,X(I56),X(I57),
     *            X(I58),X(I59),DIISTOL,DIISON,X(I44),X(I86),X(I87),
     *            X(I88),X(I89))
C
      CALL RETFM(NEED2)
      IF(MASWRK) WRITE(IW,9000)
      CALL TIMIT(1)
      RETURN
 9000 FORMAT(1X,'.... DONE WITH SCF-MI CALCULATION ....')
      END
C
C*MODULE SCFMI   *DECK MIINIT
      SUBROUTINE MIINIT(T,TTOT,RIJ,IORD,IIR,MAPP,TRAN,NOROT,
     *                  QPIX,QPJX,SAOTR,N,M,IGUESS,NF,MF,NFRAGS,
     *                  NAT,GROSS,NITM,SOGLIA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXAO=8192)
C
      DIMENSION T(M,N),TTOT(M,M),NOROT(N,N),NF(*),MF(*),MAPP(*),
     +          GROSS(NFRAGS,N),RIJ(*)
C
      COMMON /IOFILE/ IR,IW,IP,IJKAO,IJKMO,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
C
      DATA IFIRST/1/
      SAVE IFIRST
C
      IF(IFIRST.NE.1) GOTO 123
      IFIRST = 0
      IF (IGUESS.NE.-999) THEN
      IF(MASWRK) WRITE(IW,*)' A PROPER GUESS FOR THE SCFMI',
     +           ' WAVEFUNCTION IS BUILT UP'
         L1 = M
         L2 = M * (M+1)/2
         N1 = N
         N2 = N * (N+1)/2
         NOUT=0
         IPASS=0
         DO I=1,N
         DO J=1,N
           NOROT(J,I) = 0
         END DO
         END DO
C
         NATSAVE = NAT
         NAT = NFRAGS
       CVGLOCSAVE=CVGLOC
       CVGLOC=SOGLIA
       MAXLOCSAVE=MAXLOC
       MAXLOC=NITM
C
C
C
         ICONT = 0
         DO I=1,NFRAGS
            DO J=1,MF(I)
               ICONT = ICONT + 1
               MAPP(ICONT) = I
            END DO
         END DO
         CALL POPLOC2(IPASS,TTOT,TRAN,T,IORD,IIR,MAPP,RIJ,QPIX,QPJX,
     +               NOROT,SAOTR,L1,L2,N1,N2,NAT,NOUT)
         CALL POPQAT(RIJ,GROSS,N2,N,NAT)
C
         NAT = NATSAVE
       CVGLOC=CVGLOCSAVE
       MAXLOC=MAXLOCSAVE
C
C NOROT(I,1)---> FRAMMENTO CUI APPARTIENE L'ORBITALE I-ESIMO
C
         DO I=1,N
             CMAX = -999.0D+00
             DO IK=1,NFRAGS
             IF(GROSS(IK,I) .GT. CMAX ) THEN
                 IKMAX = IK
                 CMAX = GROSS(IK,I)
               END IF
             END DO
           NOROT(I,1)=IKMAX
         ENDDO
C
C ASSEGNA GLI ORBITALI A TTOT IN BASE A NOROT(I,1)
C
        CALL VCLR(TTOT,1,M*M)
C
        ICOUNT=0
        ISTART=0
        DO I=1,NFRAGS
          DO II=1,N
            IF(NOROT(II,1).EQ.I)THEN
              ICOUNT=ICOUNT+1
              DO J=ISTART+1,ISTART+MF(I)
                TTOT(J,ICOUNT)=T(J,II)
              ENDDO
            ENDIF
          ENDDO
C          ICOUNT=MF(I)
C          ISTART=MF(I)
          ICOUNT=ISTART+MF(I)
          ISTART=ISTART+MF(I)
        ENDDO
C
C COSTRUISCE T DA TTOT
C
        ICOUNT=0
        ISTART=0
        DO I=1,NFRAGS
          DO II=ISTART+1,ISTART+NF(I)
             ICOUNT=ICOUNT+1
             DO J=1,M
               T(J,ICOUNT)=TTOT(J,II)
             ENDDO
          ENDDO
          ISTART=ISTART+MF(I)
        ENDDO
C
C
      ELSE
C
      IF(MASWRK) REWIND 31
      IF(MASWRK) READ(31) TTOT
        ICOUNT=0
        ISTART=0
        DO I=1,NFRAGS
          DO II=ISTART+1,ISTART+NF(I)
             ICOUNT=ICOUNT+1
             DO J=1,M
               T(J,ICOUNT)=TTOT(J,II)
             ENDDO
          ENDDO
          ISTART=ISTART+MF(I)
        ENDDO
C
C
      CALL MIGUESS2(TTOT,T,M,N,MF,NF,NFRAGS)
      END IF
      RETURN
123      CALL DAREAD(IDAF,IODA,TTOT,M*M,316,0)
      IMO = 0
C      IOFF = -MF(1)
      IOFF = 0
         DO I=1,NFRAGS
            DO II=1,NF(I)
               IMO = IMO + 1
               DO K=1,M
                  T(K,IMO) = 0.0D+00
               END DO
               DO K=IOFF+1,IOFF+MF(I)
                  T(K,IMO) = TTOT(K,IOFF+II)
               END DO
            END DO
            IOFF = IOFF + MF(I)
         END DO
      CALL MIGUESS2(TTOT,T,M,N,MF,NF,NFRAGS)
      RETURN
      END
C
C*MODULE SCFMI   *DECK MIITER
      SUBROUTINE MIITER(T,SAO,HAO,P,SOCC,SWORK,TOTG,ROOT,ZMAT,U,
     *                  PREV,TTOT,ERN,DTOL,MMAS,ITERSCFMI,N,M,MTR,
     *                  IOPT,FFI,SSI,TVEC,MF,NF,NFRAGS,EIGEN,
     *                  EIGENWORK,INDX,TOTWORK,FDUM,FMO,SMO,EMO,
     *                  TMO,TNEW,DUM,WORK,SINV,ALPHA,DMAT,DS,VTTS,
     *                  TVDBVTTS,PFOCK,FP,FEFF,DBS,STVDOTVTTS,SDS,
     *                  SEFF,RIJ,ICONT,IORD,IIR,MAPP,TRAN,AOMULL,
     *                  NOROT,NITM,SOGLIA,IGUESS,SOCC1,ILELE3,
     *                  MAXDIIS,OLDFFI,OLDSSI,OLDERR,BMAT,DIISTOL,
     *                  DIIS,WORK2,QPIX,QPJX,SAOTR,GROSS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL PACK2E,OUT,DIRSCF,FDIFF,SCHWRZ
      LOGICAL DODIIS,DEM
      LOGICAL GOPARR,DSKWRK,MASWRK,TDSKWRK
      LOGICAL ABINIT,MINMEM,IZRFLS
      LOGICAL DOMBAS,SOSCF,FT15OP
      LOGICAL DIISON,DIIS
C
      DOUBLE PRECISION MOROKM
C
      DIMENSION DMAT(M,M)
      DIMENSION DS(M,M)
      DIMENSION VTTS(N,M)
      DIMENSION TVDBVTTS(M,M)
      DIMENSION PFOCK(M,M)
      DIMENSION FP(M,M)
      DIMENSION FEFF(M,M)
      DIMENSION DBS(M,M)
      DIMENSION STVDOTVTTS(M,M)
      DIMENSION SDS(M,M)
      DIMENSION SEFF(M,M)
      DIMENSION WORK2(*)
      DIMENSION FFI(*),SSI(*),TVEC(*)
      DIMENSION MF(NFRAGS),NF(NFRAGS)
      DIMENSION T(M,N),SAO(M,M),HAO(MTR),P(MTR),
     *          SOCC(N,N),SWORK(M,*),TOTG(M*M),ZMAT(MMAS,MMAS),
     *          ROOT(MMAS),U(MMAS,MMAS+2),PREV(MTR),TTOT(M,M)
C... GRADIENT'S VECTOR
      DIMENSION EIGEN(M),EIGENWORK(M),INDX(M),TOTWORK(M,M)
      DIMENSION FDUM(M,M),FMO(N,N)
      DIMENSION SMO(N,N)
      DIMENSION EMO(N),TMO(N,N+2)
      DIMENSION TNEW(M,M)
      DIMENSION DUM(MTR)
      DIMENSION WORK(N,N),SINV(N,N)
      DIMENSION AOMULL(M,M)
      DIMENSION SOCC1(*)
      DIMENSION OLDFFI(*),OLDSSI(*),OLDERR(*)
      DIMENSION BMAT(MAXDIIS+1,MAXDIIS+1)
C
C ***** FROM GAMESS
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXGSH=30, MXG2=MXGSH*MXGSH,
     *           MXATM=2000, MXFRG=1050, MXAO=8192,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /ACONV / RRSHFT,EXTTOL,DMPTOL,VSHTOL,IEXTIN
      COMMON /CONV  / DENTOL,EN,ETOT,EHF,EHF0,DIFF,ITER,ICALP,ICBET
      COMMON /DIISSO/ SOGTOL,ETHRSH,MAXDII,IRAF
      COMMON /DMPING/ SHIFTO,SHIFTV,DMPCUT,SWDIIS,DIRTHR
      COMMON /FMCOM / X(1)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FUNCT / ENERGY,EGRAD(3,MXATM)
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INDDIP/ TMUX,TMUY,TMUZ,TMUXD,TMUYD,TMUZD,MINMEM
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IJKAO,IJKMO,IDAF,NAV,IODA(950)
      COMMON /N2ELCT/ N2EL
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /TIMING/ CPU,WALL
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
C
      PARAMETER (ZERO=0.0D+00, HALF=0.5D+00, TWO=2.0D+00)
C
      SAVE FT15OP
C
      DATA FT15OP/.FALSE./
C
      DATA DEBUG/8HDEBUG   /, DBUGME/8HRHFCL   /
      DATA RHF/8HRHF     /
      DATA MOROKM/8HMOROKUMA/, OPTMZE/8HOPTIMIZE/
C
      DATA ZEROSCFMI/0.0D+00/,HALFSCFMI/0.5D+00/
      DATA TWOSCFMI/2.0D+00/
      DATA NONE/4HNONE/
C
C     ----- CLOSED SHELL HF-SCF CALCULATION -----
C     C.C.J.ROOTHAAN, REV.MOD.PHYS.  23,69-89(1951)
C
      ABINIT = MPCTYP.EQ.NONE
C
C         DIFFERENTIAL FOCK FORMATION IS NOT SUPPORTED
C
      IF(DIRSCF) FDIFF=.FALSE.
C
      IZRFLS = IZRF.NE.0
C
      OUT    = .FALSE.
C     PRDENS = .FALSE.
      IF(MASWRK) THEN
        OUT    = NPRINT.EQ.5  .OR. EXETYP.EQ.DEBUG .OR. EXETYP.EQ.DBUGME
      END IF
C
C     ----- START THE CLOCKS -----
C
      TIM0 = ZERO
      CALL TSECND(TIM0)
C
C     ----- SELECT CONVERGENCE METHOD(S) -----
C
      DODIIS = MOD(MCONV,32) .GE.16
      DEM    = MOD(MCONV,64) .GE.32
      SOSCF  = MOD(MCONV,128).GE.64
C
C     ----- SET POINTERS FOR PARTITIONING OF MEMORY -----
C
      IF (RUNTYP.EQ.MOROKM) THEN
         CALL STRHF(NUM,NA,NB,NE,NST,LL1,LL2,LL3,DOMBAS)
      ELSE
         NST = 1
         LL1 = NUM
         LL2 = (LL1*LL1+LL1)/2
         LL3 = LL1*LL1
         DOMBAS=.FALSE.
      END IF
C
      L0 = NQMT
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
C
      MAXIT2=(MAXIT*MAXIT+MAXIT)/2
      MAXIT3=MAXIT*MAXIT
C
      NSH2 = (NSHELL*NSHELL+NSHELL)/2
C
      CALL VALFM(LOADFM)
      CALL GOTFM(NGOTMX)
      IF(IOPT.GT.1) THEN
         WRITE(IW,*)'BEFORE ALLOCATION'
         WRITE(IW,*)'LOADFM',LOADFM
         WRITE(IW,*)'NGOTMX',NGOTMX
      END IF
C
      LSCR  = LOADFM + 1
      LIWRK = LSCR   + L1*8
      LWRK  = LIWRK  + L1
      LWRK1 = LWRK   + L1
      LWRK2 = LWRK1  + LL3
      LWRK3 = LWRK2  + LL3
      LEIG  = LWRK3  + LL3
      LOCC  = LEIG   + LL1
      LFAO  = LOCC   + L1
      LD    = LFAO   + LL2
      LVEC  = LD     + LL2
      LAST  = LVEC   + LL3
C
      CALL BASCHK(LMAX)
                    NANGM =  4
      IF(LMAX.EQ.2) NANGM =  6
      IF(LMAX.EQ.3) NANGM = 10
      IF(LMAX.EQ.4) NANGM = 15
      IF(LMAX.EQ.5) NANGM = 21
      IF(LMAX.EQ.6) NANGM = 28
      MAXG = NANGM**4
C
      IF(ABINIT) THEN
         IF(DIRSCF) THEN
            LBUF  = LAST
            LIBUF = LAST
            LXINTS= LAST
            LGHOND= LXINTS + NSH2
            LDSH  = LGHOND + MAXG
            LDDIJ = LDSH   + NSH2
            LAST  = LDDIJ  + 49*MXG2
         ELSE
            LXINTS= LAST
            LGHOND= LAST
            LDSH  = LAST
            LDDIJ = LAST
            LBUF  = LAST
            LIBUF = LBUF   + NINTMX
            LAST  = LIBUF  + NINTMX
         END IF
         LTWOEI= LAST
         LPTOT2= LAST
      ELSE
         LXINTS= LAST
         LGHOND= LAST
         LDSH  = LAST
         LDDIJ = LAST
         LBUF  = LAST
         LIBUF = LAST
         LTWOEI= LAST
         LPTOT2= LTWOEI + N2EL
         LAST  = LPTOT2 + NAT*16
      END IF
C
      IF(DODIIS) THEN
         LADIIS = LAST
         LBDIIS = LADIIS + MAXIT3
         LPDIIS = LBDIIS + MAXIT2
         LXDIIS = LPDIIS + MAXIT
         LIODII = LXDIIS + MAXIT
         LAST   = LIODII + 2*MAXDII
      ELSE
         LADIIS = LAST
         LBDIIS = LAST
         LPDIIS = LAST
         LXDIIS = LAST
         LIODII = LAST
      END IF
C
      IF(DEM) THEN
         LWRK4 = LAST
         LWRK5 = LWRK4 + L3
         LAST  = LWRK5 + L3
      ELSE
         LWRK4 = LAST
         LWRK5 = LAST
      END IF
C
      IF (IZRFLS .OR. IEFP.EQ.1) THEN
         LXDINT = LAST
         LYDINT = LXDINT  + L2
         LZDINT = LYDINT  + L2
         LAST   = LZDINT  + L2
      ELSE
         LXDINT = LAST
         LYDINT = LAST
         LZDINT = LAST
      END IF
      IF((IZRFLS .OR. IEFP.EQ.1) .AND. (DIRSCF.AND.FDIFF)) THEN
         NDIPER = LAST
         LDIPER = NDIPER  + L2
         LAST   = LDIPER  + L2
      ELSE
         NDIPER = LAST
         LDIPER = LAST
      END IF
      IF(IEFP.EQ.1)THEN
         LEFLD  = LAST
         MADD   = LEFLD  + 3*NPTTPT
         LABFLD = MADD   + 3*NPTTPT
         LAST   = LABFLD + 3*NPTTPT
      ELSE
         LEFLD  = LAST
         MADD   = LAST
         LABFLD = LAST
      END IF
C
      NFT15=15
      NPR=(NUM-NA)*NA
      IF (SOSCF) THEN
         LGRAD  = LAST
         LHESS  = LGRAD  + NPR
         LPGRAD = LHESS  + NPR
         LDISPL = LPGRAD + NPR
         LAST   = LDISPL + NPR
      ELSE
         LGRAD  = LAST
         LHESS  = LAST
         LPGRAD = LAST
         LDISPL = LAST
      END IF
C
      NEED = LAST - LOADFM - 1
      MINMEM = .FALSE.
      LH1   = LAST
      LQ    = LH1    + LL2
      LS    = LQ     + LL3
      LAST  = LS     + LL2
      NEED  = LAST -LOADFM -1
C
C     ----- GET MEMORY -----
C
      CALL GETFM(NEED)
      CALL VALFM(LOADFM)
      CALL GOTFM(NGOTMX)
C
C     ----- INITIALIZE OCCUPATION NUMBERS -----
C
      CALL VCLR(X(LOCC),1,L1)
      DO 120 I = NST,NA+NST-1
         X(LOCC-1+I) = TWO
  120 CONTINUE
C
C     ----- INITIALIZE VARIABLES -----
C
      IF(MAXIT.LE.0) MAXIT = 30
C
      IF(DIRSCF) THEN
         DENTOL = TWO * CONVHF
         ENGTOL = 1.0D-08
C        DIITOL = 1.0D-06
      ELSE
         DENTOL = CONVHF
         ENGTOL = 1.0D-09
C        DIITOL = 1.0D-07
      END IF
      IF(RUNTYP.EQ.OPTMZE  .AND.  NEVALS.GT.0) THEN
         GRMS = DDOT(3*NAT,EGRAD,1,EGRAD,1)
         GRMS = SQRT(GRMS/(3*NAT))
                               CFACT =  1.0D+00
         IF(GRMS.GT.0.005D+00) CFACT =  5.0D+00
         IF(GRMS.GT.0.020D+00) CFACT = 20.0D+00
         IF(GRMS.GT.0.100D+00) CFACT = 50.0D+00
         DENTOL = CFACT*DENTOL
         ENGTOL = CFACT*ENGTOL
C        DIITOL = CFACT*ENGTOL
      END IF
C
      EHF = ZERO
      ETOT= ZERO
      ICALP = 0
C     DAMP = ZERO
      IF (DMPCUT .LE. ZERO) DMPCUT = ZERO
      RRSHFT = ZERO
      DIFF = ZERO
C     IF(DAMPH .OR. VSHIFT) DAMP = ONE
C
      II = 0
      DO 140 I = 1,L1
         IA(I) = II
         II = II + I
  140 CONTINUE
C
C     ----- PRINT SCF OPTIONS -----
C
C     IF(SOME) THEN
C        WRITE(IW,9010) EN,MAXIT,MPUNCH,EXTRAH,DAMPH,
C    *                  VSHIFT,RSTRCT,DODIIS,DEM,SOSCF,DENTOL
C        IF(VTSCAL) WRITE(IW,9015) VTSCAL,MAXVT,VTCONV
C        IF(SOSCF) WRITE(IW,9017) NPR,SOGTOL
C        WRITE(IW,9020) NEED
C     END IF
C CHANGE...
C     IF(EXETYP.EQ.CHECK) GO TO 500
C
C
C     ----- EXCHANGE INTEGRALS FOR DIRECT SCF THRESHOLD TESTS -----
C
      NINT=0
      NSCHWZ=0
      IF(DIRSCF.AND.ABINIT) THEN
         SCHWRZ = ISCHWZ.GT.0
         IF (MASWRK) WRITE(IW,9030) SCHWRZ,FDIFF
         IF(SCHWRZ) CALL EXCHNG(X(LXINTS),X(LGHOND),X(LDDIJ),
     *                          NSH2,MAXG,INTTYP)
      ELSE
         SCHWRZ=.FALSE.
      END IF
C
C     ----- READ IN MOPAC 2-ELECTRON INTEGRALS -----
C
      IF(MPCTYP.NE.NONE) CALL DAREAD(IDAF,IODA,X(LTWOEI),N2EL,53,0)
C
C     ----- READ BARE NUCLEUS HAMILTONIAN -----
C     ----- READ OVERLAP INTEGRALS -----
C     ----- READ INITIAL DENSITY MATRIX -----
C     ----- READ SYMMETRY ADAPTED Q MATRIX -----
C     L0 = NUMBER OF CANONICAL ORTHONORMAL VECTORS KEPT.
C
      IF(.NOT.MINMEM) THEN
         CALL DAREAD(IDAF,IODA,X(LH1),LL2,11,0)
         CALL DAREAD(IDAF,IODA,X(LS) ,LL2,12,0)
         CALL DAREAD(IDAF,IODA,X(LQ) ,LL3,45,0)
      END IF
      CALL DAREAD(IDAF,IODA,X(LD),LL2,16,0)
C
      IF(SOSCF) THEN
         CALL DAREAD(IDAF,IODA,X(LVEC),LL3,15,0)
         CALL DAREAD(IDAF,IODA,X(LEIG),LL1,17,0)
         TDSKWRK = DSKWRK
         DSKWRK = .FALSE.
         IF(.NOT.FT15OP) THEN
            FT15OP = .TRUE.
            CALL SEQOPN(NFT15,'WORK15','UNKNOWN',.FALSE.,'UNFORMATTED')
         ELSE
            CALL SEQREW(NFT15)
         END IF
         DSKWRK = TDSKWRK
      END IF
C
C     ----- PREVIOUS DENSITY IS ZERO, PREVIOUS FOCK MAT. IS HCORE -----
C
      IF(DIRSCF  .AND.  FDIFF) THEN
         IF(MINMEM) CALL DAREAD(IDAF,IODA,X(LH1),LL2,11,0)
         CALL DAWRIT(IDAF,IODA,X(LH1),LL2,14,0)
         CALL VCLR(X(LWRK1),1,L2)
      END IF
C
C     ----- IF SCRF SAVE PRISTINE ONE ELEC. HAM. -----
C     -----    AND READ IN DIPOLE INTEGRALS      -----
C
      IF(IZRFLS) THEN
         CALL DAWRIT(IDAF,IODA,X(LH1),LL2,87,0)
         CALL DAREAD(IDAF,IODA,X(LXDINT),LL2,95,0)
         CALL DAREAD(IDAF,IODA,X(LYDINT),LL2,96,0)
         CALL DAREAD(IDAF,IODA,X(LZDINT),LL2,97,0)
      END IF
      IF(IEFP.EQ.1) CALL DAWRIT(IDAF,IODA,X(LH1),LL2,87,0)
C
C
C     ----- PRINT ITERATION HEADER -----
C
C     IF (MASWRK) THEN
C     IF(DAMPH  .OR.  VSHIFT) THEN
C        IF(SOSCF) THEN
C           IF(.NOT.DIRSCF) WRITE(IW,9040)
C           IF(     DIRSCF) WRITE(IW,9050)
C        ELSE
C           IF(.NOT.DIRSCF) WRITE(IW,9045)
C           IF(     DIRSCF) WRITE(IW,9055)
C        END IF
C     ELSE
C        IF(SOSCF) THEN
C           IF(.NOT.DIRSCF) WRITE(IW,9060)
C           IF(     DIRSCF) WRITE(IW,9070)
C        ELSE
C           IF(.NOT.DIRSCF) WRITE(IW,9065)
C           IF(     DIRSCF) WRITE(IW,9075)
C        END IF
C     END IF
C     END IF
C
C  *****
C
C
C
C     SCF-MI GUESS ORBITALS
C
      CALL DAREAD(IDAF,IODA,TTOT,LL3,15,0)
C---  WRITE(IW,6300) 'X(LVEC) COEFFICIENTS'
C---  CALL PRSQL(TTOT,M,M,M)
C---  WRITE(IW,6300) 'INITIAL COEFFICIENTS'
C---  CALL PRSQL(T,N,M,M)
C
C COPY MATRIX H INTO HAO X(LH1)-->HAO
C     CALL DCOPY(MTR,X(LH1),1,HAO,1)
      CALL DCOPY(MTR,X(LH1),1,HAO,1)
C
C TRASFORMATION X(LS)-->SAO(M,M)
      CALL CPYTSQ(X(LS),SAO,M,1)
C COPY NUCLEAR ENERGY EN --> ERN
      ERN=EN
C
C     SCF-MI GUESS ORBITALS
C
            CALL DAREAD(IDAF,IODA,TTOT,LL3,15,0)
C---  WRITE(IW,6300) 'X(LVEC) COEFFICIENTS'
C---  CALL PRSQL(TTOT,M,M,M)
      CALL MIINIT(T,TTOT,RIJ,IORD,IIR,MAPP,TRAN,NOROT,QPIX,QPJX,SAOTR,N,
     +M,IGUESS,NF,MF,NFRAGS,NAT,GROSS,NITM,SOGLIA)
C
C      WRITE(IW,*)'T-MATRIX'
C      CALL PRSQL(T,N,M,M)
C      WRITE(IW,*)'TTOT-MATRIX'
C      CALL PRSQL(TTOT,M,M,M)
C
C...
C INIZIALIZE ENERGY VARIABLES
C..
      ENER=ZEROSCFMI
      ENEROLD=ZEROSCFMI
C
C...
C FLAG FOR CONVERGENCE OF SCF-MI
C...
      ICONVSCFMI=0
      CALL VCLR(PREV,1,MTR)
      DIISON=.FALSE.
C LEVEL SHIFT OPTION
C...
      IF (ALPHA.GT.ZEROSCFMI) THEN
        IF(MASWRK) WRITE(IW,*)
        IF(MASWRK) WRITE(IW,*)'LEVEL SHIFT IS ACTIVE'
        IF(MASWRK) WRITE(IW,*)
      END IF
      CALL TSECND(TIM)
      TTIME=TIM
      WWALL=WALL
      IF(MASWRK) WRITE(IW,6131)
6131  FORMAT(1X,'START SCF-MI CALCULATION...')
C
      IF(MASWRK) WRITE(IW,*)
      IF(MASWRK) WRITE(IW,*) ' ITER   SCF-MI ENERGY    ENERGY'
     +,' CHANGE     DENSITY CHANGE    CPU TIME  WALL CL'
C.....
C SCFMI ITERATIONS
C.....
      DO 3000 IT=1,ITERSCFMI
      CALL FLSHBF(IW)
C
C     IF(IGUESS.EQ.-99.OR.IGUESS.EQ.-999)THEN
C       MOUT=30
C       IF(MASWRK) REWIND MOUT
C       IF(MASWRK) WRITE(MOUT) TTOT
C       IF(MASWRK) REWIND MOUT
C     ENDIF
C
C....
C OVERLAP MATRIX BETWEEN OCCUPIED OR TOTAL MOS
C.....
C SWORK=(TTS)T
      KK=1
      II=1
      KKK=MF(1)
      III=NF(1)
      DO L=1,NFRAGS
         DO J=1,M
            DO I=II,III
               SWORK(J,I)=0.0D+00
               DO K=KK,KKK
                  SWORK(J,I)=SWORK(J,I)+T(K,I)*SAO(K,J)
               END DO
            END DO
         END DO
         KK=KKK+1
         II=III+1
         KKK=KKK+MF(L+1)
         III=III+NF(L+1)
       END DO
C
C SOCC1=(TTST)TRIANG
C
       KK=1
       JJ=1
       KKK=MF(1)
       JJJ=NF(1)
       ILELE=1
       DO L=1,NFRAGS
          DO J=JJ,JJJ
             DO 1111 I=1,J
                  SOCC1(ILELE)=0.0D+00
                  DO K=KK,KKK
                     SOCC1(ILELE)=SOCC1(ILELE)+SWORK(K,I)*T(K,J)
                  END DO
                  ILELE=ILELE+1
1111         CONTINUE
          END DO
          KK=KKK+1
          JJ=JJJ+1
          KKK=KKK+MF(L+1)
          JJJ=JJJ+NF(L+1)
        END DO
      CALL CPYTSQ(SOCC1,SOCC,N,1)
      CALL MIMXINV(SOCC,SWORK,N)
C CALCOLA T(TTST)-1 TRASPOSTO
C
      KK=1
      II=1
      KKK=NF(1)
      III=MF(1)
      DO L=1,NFRAGS
         DO J=1,N
            DO I=II,III
               SWORK(J,I)=0.0D+00
               DO K=KK,KKK
                  SWORK(J,I)=SWORK(J,I)+T(I,K)*SOCC(K,J)
               END DO
            END DO
         END DO
         KK=KKK+1
         II=III+1
         KKK=KKK+NF(L+1)
         III=III+MF(L+1)
      END DO
C...
C DENSITY MATRIX FOR SCFMI ALGORITHM P(MTR)
C DENSITY D IN VECTOR P(MTR)
C P-->X(LD)
C...
      IPUT=0
      DO IQ=1,M
         DO IPP=1,IQ
            IPUT=IPUT+1
            P(IPUT)=ZEROSCFMI
            X(LD+IPUT-1)=ZEROSCFMI
            DO K=1,N
               P(IPUT)=P(IPUT)+T(IQ,K)*SWORK(K,IPP)
            END DO
            X(LD+IPUT-1)=P(IPUT)*2.0D+00
            DUM(IPUT)=X(LD+IPUT-1)
         END DO
         P(IPUT)=P(IPUT)*HALFSCFMI
      END DO
      CALL CPYTSQ(P,DMAT,M,1)
      DO ININO=1,M
         DMAT(ININO,ININO)=DMAT(ININO,ININO)*TWOSCFMI
      END DO
      DDIFFSCFMI=ZEROSCFMI
      DO 4700 I=1,MTR
4700  DDIFFSCFMI=DDIFFSCFMI+(P(I)-PREV(I))**2
      DDIFFSCFMI=SQRT(DDIFFSCFMI/MTR)
C     WRITE(IW,*) 'DENSITY MATRIX CONVERGENCE:',DDIFFSCFMI
      CALL DCOPY(MTR,P,1,PREV,1)
C
C
C
C999  CALL GETIME('START FOCK OPERATOR CALCULATION')
C
C
C **** GAMESS
C GAMESS CALCULATION OF  FOCK MATRIX
C                   X(LFAO)
C
C           ********************
C     ***** START RHF ITERATIONS *****
C           ********************
C
C     DO 300 ITER=1,MAXIT
      CALL TSECND(TITER0)
C
C     ----- CONSTRUCT SKELETON 2 ELECTRON FOCK MATRIX -----
C     FOR DIRECT SCF,
C     EITHER -LWRK1- IS THE CHANGE IN THE DENSITY, WITH THE CHANGE
C     TO THE FOCK MATRIX TO BE COMPUTED IN -LFAO-, OR ELSE
C     -LWRK1- IS THE CURRENT DENSITY, WITH ALL 2 E- TERMS IN -LFAO-.
C
      IF (ABINIT) THEN
         IF(DIRSCF) THEN
            IF(FDIFF) THEN
               CALL VSUB(X(LWRK1),1,X(LD),1,X(LWRK1),1,L2)
            ELSE
               CALL DCOPY(L2,X(LD),1,X(LWRK1),1)
            END IF
            DUMMY = 0.0D+00
            IF(SCHWRZ) THEN
               DUMMY = 0.0D+00
               CALL SHLDEN(RHF,X(LWRK1),DUMMY,DUMMY,X(LDSH),
     *                     IA,L1,L2,NSH2,1)
            END IF
C
            CALL VCLR(X(LFAO),1,L2)
            CALL TWOEI(RHF,DIRSCF,.FALSE.,.FALSE.,.FALSE.,
     *                 INTTYP,SCHWRZ,NINT,NSCHWZ,L1,L2,
     *                 DUMMY,DUMMY,IDUMMY,1,
     *                 X(LXINTS),NSH2,X(LGHOND),MAXG,X(LDDIJ),
     *                 IA,X(LWRK1),X(LFAO),DUMMY,DUMMY,X(LDSH),
     *                 DUMMY,DUMMY,1)
C
C
C   --- DIRSCF OFF DIAGONAL ELEMENTS ARE DOUBLE THE CORRECT VALUE ---
C
            CALL DSCAL(L2,HALF,X(LFAO),1)
            II=LFAO-1
            DO 210 I=1,L1
               II = II+I
               X(II) = X(II) + X(II)
  210       CONTINUE
         ELSE
            TDSKWRK = DSKWRK
            DSKWRK  = .TRUE.
            CALL SEQREW(IJKAO)
            CALL HSTAR(X(LD),X(LFAO),X(LBUF),X(LIBUF),
     *                 NINTMX,IA,NOPK,.FALSE.)
C
            DSKWRK  = TDSKWRK
         END IF
C
C   --- SUM UP PARTIAL FOCK MATRICES ---
C
         IF(GOPARR) THEN
            CALL DDI_GSUMF(1000,X(LFAO),L2)
            CALL DDI_GSUMI(1001,NINT   ,1)
            CALL DDI_GSUMI(1002,NSCHWZ ,1)
         END IF
      ELSE
         HFCO = -HALF
         CALL MPCG(X(LFAO),X(LFAO),X(LD),DUMMY,DUMMY,X(LPTOT2),
     *             X(LTWOEI),L2,HFCO)
      END IF
C
C
      IF(OUT) THEN
         WRITE(IW,*) '*** START OF DEBUG OUTPUT FOR ITERATION',ITER
         WRITE(IW,*) 'SKELETON FOCK MATRIX'
         CALL PRTRIL(X(LFAO),L1)
      END IF
C
C     ----- SYMMETRIZE SKELETON FOCK MATRIX -----
C     SCRATCH AREA AT LWRK1 IS USED AS A SYMMETRIC MATRIX
C
      CALL SYMH(X(LFAO),X(LWRK1),IA)
      IF(OUT) THEN
         WRITE(IW,*) 'SYMMETRIZED FOCK OPERATOR'
         CALL PRTRIL(X(LFAO),L1)
      END IF
C
C    -- IF SCRF GET PRISTINE ONE ELEC. HAM., CALC. INDUCED DIPOLES, --
C                   ----- ADD PERTURBATION. -----
C
      IF(IZRFLS) THEN
         CALL DAREAD(IDAF,IODA,X(LH1),LL2,87,0)
         IF(DIRSCF .AND. FDIFF) THEN
            CALL RCNFLD(X(LH1),X(LFAO),DUMMY,X(LD),X(NDIPER),
     *                  X(LDIPER),X(LXDINT),X(LYDINT),X(LZDINT),L2,L1,0)
         ELSE
            CALL RCNFLD(X(LH1),DUMMY,DUMMY,X(LD),DUMMY,
     *                  DUMMY,X(LXDINT),X(LYDINT),X(LZDINT),L2,L1,0)
         END IF
      END IF
C
C    -- IF FRAG GET PRISTINE ONE ELEC. HAM., CALC. INDUCED DIPOLES, --
C                   ----- ADD PERTURBATION. -----
C
      IF(IEFP.EQ.1) THEN
         CALL DAREAD(IDAF,IODA,X(LH1),LL2,87,0)
         IF(DIRSCF .AND. FDIFF) THEN
            CALL EFPCMP(X(LH1),X(LFAO),DUMMY,X(LD),X(NDIPER),
     *                  X(LDIPER),X(LXDINT),X(LYDINT),X(LZDINT),
     *                  X(LEFLD),X(MADD),X(LABFLD),L2,L1,0)
         ELSE
            CALL EFPCMP(X(LH1),DUMMY,DUMMY,X(LD),DUMMY,DUMMY,
     *                  X(LXDINT),X(LYDINT),X(LZDINT),X(LEFLD),
     *                  X(MADD),X(LABFLD),L2,L1,0)
         END IF
      END IF
C
C     ----- COMPLETE THE FOCK MATRIX -----
C     WE EITHER ADD THE FOCK OPERATOR OF THE PREVIOUS ITERATION TO
C     THE CHANGE IN THE FOCK MATRIX, OR WE ADD THE ONE ELECTRON
C     HAMILTONIAN TO THE FULLY COMPUTED TWO ELECTRON OPERATOR.
C
      IF(ABINIT  .AND.  DIRSCF  .AND.  FDIFF) THEN
         CALL DAREAD(IDAF,IODA,X(LWRK1),LL2,14,0)
         CALL VADD(X(LFAO),1,X(LWRK1),1,X(LFAO),1,L2)
         IF(MINMEM) CALL DAREAD(IDAF,IODA,X(LH1),LL2,11,0)
      ELSE
         IF(DOMBAS) CALL CTMOTR(X(LWRK1),X(LFAO),X(LWRK2),X(LSCR),LL1)
         IF(MINMEM) CALL DAREAD(IDAF,IODA,X(LH1),LL2,11,0)
         CALL VADD(X(LFAO),1,X(LH1),1,X(LFAO),1,L2)
      END IF
      CALL DAWRIT(IDAF,IODA,X(LFAO),LL2,14,0)
C
C
C
C     ----- DEBUG PRINT OF ENERGY, FOCK MATRIX, AND LAGRANGIAN -----
C     -FV- AND -EIJ- AT LWRK1 AND LWRK2 ARE SQUARE MATRICES
C
      IF(OUT) THEN
C        WRITE(IW,9080) EHF1,EHF2,EHF
         WRITE(IW,*) 'TOTAL FOCK OPERATOR'
         CALL PRTRIL(X(LFAO),L1)
      END IF
C     NPRINT=5
      IF(NPRINT.EQ.5 .OR. EXETYP.EQ.DEBUG .OR. EXETYP.EQ.DBUGME)
     *    CALL DAREAD(IDAF,IODA,X(LVEC),LL3,15,0)
      IF(OUT) THEN
         CALL MAKFV(X(LFAO),X(LVEC),X(LWRK1),X(LWRK),IA,1,NA,L1,L1)
         CALL MAKEIJ(X(LVEC),X(LWRK1),X(LWRK2),X(LWRK),L0,L0,L1,L0,L1)
         WRITE(IW,*) 'LAGRANGE MULTIPLIER MATRIX'
         CALL PRSQ(X(LWRK2),NA,L0,L1)
      END IF
C
C  ****
C
C COPY X(LFAO) --> TOTG
      CALL DCOPY(MTR,X(LFAO),1,TOTG,1)
C
C
C     CALL GETIME('END OF FOCK OPERATOR CALCULATION ')
C...
C COMPUTE THE ENERGY
C...
      ENER=ZEROSCFMI
      DO 5100 I=1,MTR
5100  ENER=ENER+(TOTG(I)+HAO(I))*P(I)
      ENER=ERN+ENER+ENER
C
      EDIFFOUT=ENER-ENEROLD
      ENEROLD=ENER
C
      CALL TSECND(TIM)
      TTIME1=-TTIME+TIM
      TTIME=TIM
      WWALL1=-WWALL+WALL
      WWALL=WALL
      IF(MASWRK) WRITE(IW,6789) IT,ENER,EDIFFOUT,DDIFFSCFMI,
     *                          TTIME1,WWALL1
C
 6789 FORMAT(1X,I3,2X,F18.10,2X,F16.10,2X,F15.10,3X,F7.1,3X,F7.1)
 6790 FORMAT(1X,'FINAL SCF-MI ENERGY=',3X,F20.10)
C
      IF(DIIS)THEN
        IF (DIISON) THEN
        ELSE
        IF((DDIFFSCFMI.LE.DIISTOL).AND.(ITERSCFMI.NE.1))THEN
           DIISON=.TRUE.
           NOOLD=0
           IF (NOOLD.EQ.0) THEN
           IF(MASWRK) WRITE(IW,*)'START DIIS PROCEDURE'
           END IF
        ENDIF
        END IF
      ENDIF
      IF((DDIFFSCFMI.LE.DTOL).AND.(ITERSCFMI.NE.1)) THEN
        IF(MASWRK) WRITE(IW,*)
        IF(MASWRK) WRITE(IW,*)
        IF(MASWRK) WRITE(IW,6790) ENER
        IF(MASWRK) WRITE(IW,*)
        ETOT=ENER
C       IF(IOPT.NE.0) CALL GETIME('END OF THE ITERATIONS')
        ICONVSCFMI=1
      ENDIF
      IF(ICONVSCFMI.EQ.1) THEN
          GO TO 3999
      END IF
C.....
C BUILD AND SOLVE GENERALIZED EIGENVALUE PROBLEMS
C.....
C     IF(LMAP.EQ.M) CALL REMAP(TOTG,U,MAP,LMAP,INDXMAP,MTR,1)
      CALL VRSQUARE(TOTG,TOTG,M)
C     CALL PRINTSCFMI(TOTG,M,M,0)
C
C.....
C BUILD PFOCK MATRIX PFOCK= I - DS + TV DB (VT+S)
C.....
C...  DXS
C     WRITE(IW,*)
C     WRITE(IW,6300) 'DXS'
      CALL MXATB(DMAT,SAO,M,M,M,DS)
C.... CALCOLO DB(VTTS)
      CALL VCLR(VTTS,1,N*M)
      II=1
      JJ=1
      III=NF(1)
      JJJ=MF(1)
      DO L=1,NFRAGS
         DO J=JJ,JJJ
            DO I=II,III
               DO K=1,M
                  VTTS(I,J)=VTTS(I,J)+SWORK(I,K)*SAO(K,J)
               END DO
            END DO
         END DO
         II=III+1
         JJ=JJJ+1
         III=III+NF(L+1)
         JJJ=JJJ+MF(L+1)
       END DO
C... CALCOLO (TV)  X  DB{V X [(T+) X S]}
C
      KK=1
      JJ=1
      KKK=NF(1)
      JJJ=MF(1)
      DO L=1,NFRAGS
         DO J=JJ,JJJ
            DO I=1,M
               TVDBVTTS(I,J)=0.0D+00
               DO K=KK,KKK
                  TVDBVTTS(I,J)=TVDBVTTS(I,J)+SWORK(K,I)*
     +                                        VTTS(K,J)
               END DO
            END DO
         END DO
         KK=KKK+1
         JJ=JJJ+1
         KKK=KKK+NF(L+1)
         JJJ=JJJ+MF(L+1)
       END DO
C...  PFOCK= I - DS + TV DB (VT+S)
      CALL VCLR(PFOCK,1,M*M)
      CALL VSUB(DS,1,TVDBVTTS,1,PFOCK,1,M*M)
      DO ININO=1,M
         PFOCK(ININO,ININO)=PFOCK(ININO,ININO)+1.0D+00
      END DO
C CALCOLO PTF TRASPOSTA (IN FP)
      CALL DCOPY(M*M,TOTG,1,TVDBVTTS,1)
      DO I=1,M
         DO J=1,M
            FP(J,I)=0.0D+00
            DO K=1,M
               FP(J,I)=FP(J,I)+PFOCK(K,I)*TVDBVTTS(K,J)
            END DO
         END DO
      END DO
C CALCOLO DB(PTFP)
      JJ=1
      JJJ=MF(1)
      DO L=1,NFRAGS
         DO J=JJ,JJJ
            DO I=JJ,JJJ
               FEFF(I,J)=0.0D+00
               DO K=1,M
                  FEFF(I,J)=FEFF(I,J)+FP(K,I)*
     +                                    PFOCK(K,J)
               END DO
            END DO
         END DO
         JJ=JJJ+1
         JJJ=JJJ+MF(L+1)
       END DO
C     S' = DB{S} + DB{STV}XDB{V(T+)S} - DB{SDS}
C...
C     DB{S}
C REMOVE ANY OFF DIAGONAL BLOCK COMPONENTS IN SAO
      CALL DCOPY(M*M,SAO,1,DBS,1)
      CALL DBMAKE(DBS,M,M,MF,MF,NFRAGS)
C... CALCOLO DB((VTTS)T) X DB(VTTS)
C
      KK=1
      II=1
      KKK=NF(1)
      III=MF(1)
      DO L=1,NFRAGS
         DO J=II,III
            DO I=II,III
               STVDOTVTTS(I,J)=0.0D+00
               DO K=KK,KKK
                  STVDOTVTTS(I,J)=STVDOTVTTS(I,J)+VTTS(K,I)*
     +                                            VTTS(K,J)
               END DO
            END DO
         END DO
         KK=KKK+1
         II=III+1
         KKK=KKK+NF(L+1)
         III=III+MF(L+1)
      END DO
C  CALCOLO DB(SDS)
      II=1
      III=MF(1)
      DO L=1,NFRAGS
         DO J=II,III
            DO I=II,III
               SDS(I,J)=0.0D+00
               DO K=1,M
                  SDS(I,J)=SDS(I,J)+SAO(I,K)*DS(K,J)
               END DO
            END DO
         END DO
         II=III+1
         III=III+MF(L+1)
      END DO
C     S' = DB{S} + DB{STV}XDB{V(T+)S} - DB{SDS}
      CALL VCLR(SEFF,1,M*M)
      CALL VADD(DBS,1,STVDOTVTTS,1,SEFF,1,M*M)
      CALL VSUB(SDS,1,SEFF,1,SEFF,1,M*M)
      IPUNT=0
      IOFFSET=1
      DO IFI=1,NFRAGS
         CALL SUBCOPY(FFI(IOFFSET),FEFF,MF(IFI),M,IPUNT)
         CALL SUBCOPY(SSI(IOFFSET),SEFF,MF(IFI),M,IPUNT)
         CALL SUBCOPY(TVEC(IOFFSET),TTOT,MF(IFI),M,IPUNT)
         IPUNT=IPUNT+MF(IFI)
         IOFFSET=IOFFSET+MF(IFI)*MF(IFI)
      END DO
C
      AMAXERR=0.0D+00
C
      IF (DIISON)THEN
C
         NOOLD=NOOLD+1
         IF(NOOLD.GT.MAXDIIS)THEN
              NOOLD=MAXDIIS
              DO I=1,(MAXDIIS-1)*ILELE3
                   OLDFFI(I)=OLDFFI(I+ILELE3)
                   OLDSSI(I)=OLDSSI(I+ILELE3)
                   OLDERR(I)=OLDERR(I+ILELE3)
              END DO
         END IF
C
         DO I=1,ILELE3
            J=I+(NOOLD-1)*ILELE3
            OLDFFI(J)=FFI(I)
            OLDSSI(J)=SSI(I)
         END DO
      IOFFSET=1
      IFROM=1+ILELE3*(NOOLD-1)
      DO IFI=1,NFRAGS
         CALL ERRORE(TVEC(IOFFSET),SSI(IOFFSET),FFI(IOFFSET),
     +          AOMULL,TOTWORK,OLDERR(IFROM),MF(IFI),NF(IFI),
     +           AMAXERR)
C        IPUNT=IPUNT+MF(IFI)
         IOFFSET=IOFFSET+MF(IFI)*MF(IFI)
         IFROM=IFROM+MF(IFI)*MF(IFI)
      END DO
C
4321  IF(NOOLD.GE.2)THEN
         BMAT(1,1)=0.0D+00
         DO I=2,NOOLD+1
            BMAT(1,I)=1.0D+00
            BMAT(I,1)=1.0D+00
            DO J=2,NOOLD+1
               BMAT(I,J)=0.0D+00
            ENDDO
         ENDDO
C
C
      IPUNT=1
      DO IFI=1,NFRAGS
          DO I=2,NOOLD+1
            IFROM1=IPUNT+(I-2)*ILELE3
            DO J=2,NOOLD+1
               IFROM2=IPUNT+(J-2)*ILELE3
               CALL CREABMAT(BMAT(I,J),OLDERR(IFROM1),
     +              OLDERR(IFROM2),MF(IFI),NF(IFI))
            END DO
          END DO
          IPUNT=IPUNT+MF(IFI)*MF(IFI)
      END DO
C
        CALL SUBCOPY(AOMULL,BMAT,NOOLD+1,MAXDIIS+1,0)
        CALL DGEFA(AOMULL,NOOLD+1,NOOLD+1,SWORK,INFO)
C
        CALL SUBCOPY(AOMULL,BMAT,NOOLD+1,MAXDIIS+1,0)
        IF ((INFO.NE.0))THEN
           DO I=1,(NOOLD-1)*ILELE3
              OLDFFI(I)=OLDFFI(I+ILELE3)
              OLDSSI(I)=OLDSSI(I+ILELE3)
           END DO
           NOOLD=NOOLD-1
           GO TO 4321
        END IF
        CALL MIMXINV(AOMULL,SWORK,NOOLD+1)
        DO I=1,ILELE3
           FFI(I)=0.0D+00
           SSI(I)=0.0D+00
         END DO
C
       IPUNT=1
        DO IFI=1,NFRAGS
         DO II=1,NOOLD
          CO=AOMULL(II+1,1)
          IFROM=IPUNT+(II-1)*ILELE3
C
          CALL CREAFFSS(FFI(IPUNT),SSI(IPUNT),OLDFFI(IFROM),
     +     OLDSSI(IFROM),CO,MF(IFI))
        END DO
C
         IPUNT=IPUNT+MF(IFI)*MF(IFI)
       END DO
       END IF
       END IF
C
      IPUNT=0
      IOFFSET=1
      DO IFI=1,NFRAGS
         CALL MIGSOLVE(FFI(IOFFSET),SSI(IOFFSET),TVEC(IOFFSET),ROOT,
     +          SWORK,ZMAT,U,MF(IFI),TTOT,EIGEN,M,IPUNT,IOPT,ALPHA,
     +          DIISON,WORK2,NF(IFI))
         IPUNT=IPUNT+MF(IFI)
         IOFFSET=IOFFSET+MF(IFI)*MF(IFI)
      END DO
C.....
C BUILD OCCUPIED MOS MATRIX T
C.....
      KCOL1=0
      KCOL2=0
      DO IFI=1,NFRAGS
         DO J=1,NF(IFI)
            DO I=1,M
               T(I,KCOL1+J)=TTOT(I,KCOL2+J)
            END DO
         END DO
         KCOL1=KCOL1+NF(IFI)
         KCOL2=KCOL2+MF(IFI)
       END DO
      IF(IOPT.GT.1) THEN
        IF(MASWRK) WRITE(IW,6300) 'UPDATED COEFFICIENTS'
        IF(MASWRK) CALL PRSQL(T,N,M,M)
      ENDIF
C....
      CALL DAWRIT(IDAF,IODA,TTOT,M*M,316,0)
C
C     UPDATED COEFFICIENTS IN CHANNEL 2
C      IWRITE=2
C      REWIND(IWRITE)
C      WRITE(IWRITE)TTOT
C      CALL FLSHBF(IWRITE)
C
3000  CONTINUE
      IF(MASWRK) WRITE(IW,*)
      IF(MASWRK) WRITE(IW,*) 'EXCEEDED MAXIMUM NUMBER OF ITERATIONS '
      IF(MASWRK) WRITE(IW,*)
C
C
      IF(IOPT.GT.0)THEN
      IF(MASWRK) CALL PRINTSCFMI(T,M,N,0)
      ENDIF
3999  CONTINUE
C     WRITE(IW,6300) 'FINAL COEFFICIENTS'
      IF(IOPT.GT.0)THEN
        IF(MASWRK) WRITE(IW,6300) 'CONVERGED ORBITALS'
        IF(MASWRK) WRITE(IW,*)
        IF(MASWRK) CALL PRSQL(T,N,M,M)
      ENDIF
C
C
C
C       ----- PUNCH THE SCFMI NON ORTHOGONAL ORBITALS -----
C                       FILE *.DAT
CC      IF(MPUNCH.GT.0 .AND. MASWRK) THEN
CC         CALL TMDATE(TIMSTR)
C          WRITE(IP,8000) TIMSTR,TITLE,ETOT,EN,ITER
CC          LPUN = NA
CC         IF (MPUNCH .EQ. 2) LPUN = L1
CC          WRITE(IP, FMT='(5H $VEC)')
C         CALL PUSQL(X(LVEC),LPUN,L1,L1)
C          CALL PUSQL(TTOT,LPUN,L1,L1)
CC          CALL PUSQL(TTOT,M,M,M)
C          CALL PUSQL(T,N,M,M)
CC          WRITE(IP, FMT='(5H $END)')
CC      END IF
C     OPEN(FILE='MOSCFMI',UNIT=2,FORM='UNFORMATTED',
C    +     ACCESS='SEQUENTIAL')
C     REWIND(2)
C     WRITE(UNIT=2) TTOT
C     CLOSE(UNIT=2)
C     WRITE(IW,*)' TERMINATA SCRITTURA '
C     STOP
C...
C  START OF SOME TRASFORMATIONS FOR GRADIENT...
C...
C
C...
C BUILD INDX
C...
      CALL MIINDX (M,EIGEN,INDX)
C...
C ORDERING OF EIGENVALUES AND EIGENVECTORS
C...
      CALL ORDMATVECT(TTOT,TOTWORK,EIGEN,EIGENWORK,INDX,M)
C
C...
C FROM LAST FOCK (FAO) BUILD FMO
C...
      CALL CPYTSQ(TOTG,FDUM,M,1)
      ICONT = 0
      DO I=1,N
         DO J=1,N
            ICONT = ICONT + 1
            FMO(I,J)=ZEROSCFMI
            DO L=1,M
               DO IK=1,M
                   FMO(I,J)=FMO(I,J)+TTOT(L,I)*TTOT(IK,J)*FDUM(L,IK)
               END DO
            END DO
          END DO
       END DO
C
C
C...
C TRANSORMATION SAO -> SMO
C...
      DO I=1,N
         DO J=1,N
            SMO(I,J)=ZEROSCFMI
            DO L=1,M
               DO IK=1,M
                   SMO(I,J)=SMO(I,J)+TTOT(L,I)*TTOT(IK,J)*SAO(L,IK)
               END DO
            END DO
          END DO
       END DO
C
C
C **NOTE**: SCFMI ORBITALS ARE NOW ORDERED  WITH RESPECT
C           TO THEIR EIGENVALUES
C .....
C...
C DIAGONALIZATION OF FMO TMO =  SMO TMO EMO
C...
      CALL MISOLVE(FMO,SMO,EMO,TMO,IOPT,WORK,SINV,N,WORK2)
C
       IF(IOPT.GT.1)THEN
      IF(MASWRK) WRITE(IW,*)'AUTOVALORI FMO'
      DO I=1,N
         IF(MASWRK) WRITE(IW,*)I,EMO(I)
      END DO
      IF(MASWRK) WRITE(IW,*)'AUTOVETTORI FMO'
      IF(MASWRK) CALL PRINTSCFMI(TMO,N,N,0)
      ENDIF
C
C...
C SAVE FINAL EIGENVALUES INTO EIGEN(I)
C...
      DO II=1,N
          EIGEN(II)=EMO(II)
      END DO
C
C
C...
C BUILD NEW MOS (ORTHONORMAL IN THE OCCUPIED PART)
C...
      CALL NEWORB (TTOT,TNEW,TMO,M,N)
C
C     ----- SAVE MO-S + DENSITY + ORBITAL ENERGIES -----
C
      NDAF = 15
C     CALL SCFSAV(X(LVEC),X(LD),X(LEIG),NDAF,LL1,LL2,LL3)
      CALL SCFSAV(TNEW,DUM,EIGEN,NDAF,LL1,LL2,LL3)
      IF (IOPT.GT.1)THEN
        IF(MASWRK) WRITE(IW,*)'NUOVI MOS'
        IF(MASWRK) CALL PRINTSCFMI(TNEW,M,M,0)
      END IF
C...
C CLEAN THE MEMORY USED IN ITERATE
C..
      CALL RETFM(NEED)
C
      CALL VALFM(LOADFM)
      CALL GOTFM(NGOTMX)
      RETURN
C
 6300 FORMAT(/,1X,'========>',1X,A)
 9030 FORMAT(/1X,'DIRECT SCF CALCULATION, SCHWRZ=',L1,'   FDIFF=',L1)
      END
C
C*MODULE SCFMI   *DECK MIGSOLVE
      SUBROUTINE MIGSOLVE(GMATA,SMAT,TI,ROOT,WORK,ZMAT,U,NDIM,
     + TTOT,EIGEN,M,IPUNT,IOPT,ALPHA,DIISON,WORK2,NFDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DIISON
C
      DIMENSION TTOT(M,M),WORK(M,M),GMATA(NDIM,NDIM),EIGEN(M),
     *          SMAT(NDIM,NDIM),ZMAT(NDIM,NDIM),TI(NDIM,NDIM),
     *          ROOT(NDIM),U(NDIM,NDIM+2),WORK2(*)
C
C     IF(IOPT.GT.1) CALL GETIME('SOLVE SECULAR PROBLEM')
C     CALL VCLR(RMAT,1,NDIM*M)
C
C...
C LEVEL SHIFT
C...
      TRESH=1.0D-14
*     WRITE(6,*)'ALPHA-->', ALPHA
      IF (.NOT.DIISON)THEN
        IF (ALPHA.GT.TRESH) THEN
         CALL VCLR(WORK,1,M*M)
         CALL LEVELSHIFT(GMATA,SMAT,TI,NDIM,ALPHA,
     +           WORK(2*NDIM+1,2*NDIM+1),NFDIM)
         CALL VCLR(WORK,1,M*M)
        END IF
      END IF
C.....
C  SOLVE EACH SECULAR PROBLEM
C ...
C
C       CALL PRINTSCFMI(GMATA,NDIM,NDIM,0)
        CALL VBSOLVE(GMATA,SMAT,ROOT,WORK,ZMAT,U,NDIM,IOPT,
     +  WORK2)
        DO J=1,NDIM
           DO I=1,NDIM
              TTOT(IPUNT+I,IPUNT+J)=U(I,J)
           END DO
        END DO
C       CALL PRINTSCFMI(U,NDIM,NDIM,0)
C
        DO I=1,NDIM
           EIGEN(IPUNT+I)=ROOT(I)
C          WRITE(6,*)IPUNT+I,EIGEN(IPUNT+I)
        END DO
C
      RETURN
      END
C
C*MODULE SCFMI   *DECK VBSOLVE
      SUBROUTINE VBSOLVE(G,S,ROOT,WORK,SINV,U,NTSP,IOPT,WORK2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION ROOT(NTSP),WORK(NTSP*NTSP),WORK2(*),
     *          SINV(NTSP,NTSP),G(NTSP,NTSP),S(NTSP,NTSP),
     *          U(NTSP,NTSP+2)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA ZERO/0.0D+00/,ONE/1.0D+00/
C
      CALL MXEIG(S,ROOT,U,WORK,NTSP,WORK2)
      IF(IOPT.GT.1) WRITE(6,6100) NTSP
      DO 15 I=1,NTSP
15    ROOT(I)=ONE/SQRT(ROOT(I))
C
C GENERATE S(-1/2)
C
      DO 2200 L=1,NTSP
      DO 2200 K=1,NTSP
      SINV(K,L)=ZERO
      DO 2200 IRUN=1,NTSP
2200  SINV(K,L)=SINV(K,L) + U(K,IRUN)*U(L,IRUN)*ROOT(IRUN)
C
C GENERATE S(-1/2) * H * S(-1/2)
      CALL MXATB(SINV,G,NTSP,NTSP,NTSP,U)
      CALL MXATB(U,SINV,NTSP,NTSP,NTSP,G)
      CALL MXEIG(G,ROOT,U,WORK,NTSP,WORK2)
      CALL DCOPY(NTSP*NTSP,U,1,G,1)
      CALL MXATB(SINV,G,NTSP,NTSP,NTSP,U)
      IF(IOPT.GT.1)THEN
        IF(MASWRK)  WRITE(6,'(4G20.12)') ROOT
      ENDIF
      RETURN
6100  FORMAT(1X,'ORDER OF THE MATRIX',I5,/,' EIGENVALUES OF H-E*S')
      END
C
C*MODULE SCFMI   *DECK MIINDX
      SUBROUTINE MIINDX(N,ARRIN,INDX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION ARRIN(N),INDX(N)
C
      DO  J=1,N
         INDX(J)=J
      END DO
      IF(N.EQ.1) RETURN
      L=N/2+1
      IR=N
 10   CONTINUE
          IF (L.GT.1) THEN
              L=L-1
              INDXT=INDX(L)
              Q=ARRIN(INDXT)
          ELSE
              INDXT=INDX(IR)
              Q=ARRIN(INDXT)
              INDX(IR)=INDX(1)
              IR=IR-1
              IF (IR.EQ.1) THEN
                 INDX(1)=INDXT
                 RETURN
              END IF
          END IF
          I=L
          J=L+L
 20       IF (J.LE.IR) THEN
              IF(J.LT.IR) THEN
                 IF (ARRIN(INDX(J)).LT.ARRIN(INDX(J+1))) J=J+1
              END IF
              IF (Q.LT.ARRIN(INDX(J))) THEN
                  INDX(I)=INDX(J)
                  I=J
                  J=J+J
              ELSE
                  J=IR+1
              END IF
          GO TO 20
          END IF
          INDX(I)=INDXT
      GO TO 10
      END
C
C*MODULE SCFMI   *DECK ORDMATVECT
      SUBROUTINE ORDMATVECT(TTOT,TWORK,EIGEN,EIGENWORK,INDX,M)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION TTOT(M,M),TWORK(M,M),EIGEN(M),EIGENWORK(M),INDX(M)
C
C ORDERING OF EIGENVECTORS WITH RESPECT TO THE EIGENVALUES
C
      DO I=1,M
         EIGENWORK(I)=EIGEN(I)
      END DO
C
      DO I=1,M
         JINDX=INDX(I)
         EIGEN(I)= EIGENWORK(JINDX)
      END DO
C
      DO J=1,M
         DO I=1,M
            TWORK(I,J)=TTOT(I,J)
         END DO
      END DO
C
      DO J=1,M
         JINDX=INDX(J)
         DO I=1,M
            TTOT(I,J)=TWORK(I,JINDX)
         END DO
      END DO
C
      RETURN
      END
C
C*MODULE SCFMI   *DECK MISOLVE
      SUBROUTINE MISOLVE(G,S,ROOT,U,IOPT,WORK,SINV,NTSP,WORK2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION ROOT(NTSP),WORK(NTSP*NTSP),WORK2(*),
     *          SINV(NTSP,NTSP),G(NTSP,NTSP),S(NTSP,NTSP),
     *          U(NTSP,NTSP+2)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA ZERO/0.0D+00/,ONE/1.0D+00/
C
      CALL MXEIG(S,ROOT,U,WORK,NTSP,WORK2)
      DO 15 I=1,NTSP
15    ROOT(I)=ONE/SQRT(ROOT(I))
C
C GENERATE S(-1/2)
      DO 2200 L=1,NTSP
      DO 2200 K=1,NTSP
      SINV(K,L)=ZERO
      DO 2200 IRUN=1,NTSP
2200  SINV(K,L)=SINV(K,L) + U(K,IRUN)*U(L,IRUN)*ROOT(IRUN)
C
C GENERATE S(-1/2) * H * S(-1/2)
      CALL MXATB(SINV,G,NTSP,NTSP,NTSP,U)
      CALL MXATB(U,SINV,NTSP,NTSP,NTSP,G)
      CALL MXEIG(G,ROOT,U,WORK,NTSP,WORK2)
      CALL DCOPY(NTSP*NTSP,U,1,G,1)
      CALL MXATB(SINV,G,NTSP,NTSP,NTSP,U)
C
      IF(IOPT.GT.1)THEN
        IF(MASWRK)  WRITE(6,'(4G20.12)') ROOT
      ENDIF
      RETURN
      END
C
C*MODULE SCFMI   *DECK NEWORB
      SUBROUTINE NEWORB (T,TNEW,AUTOV,M,N)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION T(M,M)
      DIMENSION TNEW(M,M)
      DIMENSION AUTOV(N,N)
C
      ZERO=0.0D+00
C
      DO J=1,N
         DO I=1,M
            TNEW(I,J)=ZERO
         END DO
      END DO
C
      DO J=N+1,M
         DO I=1,M
            TNEW(I,J)=T(I,J)
         END DO
      END DO
      DO IK=1,N
          DO L=1,N
             DO I=1,M
                 TNEW(I,IK)=TNEW(I,IK)+AUTOV(L,IK)*T(I,L)
             END DO
           END DO
       END DO
C
      RETURN
      END
C
C*MODULE SCFMI   *DECK MXEIG
      SUBROUTINE MXEIG(FULL, VAL, VEX, WORK, NDIM, VFULL)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DIISON,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MAXNFRAGS=10)
C
      DIMENSION VFULL(NDIM*(NDIM+1)/2),FULL(NDIM*NDIM),VAL(NDIM),
     *          VEX(NDIM,NDIM),WORK(NDIM*NDIM)
C
      COMMON /IOFILE/ IR,IW,IP,IJKAO,IJKMO,IDAF,NAV,IODA(950)
      COMMON /MIOPT / DTOL,CNVLOC,DIISTL,ALPHA,
     *                NF(MAXNFRAGS),MF(MAXNFRAGS),NFRAGS,MSHIFT,ITER,
     *                ITLOC,IGUESS,MXDIIS,
     *                IOPT,MIFLAG,DIISON
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     MATRIX FULL (NDIMXNDIM) BECOMES TRIANGULAR IN
C     VFULL(NDIM*(NDIM+1)/2)
C
      CALL CPYSQT(FULL,VFULL,NDIM,1)
C
C     FULL EIGENVALUES (VAL) AND EIGENVECTORS (VEX)
C
      CALL GLDIAG(NDIM,NDIM,NDIM,VFULL,WORK,VAL,VEX,
     *            IGERR,WORK(8*NDIM+1))
      IF(IGERR.NE.0) THEN
        IF(MASWRK) WRITE(IW,*) ' NO CONVERGENCE OF GLDIAG IN MXEIG'
        CALL ABRT
      ENDIF
C
      IF(IOPT.GT.1) THEN
        IF(MASWRK) WRITE(6,*) ' MATRIX DIAGONALIZED (GLDIAG)'
      ENDIF
C
      RETURN
      END
C
C*MODULE SCFMI   *DECK MIMXINV
      SUBROUTINE MIMXINV(FULL, FWORK, NDIM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,DIISON
C
      PARAMETER (MAXNFRAGS=10)
C
      DIMENSION FULL(NDIM,NDIM), FWORK(NDIM,*), DET(2)
C
      COMMON /IOFILE/ IR,IW,IP,IJKAO,IJKMO,IDAF,NAV,IODA(950)
      COMMON /MIOPT / DTOL,CNVLOC,DIISTL,ALPHA,
     *                NF(MAXNFRAGS),MF(MAXNFRAGS),NFRAGS,MSHIFT,ITER,
     *                ITLOC,IGUESS,MXDIIS,
     *                IOPT,MIFLAG,DIISON
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C          DRIVER FOR MATRIX INVERSION VIA LINPACK ROUTINES
C
      CALL DGEFA(FULL,NDIM,NDIM,FWORK,INFO)
      IF(INFO.EQ.0) THEN
         JOB=1
         CALL DGEDI(FULL,NDIM,NDIM,FWORK,DET,FWORK(1,2),JOB)
C
      ELSE
        IF(MASWRK)  WRITE(IW,*) ' SINGULAR MATRIX IN MIMXINV'
         CALL ABRT
      ENDIF
      IF(IOPT.GT.1) THEN
        IF(MASWRK) WRITE(6,*) ' MATRIX INVERTED (DGEDI)'
      ENDIF
      RETURN
      END
C
C*MODULE SCFMI   *DECK PRINTSCFMI
      SUBROUTINE PRINTSCFMI(RMAT,NROW,NCOL,ITYPE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (NBUF=8)
      LOGICAL GOPARR,DSKWRK,MASWRK
      DIMENSION BUFFER(NBUF),IBUF(NBUF),RMAT(*)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C PRINTS MATRIX RMAT, STORED ACCORDING TO ITYPE
C
      JIN=1
1100  JEND=JIN+NBUF-1
      IF(NCOL.LE.NBUF) JEND=NCOL
      IF(JEND.GT.NCOL+NBUF-1) RETURN
      JEND=MIN0(NCOL,JEND)
      K=0
      DO 2100 J=JIN,JEND
      K=K+1
2100  IBUF(K)=J
      IF(MASWRK) WRITE(6,6100) (IBUF(I),I=1,JEND-JIN+1)
      DO 3100 I=1,NROW
      K=0
      DO 3200 J=JIN,JEND
      K=K+1
      IF(ITYPE.EQ.0) THEN
        IND=(J-1)*NROW+I
      ELSEIF(ITYPE.EQ.1) THEN
        IF(I.GE.J) THEN
          IND=I*(I-1)/2 + J
        ELSE
          IND=J*(J-1)/2 + I
        ENDIF
      ELSE
        IND=(I-1)*NCOL+J
      ENDIF
3200  BUFFER(K)=RMAT(IND)
3100  IF(MASWRK) WRITE(6,6200) I,(BUFFER(II),II=1,JEND-JIN+1)
      JIN=JEND+1
      IF(NCOL.GT.NBUF) GOTO 1100
6100  FORMAT(/,4X,8(3X,I3,3X))
6200  FORMAT(1X,I4,8F9.5)
      END
C
C*MODULE SCFMI   *DECK MXATB
      SUBROUTINE MXATB(A,B,N1,N2,N3,C)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(N1,N2), B(N2,N3), C(N1,N3)
      DATA ZERO/0.0D+00/
C
      DO 1100 J=1,N3
      DO 1100 I=1,N1
      C(I,J)=ZERO
      DO 1100 IRUN=1,N2
1100  C(I,J)=C(I,J) + A(I,IRUN)*B(IRUN,J)
      RETURN
      END
C
C*MODULE SCFMI   *DECK DBMAKE
      SUBROUTINE DBMAKE(A,N1,N2,MVECT,NVECT,NBLOCK)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(N1,N2),MVECT(NBLOCK),NVECT(NBLOCK)
C...
C TRASFORMA LA MATRICE A IN DB{A}
C...             A -> DB{A}
C
      MTOT=0
      DO I=1,NBLOCK
         MTOT=MTOT+MVECT(I)
      END DO
C
      KCOL=0
      KROW=0
      DO IBLOCK=1,NBLOCK
         DO J=KCOL+1,KCOL+NVECT(IBLOCK)
            DO I=1,KROW
               A(I,J)=0.0D+00
            END DO
            DO I=KROW+MVECT(IBLOCK)+1,MTOT
               A(I,J)=0.0D+00
            END DO
         END DO
         KCOL=KCOL+NVECT(IBLOCK)
         KROW=KROW+MVECT(IBLOCK)
      END DO
      RETURN
      END
C
C*MODULE SCFMI   *DECK SUBCOPY
      SUBROUTINE SUBCOPY(A,ATOT,MF,M,IPUNT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(MF,MF),ATOT(M,M)
C
      DO J=1,MF
         DO I=1,MF
            A(I,J)=ATOT(IPUNT+I,IPUNT+J)
         END DO
      END DO
      RETURN
      END
C
C*MODULE SCFMI   *DECK LEVELSHIFT
      SUBROUTINE LEVELSHIFT(GMATA,SMAT,TI,NDIM,ALPHA,ST,NFDIM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION GMATA(NDIM,*),SMAT(NDIM,*),TI(NDIM,*),ST(NDIM,*)
C...
C LEVEL SHIFT
C...
C...  S X T
      CALL MXATB(SMAT,TI,NDIM,NDIM,NFDIM,ST)
C
C...  F=F + ALPHA * [S - (S T * TT S)]
      DO I=1,NDIM
        DO J=1,NDIM
          GMATA(I,J)=GMATA(I,J)/ALPHA+SMAT(I,J)
          DO K=1,NFDIM
            GMATA(I,J)=GMATA(I,J)-ST(I,K)*ST(J,K)
          ENDDO
          GMATA(I,J)=ALPHA*GMATA(I,J)
        ENDDO
      ENDDO
C
      RETURN
      END
C
C*MODULE SCFMI   *DECK ERRORE
      SUBROUTINE ERRORE(T,S,F,WORK,WORK1,Z,M,N,AMAXERR)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION T(M,M),S(M,M),F(M,M),WORK(M,M),WORK1(M,M),Z(M,M)
C
C WORK = TTT
C
      DO I=1,M
        DO J=1,M
           WORK(I,J)=0.0D+00
           DO K=1,N
              WORK(I,J)=WORK(I,J)+T(I,K)*T(J,K)
           END DO
         END DO
       END DO
C
C WORK1= STTT
C
      DO I=1,M
        DO J=1,M
           WORK1(I,J)=0.0D+00
           DO K=1,M
              WORK1(I,J)=WORK1(I,J)-S(I,K)*WORK(K,J)
           END DO
         END DO
       END DO
C
       DO I=1,M
          WORK1(I,I)=WORK1(I,I)+1.0D+00
       END DO
C
C WORK = (I - STTT)F
C
        DO I=1,M
        DO J=1,M
           WORK(I,J)=0.0D+00
           DO K=1,M
              WORK(I,J)=WORK(I,J)+WORK1(I,K)*F(K,J)
           END DO
         END DO
       END DO
C
C WORK1 = (I -STTT)FT
C
      DO I=1,M
        DO J=1,N
           Z(I,J)=0.0D+00
           DO K=1,M
              Z(I,J)=Z(I,J)+WORK(I,K)*T(K,J)
           END DO
           IF(Z(I,J).GT.AMAXERR)THEN
              AMAXERR=Z(I,J)
           ENDIF
         END DO
       END DO
       RETURN
       END
C
C*MODULE SCFMI   *DECK CREABMAT
      SUBROUTINE CREABMAT(B,Z1,Z2,M,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Z1(M,M),Z2(M,M)
C
      DO K=1,M
         DO J=1,N
           B=B+Z1(K,J)*Z2(K,J)
         END DO
      END DO
      RETURN
      END
C
C*MODULE SCFMI   *DECK CREAFFSS
      SUBROUTINE CREAFFSS(F,S,FOLD,SOLD,C,M)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION F(M,M),S(M,M),FOLD(M,M),SOLD(M,M)
C
      DO I=1,M
         DO J=1,M
            F(I,J)=F(I,J)+C*FOLD(I,J)
            S(I,J)=S(I,J)+C*SOLD(I,J)
         END DO
      END DO
      RETURN
      END
C
C*MODULE SCFMI   *DECK POPLOC2
      SUBROUTINE POPLOC2(IPASS,VECIN,TRAN,VECOUT,IORD,IIR,
     *                  MAP,RIJ,QPIX,QPJX,NOROT,SAO,
     *                  L1,L2,N1,N2,NATOMS,NOUT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION LMORND
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXAO=8192)
C
      DIMENSION VECIN(L1,N1),TRAN(N1,N1),SAO(L2),NOROT(N1,N1),
     *          VECOUT(L1,N1),RIJ(N2,NATOMS),
     *          IIR(L1),MAP(L1),IORD(L1),QPIX(L1),QPJX(L1)
C
      COMMON /IOFILE/ IR,IW,IP,IJKAO,IJKMO,IDAF,NAV,IODA(950)
      COMMON /OPTLOC/ CVGLOC,MAXLOC,IPRTLO,ISYMLO,IFCORE,NOUTA,NOUTB,
     *                MOOUTA(MXAO),MOOUTB(MXAO)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00, FOUR=4.0D+00,
     *           TENM3=1.0D-03, TENM8=1.0D-08, TENM10=1.0D-10)
C
      INDX(I,J) = ((MAX(I,J)*(MAX(I,J)-1))/2+MIN(I,J))
C
C     POPLOC CALCULATES THE LOCALIZED MOLECULAR ORBITALS BY THE
C     POPULATION METHOD OF J. PIPEK AND P. G. MEZEY,
C     J. CHEM. PHYS. 90, 4916 (1989),
C     THIS IMPLEMENTATION IS DERIVED FROM QCPE PROGRAM 354,
C     AND WAS WRITTEN BY JOHN MONTGOMERY, JULY 1991
C
      NORB  = N1
      NBASIS= L1
      MAXIT = MAXLOC
C
      CALL AOLIM
C     DO 20 K=1,NBASIS
C     DO 20 I=1,NATOMS
C     IF(LIMLOW(I).LE.K.AND.LIMSUP(I).GE.K) MAP(K)=I
C  20 CONTINUE
C
      NREDO = 0
  110 CONTINUE
      NREDO=NREDO+1
C
C          CONSTRUCT INITIAL ATOMIC POPULATIONS
C
      CALL VCLR(RIJ,1,N2*NATOMS)
      CALL DAREAD(IDAF,IODA,SAO,L2,12,0)
      IJ = 0
      DO 280 I = 1,NORB
         DO 280 J = 1,I
            IJ = IJ+1
            DO 260 K = 1,NBASIS
               KK = MAP(K)
               DO 260 L = 1,NBASIS
                  LL = MAP(L)
                  KL = INDX(K,L)
                  SUM = VECIN(K,I)*VECIN(L,J)*SAO(KL)/TWO
                  RIJ(IJ,KK) = RIJ(IJ,KK) + SUM
                  RIJ(IJ,LL) = RIJ(IJ,LL) + SUM
  260       CONTINUE
  280 CONTINUE
C
C          COMPUTE INITIAL LOCALIZATION SUM
C
      SUMRR = ZERO
      DO 320 I=1,NORB
         IF(NOROT(I,I).EQ.1) GO TO 320
         II = INDX(I,I)
         DO 310 K=1,NATOMS
            SUMRR = SUMRR + RIJ(II,K)**2
  310    CONTINUE
  320 CONTINUE
      IF (MASWRK) WRITE(IW,9010) SUMRR
      IF(IPRTLO.EQ.1) CALL POPQAT(RIJ,VECOUT,N2,NORB,NATOMS)
C
C          SEED THE RANDOM FUNCTION, INITIALIZE
C          THE LOCALIZATION TRANSFORMATION, ETC.
C
      ONEPT0 = ONE
      IF (NREDO.EQ.1) XX = LMORND(ONEPT0,VECIN,L1)
C
      CALL VCLR(TRAN,1,NORB*NORB)
      DO 340 I = 1,NORB
         TRAN(I,I) = ONE
  340 CONTINUE
      ITER = 0
      SHIFT = ATAN(ONEPT0)
C
C          BEGIN LOCALIZATION CYCLES
C
  360 CONTINUE
      CHANGE = ZERO
      ITER = ITER+1
      DO 380 I = 1,NORB
         IIR(I) = I
  380 CONTINUE
      NNN = NORB
      DO 400 I = 1,NORB
         XX = LMORND(CHANGE,VECIN,L1)
         III = INT(XX*NNN+ONE)
         IORD(I) = IIR(III)
         IIR(III) = IIR(NNN)
         NNN = NNN-1
  400 CONTINUE
C
C        FOR EACH PAIR OF ORBITALS A TWO DIMENSIONAL UNITARY
C        TRANSFORMATION IS PERFORMED. THE TRANSFORMATION IS
C
C           PSI'(I) =  COS(T)*PSI(I) + SIN(T)*PSI(J)  AND
C           PSI'(J) = -SIN(T)*PSI(I) + COS(T)*PSI(J).
C
C        LOCALIZATION REQUIRES THAT T BE SUCH AS TO MAXIMIZE
C        THE SUM OF THE SQUARES OF THE ATOMIC POPULATIONS.
C
      DO 920 III = 1,NORB
      I = IORD(III)
      IF(NOROT(I,I).EQ.1) GO TO 920
      II = INDX(I,I)
      JM = 1
      RM = ZERO
      TM = ZERO
      SM = ZERO
      CM = ONE
      DO 580 J = 1,NORB
      IF(I.EQ.J) GO TO 580
      IF(NOROT(I,J).EQ.1) GO TO 580
      IJ = INDX(I,J)
      JJ = INDX(J,J)
      T = ZERO
      TX = ZERO
      DO 480 KK = 1,NATOMS
         T = T + FOUR*RIJ(IJ,KK)**2 - RIJ(II,KK)**2 - RIJ(JJ,KK)**2
     *         + TWO*RIJ(II,KK)*RIJ(JJ,KK)
         TX = TX + RIJ(IJ,KK)*(RIJ(JJ,KK) - RIJ(II,KK))
  480 CONTINUE
      IF ((ABS(T) .LE. TENM10) .AND. (ABS(TX) .LE. TENM10)) GO TO 580
      TX = FOUR*TX
      T = ATAN2(TX,T)/FOUR
      SIGN = ONE
      IF (T .GT. ZERO) SIGN = -ONE
      T = T+SIGN*SHIFT
      ITIM = 0
  500 ITIM = ITIM+1
      S = SIN(T)
      C = COS(T)
      RIN = ZERO
      DO 520 KK = 1,NATOMS
         QPI = C*C*RIJ(II,KK)+S*S*RIJ(JJ,KK)+TWO*C*S*RIJ(IJ,KK)
         QPJ = C*C*RIJ(JJ,KK)+S*S*RIJ(II,KK)-TWO*C*S*RIJ(IJ,KK)
         RIN = RIN+QPI*QPI+QPJ*QPJ-RIJ(II,KK)**2-RIJ(JJ,KK)**2
  520 CONTINUE
      TTEST = ABS(T)-SHIFT
      IF ((ABS(T) .LE. TENM8) .OR. (ABS(TTEST) .LE. TENM8)) GO TO 560
      IF (RIN .GE. -TENM8) GO TO 560
      IF (ITIM .LE. 1) GO TO 540
      IF (MASWRK) THEN
         WRITE (IW,9020) I,J
         WRITE (IW,9030) T,S,C,RIN
      ENDIF
      RETURN
C
  540 SIGN = ONE
      IF (T .GT. ZERO) SIGN = -ONE
      T = T+SHIFT*SIGN
      GO TO 500
C
  560 IF (RIN .LE. RM) GO TO 580
      RM = RIN
      TM = T
      SM = S
      CM = C
      JM = J
  580 CONTINUE
C
      RIN = RM
      T = TM
      S = SM
      C = CM
      J = JM
      IJ = INDX(I,J)
      JJ = INDX(J,J)
      IF(NOROT(I,J).EQ.1) GO TO 920
C
C        ACCUMULATE THE 2X2 ROTATION
C
      CHANGE = CHANGE+T*T
      CALL DROT(NORB,TRAN(1,I),1,TRAN(1,J),1,C,S)
C
C        UPDATE THE ATOMIC POPULATIONS
C
      DO 880 KK = 1,NATOMS
         QPI = C*C*RIJ(II,KK)+S*S*RIJ(JJ,KK)+TWO*C*S*RIJ(IJ,KK)
         QPJ = C*C*RIJ(JJ,KK)+S*S*RIJ(II,KK)-TWO*C*S*RIJ(IJ,KK)
         QPIJ = (C*C-S*S)*RIJ(IJ,KK)+C*S*(RIJ(JJ,KK)-RIJ(II,KK))
         DO 720 K = 1,NORB
            IF (I.EQ.K.OR.J.EQ.K) GOTO 720
            IK = INDX(I,K)
            JK = INDX(J,K)
            QPIX(K) = C*RIJ(IK,KK)+S*RIJ(JK,KK)
            QPJX(K) = C*RIJ(JK,KK)-S*RIJ(IK,KK)
            RIJ(IK,KK) = QPIX(K)
            RIJ(JK,KK) = QPJX(K)
  720    CONTINUE
         RIN = RIN+QPI+QPJ-RIJ(II,KK)-RIJ(JJ,KK)
         RIJ(II,KK) = QPI
         RIJ(JJ,KK) = QPJ
         RIJ(IJ,KK) = QPIJ
  880 CONTINUE
  920 CONTINUE
C
C          TEST FOR CONVERGENCE OF LOCALIZATION PROCEDURE
C
      CHANGE = SQRT(TWO*CHANGE/(NORB*(NORB-1)))
      IF((NPRINT.EQ.-7) .AND. MASWRK) WRITE(IW,9050) ITER,CHANGE
      IF((MOD(ITER,10).EQ.0) .AND. MASWRK) WRITE(IW,9050) ITER,CHANGE
      IF(ITER.LT.MAXIT  .AND.  CHANGE.GT.TENM3*CVGLOC) GO TO 360
C     IF(ITER.LT.MAXIT  .AND.  CHANGE.GT.1.0D-14) GO TO 360
      IF(CHANGE.LE.CVGLOC) GO TO 1000
         IF((NREDO.LE.2) .AND. MASWRK) WRITE(IW,9060)
         IF(NREDO.LE.2) GO TO 110
            IF (MASWRK) WRITE(IW,9070)
            RETURN
C
C          FINISHED WITH LOCALIZATION CYCLES
C
 1000 CONTINUE
      IF (MASWRK) WRITE (IW,9080) ITER
C
C        TRANSFORM TO FINAL ORBITALS, COPY VIRTUAL SPACE
C
      CALL MRARBR(VECIN,NBASIS,NBASIS,NORB,TRAN,NORB,NORB,VECOUT,NBASIS)
      NVIRT = NBASIS - NORB
      CALL DCOPY(NVIRT*NBASIS,VECIN(1,NORB+1),1,VECOUT(1,NORB+1),1)
C
C          PRINT/PUNCH FINAL LOCALIZED ORBITALS
C
C     IF (MASWRK) WRITE(IW,9090)
C9090 FORMAT(/10X,'THE PIPEK-MEZEY POPULATION LOCALIZED ORBITALS ARE')
C     CALL PRSQL(VECOUT,NORB,L1,L1)
C
C     ----- SAVE THIS FOR EVENTUAL MP2 LOCALIZATION -----
C     WE NEED THE MATRIX WITHOUT THE CORE
C
      DO 1110 I = 1,NORB
         DO 1100 J = 1,NORB
            TRAN(I,J) = TRAN(I+NOUT,J+NOUT)
 1100    CONTINUE
 1110 CONTINUE
      NLOC = NORB - NOUT
      CALL SQZMAT(TRAN,NORB,NLOC,NORB)
      IF(IPASS.EQ.1) CALL DAWRIT(IDAF,IODA,TRAN,NLOC*NLOC,73,0)
      IF(IPASS.EQ.2) CALL DAWRIT(IDAF,IODA,TRAN,NLOC*NLOC,74,0)
C
CC      IF (MASWRK) WRITE (IW,9130) NORB,NFROZ,NBASIS
C
C      IF((IPASS.EQ.1) .AND. MASWRK) WRITE(IP,8000)
C      CALL PUSQL(VECOUT,NORB,L1,L1)
C      IF((SCFTYP.NE.UHF  .OR.  IPASS.EQ.2) .AND. MASWRK)
C     *   WRITE(IP,8010)
C8000 FORMAT('PIPEK-MEZEY POPULATION LOCALIZED ORBITALS'/' $VEC')
C8010 FORMAT(' $END')
C
C          CONSTRUCT FINAL ATOMIC POPULATIONS
C
      CALL VCLR(RIJ,1,N2*NATOMS)
      CALL DAREAD(IDAF,IODA,SAO,L2,12,0)
      IJ = 0
      DO 1280 I = 1,NORB
         DO 1280 J = 1,I
            IJ = IJ+1
            DO 1260 K = 1,NBASIS
               KK = MAP(K)
               DO 1260 L = 1,NBASIS
                  LL = MAP(L)
                  KL = INDX(K,L)
                  SUM = VECOUT(K,I)*VECOUT(L,J)*SAO(KL)/TWO
                  RIJ(IJ,KK) = RIJ(IJ,KK) + SUM
                  RIJ(IJ,LL) = RIJ(IJ,LL) + SUM
 1260       CONTINUE
 1280 CONTINUE
C
C               COMPUTE FINAL LOCALIZATION SUM (SKIP FROZEN MO-S)
C
      SUMRR = ZERO
      DO 1380 I=1,NORB
         IF(NOROT(I,I).EQ.1) GO TO 1380
         II = INDX(I,I)
         DO 1370 K = 1,NATOMS
            SUMRR = SUMRR + RIJ(II,K)**2
 1370    CONTINUE
 1380 CONTINUE
      IF (MASWRK) WRITE(IW,9140) SUMRR
C      CALL POPQAT(RIJ,VECIN,N2,NORB,NATOMS)
      RETURN
C
 9010 FORMAT(/10X,'THE INITIAL LOCALIZATION SUM IS',F14.6)
 9020 FORMAT(1H0,'NO ROTATION INCREASES ATOMIC POPULATIONS',
     *           ' --- LOCALIZATION ABORTED'/
     *           10X,2HI=,I3,5X,2HJ=,I3)
 9030 FORMAT(5X,8HTHETA = ,G20.10/5X,12HSIN(THETA)= ,F10.7,
     +     15H   COS(THETA)= ,F10.7/5X,29HTOTAL CHANGE TO THIS POINT = ,
     +     G20.10)
 9050 FORMAT(10X,'ITERATION',I4,'   ORBITAL CHANGE=',G20.10)
 9060 FORMAT (/10X,
     +     114HLOCALIZATION HAS BEEN UNSUCESSFUL. PROGRAM WILL RESTART W
     +ITH NEW RANDOM NUMBER AND ROTATION SEQUENCE FOR ORBITALS.)
 9070 FORMAT(1H0,'LOCALIZATION FAILS --- LOCALIZATION ABORTED')
 9080 FORMAT(10X,'LOCALIZATION CONVERGED IN',I4,' ITERATIONS')
 9140 FORMAT(10X,'THE FINAL LOCALIZATION SUM IS',F16.6/)
      END
C
C*MODULE SCFMI   *DECK MIGUESS2
      SUBROUTINE MIGUESS2(TTOT,T,M,N,MF,NF,NFRAGS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION TTOT(M*M),T(M*N),MF(NFRAGS),NF(NFRAGS)
C
      CALL VCLR(TTOT,1,M*M)
      DO I=1,(M*NF(1))
         TTOT(I)=T(I)
      END DO
      ISTART=M*NF(1)+1
      ISHIFT=0
      DO I=2,NFRAGS
         ISHIFT=ISHIFT+M*(MF(I-1)-NF(I-1))
         DO II=ISTART,(ISTART+(M*NF(I)-1))
              III=II+ISHIFT
              TTOT(III)=T(II)
              ISAVE=II
         END DO
         ISTART=ISAVE+1
      END DO
      RETURN
      END
C
C*MODULE SCFMI   *DECK VRSQUARE
      SUBROUTINE VRSQUARE(STR,S,N)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION S(N,N),STR(*)
C
C         NOTE THAT -STR- AND -S- MAY BE THE SAME ADDRESS
C
      IJ=N*(N+1)/2+1
      DO 5100 I=N,1,-1
      DO 5100 J=I,1,-1
      IJ=IJ-1
5100  S(J,I)=STR(IJ)
      DO 5200 I=1,N
      DO 5200 J=1,I-1
5200  S(I,J)=S(J,I)
      RETURN
      END
