      SUBROUTINE ENPART(UHF,H,ALPHA,BETA,P,Q,COORD)
C----------------------------------------------------------*
C
C     NEW SUB. ENPART,  MODIFIED BY TSUNEO HIRANO 1986/6/3/
C
C---------------------------------------------------------*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      PARAMETER (NATMS2 = (NUMATM*(NUMATM+1))/2)
      COMMON /ELEMTS/ ELEMNT(107)
      DIMENSION H(*), ALPHA(*), BETA(*), P(*), Q(*), COORD(3,*)
C--- DEFINED HERE, AND TO BE USED FOR ENPART-PRINT ONLY ---*
      DIMENSION EX(NATMS2,3), EMAT(NUMATM,NUMATM)
C--- END OF DIMENSION DEFINITION ----------------- BY TH --*
      LOGICAL UHF, MINDO3, AM1
      CHARACTER*241 KEYWRD
      CHARACTER*2 ELEMNT
C***********************************************************************
C
C *** ENERGY PARTITIONING WITHIN THE UMINDO/3 AND UMNDO SCHEME
C     ROUTINE WRITTEN BY S.OLIVELLA, BARCELONA NOV. 1979.
C     EXTENDED TO AM1 AND PM3 BY JJPS.
C
C   ON INPUT UHF     = .TRUE. IF A U.H.F. CALCULATION.
C            H       = ONE-ELECTRON MATRIX.
C            ALPHA   = ALPHA ELECTRON DENSITY.
C            BETA    = BETA ELECTRON DENSITY.
C            P       = TOTAL ELECTRON DENSITY.
C            Q       = ATOM ELECTRON DENSITIES.
C
C    NOTHING IS CHANGED ON EXIT.
C
C***********************************************************************
      COMMON /ONELEC/ USS(107), UPP(107), UDD(107)
      COMMON /CORE  / CORE(107)
      COMMON /IDEAS / FN1(107,10),FN2(107,10),FN3(107,10)
      COMMON /ALPHA3/ ALP3(153)
      COMMON /TWOEL3/ F03(107)
      COMMON /ALPHA / ALP(107)
      COMMON /TWOELE/ GSS(107),GSP(107),GPP(107),GP2(107),HSP(107)
     1                ,GSD(107),GPD(107),GDD(107)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM)
     1               ,NLAST(NUMATM),NORBS,NELECS,
     2                NALPHA,NBETA,NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /WMATRX/ W(N2ELEC*2)
      COMMON /KEYWRD/ KEYWRD
      PARAMETER (MDUMY=MAXPAR*MAXPAR-NUMATM*3-NATMS2*4)
      COMMON /SCRACH/ EA(NUMATM,2),EAT(NUMATM), E(NATMS2,4),
     1XDUMY(MDUMY)
      MINDO3=(INDEX(KEYWRD,'MINDO').NE.0)
C@ --------------------------*
C@    CLEAR EMAT(I,I)
      DO 10 I=1,NUMAT
   10 EMAT(I,I)=0.0D0
C@ --------------------------*
C *** RECALCULATE THE DENSITY MATRICES IN THE UHF SCHEME
C
      LINEAR=NORBS*(NORBS+1)/2
      IF( .NOT. UHF) THEN
         DO 20 I=1,LINEAR
   20    BETA(I)=ALPHA(I)
      ENDIF
