      SUBROUTINE LOCAL(C,MDIM,NOCC,EIG)
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      dimension C(MDIM,MDIM), EIG(MAXORB)
      common /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
C**********************************************************************
C
C   LOCALISATION SUBROUTINE
C ON INPUT
C        C = EIGENVECTORS IN AN MDIM*MDIM MATRIX
C        NOCC = NUMBER OF FILLED LEVELS
C        NORBS = NUMBER OF ORBITALS
C        NUMAT = NUMBER OF ATOMS
C        NLAST   = INTEGER ARRAY OF ATOM ORBITAL COUNTERS
C        NFIRST   = INTEGER ARRAY OF ATOM ORBITAL COUNTERS
C
C       SUBROUTINE MAXIMIZES (PSI)**4
C       REFERENCE_
C       A NEW RAPID METHOD FOR ORBITAL LOCALISATION, P.G. PERKINS AND
C       J.J.P. STEWART, J.C.S. FARADAY (II) 77, 000, (1981).
C
C       MODIFIED AND CORRECTED TO AVOID SIGMA-PI ORBITAL MIXING BY
C       JUAN CARLOS PANIAGUA, UNIVERSITY OF BARCELONA, MAY 1983.
C
C**********************************************************************
      common /SCRACH/ COLD(MAXORB,MAXORB),XDUMY(MAXPAR**2-MAXORB*MAXORB)
      dimension EIG1(MAXORB),PSI1(MAXORB),PSI2(MAXORB),
     1          CII(MAXORB), REFEIG(MAXORB),IEL(20)
      character*2 ELEMNT(99)
      save ELEMNT
      data ELEMNT/'H','HE',
     1 'LI','BE','B','C','N','O','F','NE',
     2 'NA','MG','AL','SI','P','S','CL','AR',
     3 'K','CA','SC','TI','V','CR','MN','FE','CO','NI','CU',
     4 'ZN','GA','GE','AS','SE','BR','KR',
     5 'RB','SR','Y','ZR','NB','MO','TC','RU','RH','PD','AG',
     6 'CD','IN','SN','SB','TE','I','XE',
     7 'CS','BA','LA','CE','PR','ND','PM','SM','EU','GD','TB','DY',
     8 'HO','ER','TM','YB','LU','HF','TA','W','RE','OS','IR','PT',
     9 'AU','HG','TL','PB','BI','PO','AT','RN',
     1 'FR','RA','AC','TH','PA','U','NP','PU','AM','CM','BK','CF','XX'/
      NITER=100
      EPS=1.0D-7
      do 10 I=1,NORBS
         REFEIG(I)=EIG(I)
         do 10 J=1,NORBS
   10 COLD(I,J)=C(I,J)
      ITER=0
   20 continue
      SUM=0.D0
      ITER=ITER+1
      do 80 I=1,NOCC
         do 70 J=1,NOCC
            if (J.EQ.I) go to 70
            XIJJJ=0.0D0
            XJIII=0.0D0
            XIIII=0.0D0
            XJJJJ=0.0D0
            XIJIJ=0.0D0
            XIIJJ=0.0D0
            do 30 K=1,NORBS
               PSI1(K)=C(K,I)
   30       PSI2(K)=C(K,J)
C NOW FOLLOWS THE RATE-DETERMINING STEP FOR THE CALCULATION
            do 50 K1=1,NUMAT
               KL=NFIRST(K1)
               KU=NLAST(K1)
               DIJ=0.D0
               DII=0.D0
               DJJ=0.D0