C
C *** ONE-CENTER ENERGIES
      K=0
      DO 40 I=1,NUMAT
         IA=NFIRST(I)
         IB=NLAST(I)
         NI=NAT(I)
         EA(I,1)=0.0D0
         DO 30 J=IA,IB
            K=K+J
            T=UPP(NI)
            IF(J.EQ.IA) T=USS(NI)
   30    EA(I,1)=EA(I,1)+P(K)*T
         ISS=(IA*(IA+1))/2
         EA(I,2)=0.5D0*GSS(NI)*P(ISS)*P(ISS)
     1  -0.5D0*GSS(NI)*(ALPHA(ISS)*ALPHA(ISS)+BETA(ISS)*BETA(ISS))
         IF(IA.EQ.IB) GO TO 40
         IA1=IA+1
         IA2=IA+2
         IXX=IA1*IA2/2
         IYY=IA2*IB/2
         IZZ=(IB*(IB+1))/2
         IXY=IA1+IA2*IA1/2
         IXZ=IA1+IB*IA2/2
         IYZ=IA2+IB*IA2/2
         ISX=IA+IA1*IA/2
         ISY=IA+IA2*IA1/2
         ISZ=IA+IB*IA2/2
         SS1=P(IXX)*P(IXX)+P(IYY)*P(IYY)+P(IZZ)*P(IZZ)
         SS2=P(ISS)*(P(IXX)+P(IYY)+P(IZZ))
         SS3=P(IXX)*P(IYY)+P(IXX)*P(IZZ)+P(IYY)*P(IZZ)
         SS4=P(ISX)*P(ISX)+P(ISY)*P(ISY)+P(ISZ)*P(ISZ)
         SS5=P(IXY)*P(IXY)+P(IXZ)*P(IXZ)+P(IYZ)*P(IYZ)
         TT1=ALPHA(IXX)*ALPHA(IXX)+ALPHA(IYY)*ALPHA(IYY)
     1+ALPHA(IZZ)*ALPHA(IZZ)+BETA(IXX)*BETA(IXX)
     2+BETA(IYY)*BETA(IYY)+BETA(IZZ)*BETA(IZZ)
         TT2=ALPHA(ISS)*(ALPHA(IXX)+ALPHA(IYY)+ALPHA(IZZ))
     1   +BETA(ISS)*(BETA(IXX)+BETA(IYY)+BETA(IZZ))
         TT3=ALPHA(IXX)*ALPHA(IYY)+ALPHA(IXX)*ALPHA(IZZ)
     1+ALPHA(IYY)*ALPHA(IZZ)+BETA(IXX)*BETA(IYY)
     2+BETA(IXX)*BETA(IZZ)+BETA(IYY)*BETA(IZZ)
         TT4=ALPHA(ISX)*ALPHA(ISX)+ALPHA(ISY)*ALPHA(ISY)
     1+ALPHA(ISZ)*ALPHA(ISZ)+BETA(ISX)*BETA(ISX)
     2+BETA(ISY)*BETA(ISY)+BETA(ISZ)*BETA(ISZ)
         TT5=ALPHA(IXY)*ALPHA(IXY)+ALPHA(IXZ)*ALPHA(IXZ)
     1+ALPHA(IYZ)*ALPHA(IYZ)+BETA(IXY)*BETA(IXY)
     2+BETA(IXZ)*BETA(IXZ)+BETA(IYZ)*BETA(IYZ)
         EA(I,2)=EA(I,2)+0.5D0*GPP(NI)*SS1+GSP(NI)*SS2
     1+GP2(NI)*SS3+HSP(NI)*SS4*2.0D0+0.5D0*(GPP(NI)-GP2(NI))*SS5*2.0D0
     2                -0.5D0*GPP(NI)*TT1-GSP(NI)*TT4-GP2(NI)*TT5-
     3        HSP(NI)*(TT2+TT4)-0.5D0*(GPP(NI)-GP2(NI))*(TT3+TT5)
   40 CONTINUE
      AM1=(INDEX(KEYWRD,'AM1')+INDEX(KEYWRD,'PM3').NE.0)
      IF(MINDO3) THEN
         WRITE(6,'(///,10X,''TOTAL ENERGY PARTITIONING IN MINDO/3'')')
      ELSEIF( INDEX(KEYWRD,'PM3').NE.0 ) THEN
         WRITE(6,'(///,10X,''TOTAL ENERGY PARTITIONING IN PM3'')')
      ELSEIF( INDEX(KEYWRD,'AM1').NE.0 ) THEN
         WRITE(6,'(///,10X,''TOTAL ENERGY PARTITIONING IN AM1'')')
      ELSE
         WRITE(6,'(///,10X,''TOTAL ENERGY PARTITIONING IN MNDO'')')
      ENDIF
      WRITE(6,'(/10X,''ALL ENERGIES ARE IN ELECTRON VOLTS'')')
      KL=0
   50 K=KL+1
      KL=KL+10
      KL=MIN0(KL,NUMAT)
      DO 60 I=K,KL
   60 EAT(I)=EA(I,1)+EA(I,2)
      IF(NUMAT.GT.KL) GO TO 50
   70 EAU=0.0D0
      EAE=0.0D0
      DO 80 I=1,NUMAT
         EAU=EAU+EA(I,1)
   80 EAE=EAE+EA(I,2)
      TONE=EAU+EAE
C *** TWO-CENTER ENERGIES
C     RESONANCE (E(N,1)) TERMS
      N=1
      DO 100 II=2,NUMAT
         E(N,1)=0.0D0
         IA=NFIRST(II)
         IB=NLAST(II)
         IMINUS=II-1
         ONEII=1.D0
         IF(NAT(II).EQ.102)ONEII=0.D0
         DO 90 JJ=1,IMINUS
            N=N+1
            JA=NFIRST(JJ)
            JB=NLAST(JJ)
            ONEJJ=1.D0
            IF(NAT(JJ).EQ.102)ONEJJ=0.D0
            E(N,1)=0.0D0
            DO 90 I=IA,IB
               KA=(I*(I-1))/2
               DO 90 K=JA,JB
                  IK=KA+K
   90    E(N,1)=E(N,1)+2.0D0*P(IK)*H(IK)*ONEII*ONEJJ
  100 N=N+1
C
C     THE CODE THAT FOLLOWS APPLIES ONLY TO MNDO
C
      IF(.NOT.MINDO3) THEN
C     CORE-CORE REPULSION (E(N,2)) AND CORE-ELEC. ATTRACTION (E(N,3)).
         N=1
         KK=0
         DO 180 II=2,NUMAT
            E(N,2)=0.0D0
            E(N,3)=0.0D0
            IA=NFIRST(II)
            IB=NLAST(II)
            NI=NAT(II)
            ISS=(IA*(IA+1))/2
            IMINUS=II-1
            DO 170 JJ=1,IMINUS
               N=N+1
               JA=NFIRST(JJ)
               JB=NLAST(JJ)
               NJ=NAT(JJ)
               JSS=(JA*(JA+1))/2
               KK=KK+1
               G=W(KK)
               R=SQRT((COORD(1,II)-COORD(1,JJ))**2+(COORD(2,II)-COORD(2,
     1JJ))**2+  (COORD(3,II)-COORD(3,JJ))**2)
               SCALE=1.0D0+EXP(-ALP(NI)*R)+EXP(-ALP(NJ)*R)
               NT=NI+NJ
               IF(NT.LT.8.OR.NT.GT.9) GO TO 110
               IF(NI.EQ.7.OR.NI.EQ.8) SCALE=SCALE+(R-1.0D0)*EXP(-ALP(NI)
     1*R)
               IF(NJ.EQ.7.OR.NJ.EQ.8) SCALE=SCALE+(R-1.0D0)*EXP(-ALP(NJ)
     1*R)
  110          E(N,2)=CORE(NI)*CORE(NJ)*G*SCALE
               IF( AM1 )THEN
                  SCALE=0.0D0
                  DO 120 IG=1,10
                     IF(ABS(FN1(NI,IG)).GT.0.D0)
     1SCALE=SCALE +CORE(NI)*CORE(NJ)/R*
     2FN1(NI,IG)*EXP(-FN2(NI,IG)*(R-FN3(NI,IG))**2)
                     IF(ABS(FN1(NJ,IG)).GT.0.D0)
     1SCALE=SCALE +CORE(NI)*CORE(NJ)/R*
     2FN1(NJ,IG)*EXP(-FN2(NJ,IG)*(R-FN3(NJ,IG))**2)
  120             CONTINUE
                  E(N,2)=E(N,2)+SCALE
               ENDIF
               E(N,3)=-(P(ISS)*CORE(NJ)+P(JSS)*CORE(NI))*G
               IF(NJ.LT.3) GO TO 140
               KINC=9
               JAP1=JA+1
               DO 130 K=JAP1,JB
                  KC=(K*(K-1))/2
                  DO 130 L=JA,K
                     KL=KC+L
                     BB=2.0D0
                     IF(K.EQ.L) BB=1.0D0
                     KK=KK+1
  130          E(N,3)=E(N,3)-P(KL)*CORE(NI)*BB*W(KK)
               GO TO 150
  140          KINC=0
  150          IF(NI.LT.3) GO TO 170
               IAP1=IA+1
               DO 160 I=IAP1,IB
                  KA=(I*(I-1))/2
                  DO 160 J=IA,I
                     IJ=KA+J
                     AA=2.0D0
                     IF(I.EQ.J) AA=1.0D0
                     KK=KK+1
                     E(N,3)=E(N,3)-P(IJ)*CORE(NJ)*AA*W(KK)
  160          KK=KK+KINC
  170       CONTINUE
  180    N=N+1