C$DOIT ASIS
               do 40 K=KL,KU
                  DIJ=DIJ+PSI1(K)*PSI2(K)
                  DII=DII+PSI1(K)*PSI1(K)
                  DJJ=DJJ+PSI2(K)*PSI2(K)
   40          continue
               XIJJJ=XIJJJ+DIJ*DJJ
               XJIII=XJIII+DIJ*DII
               XIIII=XIIII+DII*DII
               XJJJJ=XJJJJ+DJJ*DJJ
               XIJIJ=XIJIJ+DIJ*DIJ
               XIIJJ=XIIJJ+DII*DJJ
   50       continue
            AIJ=XIJIJ-(XIIII+XJJJJ-2.0D0*XIIJJ)/4.0D0
            BIJ=XJIII-XIJJJ
            CA=SQRT(AIJ*AIJ+BIJ*BIJ)
            SA=AIJ+CA
            if (SA.LT.1.0D-14) go to 70
            SUM=SUM+SA
            CA=-AIJ/CA
            CA=(1.0D0+SQRT((1.0D0+CA)/2.0D0))/2.0D0
            if ((2.0D0*CA-1.0D0)*BIJ.LT.0.0D0) CA=1.0D0-CA
            SA=SQRT(1.0D0-CA)
            CA=SQRT(CA)
            do 60 K=1,NORBS
               C(K,I)=CA*PSI1(K)+SA*PSI2(K)
   60       C(K,J)=-SA*PSI1(K)+CA*PSI2(K)
   70    continue
   80 continue
      SUM1=0.D0
      do 100 I=1,NOCC
         do 100 J=1,NUMAT
            IL=NFIRST(J)
            IU=NLAST(J)
            X=0.0
C$DOIT ASIS
            do 90 K=IL,IU
   90       X=X+C(K,I)**2
  100 SUM1=SUM1+X*X
      if (SUM.GT.EPS.AND.ITER.LT.NITER) go to 20
      write(MIOB,110) ITER,SUM1
  110 format(/10X,'NUMBER OF ITERATIONS =',I4/
     110X,'LOCALISATION VALUE =',F14.9,/)
      write(MIOB,120)
  120 format(3X,'NUMBER OF CENTERS',14X,'(COMPOSITION OF ORBITALS)'//)
      do 150 I=1,NOCC
         SUM=0.D0
         do 140 J=1,NOCC
            CO=0.D0
            do 130 K=1,NORBS
  130       CO=CO+COLD(K,J)*C(K,I)
  140    SUM=SUM+CO*CO*EIG(J)
  150 EIG1(I)=SUM
      do 180 I=1,NOCC
         X=100.D0
         do 160 J=I,NOCC
            if (X.LT.EIG1(J))  go to  160
            X=EIG1(J)
            I1=J
  160    continue
         EIG(I)=EIG1(I1)
         X=EIG1(I1)
         EIG1(I1)=EIG1(I)
         EIG1(I)=X
         do 170 J=1,NORBS
            X=C(J,I1)
            C(J,I1)=C(J,I)
  170    C(J,I)=X
  180 continue
      do 250 I=1,NOCC
         X=0.D0
         do 200 K1=1,NUMAT
            KL=NFIRST(K1)
            KU=NLAST(K1)
            DII=0.D0
            do 190 K=KL,KU
  190       DII=DII+C(K,I)**2
            X=X+DII*DII
  200    PSI1(K1)=DII*100.D0
         X=1/X
         do 220 II=1,NUMAT
            SUM=0.D0
            do 210 J=1,NUMAT
               if (PSI1(J).LT.SUM) go to 210
               SUM=PSI1(J)
               K=J
  210       continue
            PSI1(K)=0.D0
            CII(II)=SUM
            IEL(II)=K
            if (SUM.LT.1.D0) go to 230
  220    continue
  230    continue
         II=II-1
         write(MIOB,240) X,(ELEMNT(NAT(IEL(K))),IEL(K),CII(K),K=1,II)
  240    format(F10.4,4(5(3X,A2,I3,F6.2),/10X))
  250 continue
  260 format(//20X,20H LOCALIZED ORBITALS   ,//)
      write(MIOB,260)
      call MATOUT(C,EIG,NOCC,NORBS,MDIM)
  270 format(10F12.6)
      do 280 I=1,NOCC
         EIG(I)=REFEIG(I)
         do 280 J=1,NORBS
  280 C(J,I)=COLD(J,I)
      return
      end