C     COULOMB (E(N,4)) AND EXCHANGE (EX(N)) TERMS
         N=1
         KK=0
         DO 200 II=2,NUMAT
            E(N,4)=0.0D0
            EX(N,1)=0.0D0
            IA=NFIRST(II)
            IB=NLAST(II)
            IMINUS=II-1
            DO 190 JJ=1,IMINUS
               JA=NFIRST(JJ)
               JB=NLAST(JJ)
               N=N+1
               E(N,4)=0.0D0
               EX(N,1)=0.0D0
               DO 190 I=IA,IB
                  KA=(I*(I-1))/2
                  DO 190 J=IA,I
                     KB=(J*(J-1))/2
                     IJ=KA+J
                     AA=2.0D0
                     IF(I.EQ.J) AA=1.0D0
                     PIJ=P(IJ)
                     DO 190 K=JA,JB
                        KC=(K*(K-1))/2
                        IK=KA+K
                        JK=KB+K
                        DO 190 L=JA,K
                           IL=KA+L
                           JL=KB+L
                           KL=KC+L
                           BB=2.0D0
                           IF(K.EQ.L) BB=1.0D0
                           KK=KK+1
                           G=W(KK)
                           E(N,4)=E(N,4)+AA*BB*G*PIJ*P(KL)
  190       EX(N,1) = EX(N,1)
     1    -0.5D0*AA*BB*G*(ALPHA(IK)*ALPHA(JL)+ALPHA(IL)*ALPHA(JK)+
     2    BETA(IK)*BETA(JL)+BETA(IL)*BETA(JK))
  200    N=N+1
      ELSE
         N=1
         DO 260 I=2,NUMAT
            E(N,2)=0.0D0
            E(N,3)=0.0D0
            E(N,4)=0.0D0
            EX(N,1)=0.0D0
            IA=NFIRST(I)
            IB=NLAST(I)
            NI=NAT(I)
            IMINUS=I-1
            DO 250 J=1,IMINUS
               N=N+1
               JA=NFIRST(J)
               JB=NLAST(J)
               NJ=NAT(J)
               RIJ=(COORD(1,I)-COORD(1,J))**2+(COORD(2,I)-COORD(2,J))**2
     1+  (COORD(3,I)-COORD(3,J))**2
               GIJ=14.399D0/SQRT(RIJ+(7.1995D0/F03(NI)+7.1995D0/F03(NJ))
     1**2)
               PAB2=0.0D0
               IJ=MAX(NI,NJ)
               NBOND=(IJ*(IJ-1))/2+NI+NJ-IJ
               RIJ=SQRT(RIJ)
               IF(NBOND.EQ.22 .OR. NBOND .EQ. 29) GO TO 210
               GO TO 220
  210          SCALE=ALP3(NBOND)*EXP(-RIJ)
               GO TO 230
  220          SCALE=EXP(-ALP3(NBOND)*RIJ)
  230          CONTINUE
               E(N,2)=CORE(NI)*CORE(NJ)*GIJ+
     1     ABS(CORE(NI)*CORE(NJ)*(14.399D0/RIJ-GIJ)*SCALE)
               E(N,3)=(-Q(I)*CORE(NJ)-Q(J)*CORE(NI))*GIJ
               E(N,4)=Q(I)*Q(J)*GIJ
               DO 240 K=IA,IB
                  KK=(K*(K-1))/2
                  DO 240 L=JA,JB
                     LK=KK+L
  240          PAB2=PAB2+ALPHA(LK)*ALPHA(LK)+BETA(LK)*BETA(LK)
  250       EX(N,1) = -PAB2*GIJ
  260    N=N+1
      ENDIF
      NUMAT1=(NUMAT*(NUMAT+1))/2
      DO 270 I=1,4
  270 E(NUMAT1,I)=0.0D0
      DO 280 I=1,3
  280 EX(NUMAT1,I)=0.0D0
C@ --------------------------*
C-----PRINT OUT ONE AND TWO CENTER ENERGIES
C
C     E(I,1):     RESONANCE ENERGY
C     E(I,2):     NUCLEAR-NUCLEAR REPULSION ENERGY
C     E(I,3):     ELECTRON-NUCLEAR ATTRACTION ENERGY
C     E(I,4):     ELECTRON-ELECTRON REPULSION ENERGY
C     EX(I,1):    EXCHANGE  ENERGY
C     EX(I,2):    EXCHANGE + RESONANCE ENERGY
C#      WRITE(6,'(//,''       ONE AND TWO CENTER ENERGIES (EV) '')')
C
C#      WRITE(6,'(/,''  [RESONANCE TERM] (EV)'')')
C#      CALL VECPRT(E,NUMAT)
C
C#      WRITE(6,'(/,''  [EXCHANGE TERM] (EV)'')')
C#      CALL VECPRT(EX,NUMAT)
C
C#      WRITE(6,'(/,''  [RESONANCE + EXCHANGE] (EV)'')')
      DO 290 N=1,NUMAT1
  290 EX(N,2) =E(N,1) + EX(N,1)
C
C   ADD IN MONOCENTRIC EXCHANGE AND COULOMBIC TERM
C
      DO 300 I=1,NUMAT
  300 EX((I*(I+1))/2,2)=EA(I,2)
C#      CALL VECPRT(EX(1,2),NUMAT)
C
C#      WRITE(6,'(/,''  [ELECTRON - ELECTRON REPULSION] (EV)'')')
C#      CALL VECPRT(E(1,4),NUMAT)
C
C#      WRITE(6,'(/,''  [ELECTRON-NUCLEAR ATTRACTION] (EV)'')')
      DO 310 I=1,NUMAT
  310 E((I*(I+1))/2,3)=EA(I,1)
C#      CALL VECPRT(E(1,3),NUMAT)
C
C#      WRITE(6,'(/,''  [NUCLEAR-NUCLEAR REPULSION] (EV)'')')
C#      CALL VECPRT(E(1,2),NUMAT)
C
      DO 320 N=1,NUMAT1
  320 EX(N,3) =E(N,4) + E(N,3) + E(N,2)
C     PRINT OUT OF TOTAL COULOMB TERM
C#      WRITE(6,'(/,''  [TOTAL COULOMB TERM (E-E, E-N, AND N-N)] (EV)'')
C#      CALL VECPRT(EX(1,3),NUMAT)
C     PRINT OUT OF TWO-CENTER SUM(OFF-DIAGONAL) +
C                  ONE-CENTER SUM(DIAGONAL).
C#      WRITE(6,'(/,''  [TWO-CENTER SUM (OFF-DIAGONAL), AND  '',
C#     1''ONE-CENTER SUM (DIAGONAL)] (EV)'')')
C#      DO 340 N=1,NUMAT1
C#  340 EX(N,3)=EX(N,3)+EX(N,2)
C#      CALL VECPRT(EX(1,3),NUMAT)
      WRITE(6,'(/,8(10X,A,/))')
     1'  ONE-CENTER TERMS',' ',
     2'E-E:  ELECTRON-ELECTRON REPULSION',
     3'E-N:  ELECTRON-NUCLEAR ATTRACTION'
      WRITE(6,'(/,''   ATOM     E-E       E-N    (E-E + E-N)'')')
      DO 330 I=1,NUMAT
         J=(I*(I+1))/2
         WRITE(6,'(2X,A2,I3,1X,2F10.4,F10.4)')
     1ELEMNT(NAT(I)),I, EX(J,2), E(J,3), EX(J,2)+E(J,3)
  330 CONTINUE
      WRITE(6,'(/,8(10X,A,/))')
     1'    TWO-CENTER TERMS',' ',
     2'J:   RESONANCE ENERGY          E-E: ELECTRON-ELECTRON REPULSION',
     3'K:   EXCHANGE ENERGY           E-N: ELECTRON-NUCLEAR ATTRACTION',
     4'                               N-N: NUCLEAR-NUCLEAR REPULSION',
     5'C:   COULOMBIC INTERACTION = E-E + E-N + N-N',
     6'EE:  TOTAL OF ELECTRONIC AND NUCLEAR ENERGIES'
      WRITE(6,'(/,''   ATOM          J        K       E-'',
     1''E       E-N      N-N      C        EE'')')
      WRITE(6,'(''   PAIR'')')
      IJ=0
      DO 340 I=1,NUMAT
         DO 340 J=1,I
            IJ=IJ+1
            IF(I.NE.J) THEN
               WRITE(6,'(1X,A2,I3,1X,A2,I3,1X,2F9.4,F9.4,F10.4,F9.4,F8.4
     1,F9.4)')ELEMNT(NAT(I)),I,ELEMNT(NAT(J)),J, E(IJ,1), EX(IJ,1),
     2E(IJ,4), E(IJ,3), E(IJ,2), EX(IJ,3), EX(IJ,2)+EX(IJ,3)
            ELSE
               IF(I.LT.6.OR.I.EQ.NUMAT)THEN
                  WRITE(6,*)
               ELSE
                  WRITE(6,'(/,''   ATOM          J        K       E-'',
     1''E       E-N      N-N      C        EE'')')
                  WRITE(6,'(''   PAIR'')')
               ENDIF
            ENDIF
  340 CONTINUE
C
C     ++++   TOTALS   ++++
C
      EABR=0.0D0
      EABX=0.0D0
      EABEE=0.0D0
      EABEN=0.0D0
      EABNN=0.0D0
      DO 350 I=1,NUMAT
  350 E((I*(I+1))/2,3)=0.D0
      DO 360 I=1,NUMAT1
         EABR=EABR+E(I,1)
         EABX=EABX+EX(I,1)
         EABEE=EABEE+E(I,4)
         EABEN=EABEN+E(I,3)
         EABNN=EABNN+E(I,2)
  360 CONTINUE
      EABRX=EABR+EABX
      EABE=EABEE+EABEN+EABNN
      TTWO=EABRX+EABE
      ET=TONE+TTWO
C@ ***************************************************************
      WRITE(6,370)
  370 FORMAT(///,'***  SUMMARY OF ENERGY PARTITION  ***')
      WRITE(6,380)
  380 FORMAT(1H ,'---------------------------------------')
      WRITE(6,'(''     ONE-CENTER TERMS'')')
      WRITE(6,390) EAU
  390 FORMAT(/,' ELECTRON-NUCLEAR  (ONE-ELECTRON) ',F17.4,' EV')
      WRITE(6,400) EAE
  400 FORMAT(' ELECTRON-ELECTRON (TWO-ELECTRON) ',F17.4,' EV')
      WRITE(6,410) TONE
  410 FORMAT(/,' TOTAL OF ONE-CENTER TERMS        ',18X,F15.4,' EV')
      WRITE(6,380)
      WRITE(6,'(''     TWO-CENTER TERMS'')')
      WRITE(6,420) EABR
  420 FORMAT(/,' RESONANCE ENERGY',8X,F15.4,' EV')
      WRITE(6,430) EABX
  430 FORMAT(' EXCHANGE ENERGY ',8X,F15.4,' EV')
      WRITE(6,440) EABRX
  440 FORMAT(/,' EXCHANGE + RESONANCE ENERGY:       ',F15.4,' EV')
      WRITE(6,450) EABEE
  450 FORMAT(/,' ELECTRON-ELECTRON REPULSION',F12.4,' EV')
      WRITE(6,460) EABEN
  460 FORMAT(  ' ELECTRON-NUCLEAR ATTRACTION',F12.4,' EV')
      WRITE(6,470) EABNN
  470 FORMAT(  ' NUCLEAR-NUCLEAR REPULSION  ',F12.4,' EV')
      WRITE(6,480) EABE
  480 FORMAT(/,' TOTAL ELECTROSTATIC INTERACTION    ',F15.4,' EV',/)
      WRITE(6,490) TTWO
  490 FORMAT(' GRAND TOTAL OF TWO-CENTER TERMS   ',17X,F15.4,' EV')
      WRITE(6,380)
      WRITE(6,500) ET
  500 FORMAT(' ETOT (EONE + ETWO)   ',30X,F15.4,' EV'//)
      RETURN
      END
