C  7 APR 06 - MWS - REMOVE NON-ARGUMENT ADJUSTABLE ARRAYS
C 29 MAR 06 - FLG - LIBRARY ROUTINES FOR ELONGATION METHOD
C
C*MODULE ELGLIB  *DECK ELGINP
      SUBROUTINE ELGINP
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL LCUT,OPTELG,SKPCUT,ATOB,NSPACE,GOPARR,DSKWRK,MASWRK
C
      COMMON /ELGPMT/ NELONG,NATM,NASPIN,NCT,NBNDAB,NTMLB,IPRI,LDOS
      COMMON /ELGRST/ IRSTRT,I2EA,IGOOD
      COMMON /ELGIDX/ LCUT
      COMMON /ELGOPT/ OPTELG
      COMMON /CTSHFT/ ATOB
      COMMON /ECTSKP/ SKPCUT
      COMMON /LOCSPC/ NSPACE
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (NNAM=12)
      DIMENSION ELGNAM(NNAM),KELGNM(NNAM)
      DATA ELG /8HELG     /
      DATA ELGNAM/8HNELONG  ,8HNATM    ,8HNASPIN  ,8HNCT     ,
     *            8HNBNDAB  ,8HNTMLB   ,8HIPRI    ,8HLDOS    ,
     *            8HI2EA    ,8HSKPCUT  ,8HATOB    ,8HNSPACE  /
      DATA KELGNM /1,1,1,1,  1,1,1,1,  1,0,0,0/
C
C   ----- INITIALIZE VARIABLES FOR NAMELIST $ELG -----
C
      NELONG = 0                 ! INITIALIZE ELG CONTRL
      NATM   = 0                 ! NUMBER OF ATOMS IN A-REGION
      NASPIN = 0                 ! MULTIPLICITY A-REGION
      NBNDAB = 0                 ! NUMBER OF BONDS BETWEEN A AND B
      NTMLB  = 0                 ! NUMBER OF TERMINAL ATOMS IN B-REGION
      NCT    = 0                 ! CONTROLLER FOR AO-CUT
      IPRI   = 0                 ! PRINT LEVEL
      LDOS   = 0                 ! LDOS CALCULATION
      I2EA   = 0                 ! READ-IN 2E-INTEGRALS FOR A-REGION
      LCUT   = .FALSE.
      OPTELG = .FALSE.
      SKPCUT = .FALSE.
      ATOB = .TRUE.
      NSPACE = .FALSE.
C
      JRET = 0
      CALL NAMEIO(IR,JRET,ELG,NNAM,ELGNAM,KELGNM,
     *            NELONG,NATM,NASPIN,NCT,NBNDAB,NTMLB,IPRI,LDOS,I2EA,
     *            SKPCUT,ATOB,NSPACE,
     *            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,   0,0)
C
C
      IF(JRET .EQ. 2) THEN
         IF (MASWRK) WRITE (IW,'(''TYPING ERROR IN $ELG INPUT'')')
         CALL ABRT
      END IF
      IF (NELONG.GE.1.AND.MASWRK) THEN
         WRITE(IW,*)'READ IN ELG PARAMETERS:'
         WRITE(IW,*)'NELONG,NATM,NASPIN,NCT,NBNDAB,NTMLB,IPRI,LDOS,I2EA'
         WRITE(IW,*) NELONG,NATM,NASPIN,NCT,NBNDAB,NTMLB,IPRI,LDOS,I2EA
      END IF
      IF(NELONG.LE.1) THEN      ! SET BACK TO DEFAULTS IF NELONG<=1
         NATM   = 0
         NASPIN = 0
         NBNDAB = 0
         NTMLB  = 0
         NCT    = 0
         IPRI   = 0
         LDOS   = 0
         I2EA = 0
         LCUT = .FALSE.
         ATOB = .FALSE.
      ELSE
         DSKWRK = .TRUE.
         CALL ELGSTRT
         IF(NELONG.EQ.2) I2EA = 0
         IF(NCT.GT.0) CALL CTDONE(LCUT)
      ENDIF
      RETURN
      END
C*MODULE ELGLIB  *DECK ELGSTRT
      SUBROUTINE ELGSTRT
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /ELGFLS /NEG0,NEG1,NEG2,NEG3,NEG4,NEG5,NEG6,NEG7,NEG8,NEG9
C
      NEG0 = 70
      NEG1 = 71
      NEG2 = 72
      NEG4 = 74
      NEG5 = 75
      NEG7 = 77
      NEG8 = 78
      NEG9 = 79
       NNN = 80
C
C      CALL SEQOPN(NEG0,'ELGDOS','UNKNOWN',.FALSE.,'FORMATTED')
      CALL SEQOPN(NEG1,'ELGDAT','UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQOPN(NEG2,'ELGPAR', 'UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQOPN(NEG4,'ELGCUT', 'UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQOPN(NEG5,'ELGVEC', 'UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQOPN(NEG7,'EGINTA', 'UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQOPN(NEG8,'EGINTB', 'UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQOPN(NEG9,'EGTDHF', 'UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQOPN(NNN, 'EGTEST', 'UNKNOWN',.FALSE.,'UNFORMATTED')
C
      RETURN
      END
C
C*MODULE ELGLIB  *DECK LINEAR
      SUBROUTINE LINEAR(A,B,N,NMAX)
C
C     PLACES SYMMETRIC SQUARE ARRAY IN LINEAR FORM
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(NMAX,NMAX),B(*)
C
      K=1
      DO J=1,N
         DO I=1,J
            B(K)=A(I,J)
            K=K+1
         END DO
      END DO
      RETURN
      END
C
C*MODULE ELGLIB  *DECK WRITN4
      SUBROUTINE WRITN4(N,N1,A,M1,M2)
C
C     PRINTING SUBROUTINE OF MATRIX
C     WITH DIFFERENT NUMBERS OF COLUMN AND ROW
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(M2,M1)
C
    7 FORMAT(1H )
    8 FORMAT(1H ,I3,2X,10F10.5)
    9 FORMAT(1H ,5X,10(2X,I3,5X))
      INDEX=(N-1)/10+1
      DO 20 I=1,INDEX
      I1=10*(I-1)+1
      I2=10*I
      IF(N.LT.I2) I2=N
      WRITE(6,9) (ID,ID=I1,I2)
      WRITE(6,7)
      DO 20 J=1,N1
      WRITE(6,8) J,(A(J,ID),ID=I1,I2)
   20 CONTINUE
      RETURN
      END
C
C*MODULE ELGLIB  *DECK NATCUT
      INTEGER FUNCTION NATCUT(NBF1,NBF2)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      CHARACTER NAME*8
C
      NCUT = 1
      WRITE(NAME,'(A8)') BFLAB(NBF1)
      READ (NAME,'(2X,I3,3X)') NCUT1
      DO 10 I=NBF1+1,NBF2
         WRITE(NAME,'(A8)') BFLAB(I)
         READ (NAME,'(2X,I3,3X)') NCUT2
         IF(NCUT2.NE.NCUT1) THEN
            NCUT = NCUT + 1
            NCUT1 = NCUT2
         ENDIF
 10   CONTINUE
      NATCUT = NCUT
C
      RETURN
      END
C
C*MODULE ELGLIB  *DECK JKFRZ
      SUBROUTINE JKFRZ(MAT,LX,LP,LU,LMAX)
C
C***********************************************************************
C                 ROUTINE WRITTEN BY JACEK KORCHOWIEC
C     ALLOWS READING IN SIMPLE WAY VARIABLES TO INTEGER MATRIX
C       (E.G. FROZEN COORDINATES FOR GEOMETRY OPTIMIZATION RUN)
C   MAT  - INTEGER, ONE-DIMENSION MATRIX
C   LX   - THE NUMBER OF READ IN VARIABLES
C   LP   - FIRST POSITION IN MATRIX MAT (USUALLY ONE)
C   LU   - THE EXTERNAL READ UNIT NUMBER
C   LMAX - THE MAT MATRIX DIMENSION
C
C E.G. INSTEAD OF           1,2,3,4,5,6,7,8,9,10,11,12,24,33,34,35,36,37
C      IT CAN BE WRITTEN                1-12 24, 35-37
C RUTINE REQUIRES $JKSTPT ... $END GROUP
C***********************************************************************
C
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      INTEGER MAT(LMAX)
      CHARACTER*80 ALFA
      LOGICAL LOGIK,LFIRST,LAST
C
      NUMER=-10
      LICZ2=LP
      LICZ3=0
      LFIRST=.FALSE.
      LAST=.FALSE.
      LSTRT = 1
C
      REWIND(LU)
C
 10   CONTINUE
      READ(LU,'(A80)',ERR=500,END=600) ALFA
      K=INDEX(ALFA,'$JKSTPT')
      L=INDEX(ALFA,'$JKSTPT')
      IF(K.NE.0 .OR. L.NE.0) THEN
         IF(K.NE.0) LSTRT = K+7
         IF(L.NE.0) LSTRT = L+7
         K=INDEX(ALFA,'$END')
         L=INDEX(ALFA,'$END')
         IF(K.NE.0) THEN
            LEND = K-1
            LFIRST = .TRUE.
         ELSEIF(L.NE.0) THEN
            LEND = L-1
            LFIRST = .TRUE.
         ELSE
            LEND = 80
         ENDIF
         GOTO 20
      ENDIF
      GOTO 10
C
 20   CONTINUE
C
      DO 50 I=LSTRT,LEND
C
         LOGIK=.TRUE.
         IF(ALFA(I:I).EQ.',') LOGIK=.FALSE.
         IF(ALFA(I:I).EQ.' ') LOGIK=.FALSE.
         IF(ALFA(I:I).EQ.'-') THEN
            LOGIK=.FALSE.
            NUMER=LICZ2
            IF(ALFA(I-1:I-1).EQ.' ') NUMER=LICZ2-1
         ENDIF
C
         IF(LOGIK) THEN
            LICZ3=LICZ3+1
C
         ELSE
            IF(LICZ3.EQ.0) GO TO 40
            IF(LICZ3.EQ.1) THEN
               READ(ALFA(I-LICZ3:I-1),100,ERR=500) MAT(LICZ2)
            ELSE IF(LICZ3.EQ.2) THEN
               READ(ALFA(I-LICZ3:I-1),120,ERR=500) MAT(LICZ2)
            ELSE IF(LICZ3.EQ.3) THEN
               READ(ALFA(I-LICZ3:I-1),140,ERR=500) MAT(LICZ2)
            ELSE IF(LICZ3.EQ.4) THEN
               READ(ALFA(I-LICZ3:I-1),160,ERR=500) MAT(LICZ2)
            ENDIF
C
            IF(NUMER.EQ.(LICZ2-1)) THEN
               K1=MAT(NUMER)
               K2=MAT(LICZ2)
               LICZ2=NUMER
C
               DO 30 J=K1,K2
                  MAT(LICZ2)=J
                  LICZ2=LICZ2+1
 30            CONTINUE
               LICZ2=LICZ2-1
            ENDIF
C
            LICZ2=LICZ2+1
            LICZ3=0
 40         CONTINUE
         ENDIF
 50   CONTINUE
C
      LICZ2=LICZ2-1
C
      LX=LICZ2-LP+1
      LICZ2=LICZ2+1
      LSTRT=1
C
      IF(LAST) GOTO 600
      IF(LFIRST) THEN
         GOTO 600
      ELSE
         READ(LU,'(A80)',ERR=500) ALFA
         K=INDEX(ALFA,'$END')
         L=INDEX(ALFA,'$END')
         IF(K.NE.0 .OR. L.NE.0) THEN
            IF(K.NE.0) LEND=K-1
            IF(L.NE.0) LEND=L-1
            LAST=.TRUE.
         ENDIF
      ENDIF
      GOTO 20
C
 500  CONTINUE
      IF(MASWRK) THEN
      WRITE(6,*)' CHECK YOUR INPUT DATA FILE'
      WRITE(6,*)' ERROR OCCURED DURING READING'
      WRITE(6,*)' GROUP $JKSTPT ..... $END   '
      WRITE(6,*)'       !!!!!!!!!!!!!        '
      ENDIF
C
 600  CONTINUE
      REWIND(LU)
      RETURN
C
 100  FORMAT(I1)
 120  FORMAT(I2)
 140  FORMAT(I3)
 160  FORMAT(I4)
      END
C
C*MODULE ELGLIB  *DECK ELGZMT
      SUBROUTINE ELGZMT(NAT,NEG2,NTMLB,ZMAT,ZWRK)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL LINEAR
C
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      DIMENSION ZMAT(3*NAT),ZWRK(NAT)
C
      DATA OPT/8HOPTIMIZE/
      DATA ANGS/0.52917724924D+00/
      DATA PI/3.14159265359D+00/ , ONE80/180.0D+00/
C
      IF(RUNTYP.NE.OPT) RETURN
      IF(NZVAR.EQ.0) RETURN
      CALL DAREAD(IDAF,IODA,ZMAT,NZVAR,47,0)
C
      ZWRK(1)=ZMAT(1)*ANGS
      ZWRK(2)=ZMAT(2)*ANGS
      II=4
      DO 10 I=3,NAT-NTMLB-1
         ZWRK(I)=ZMAT(II)*ANGS
         II=II+3
 10   CONTINUE
C
      CALL SEQREW (NEG2)
      WRITE(NEG2) NAT-NTMLB-1
      WRITE(NEG2) (ZWRK(I),I=1,NAT-NTMLB-1)
C
      ZWRK(1)=ZMAT(3)*ONE80/PI
      II=5
      DO 20 I=2,NAT-NTMLB-2
         ZWRK(I)=ZMAT(II)*ONE80/PI
         II=II+3
 20   CONTINUE
      WRITE(NEG2) (ZWRK(I),I=1,NAT-NTMLB-2)
      II=6
      DO 30 I=1,NAT-NTMLB-3
         ZWRK(I)=ZMAT(II)*ONE80/PI
         II=II+3
 30   CONTINUE
      WRITE(NEG2) (ZWRK(I),I=1,NAT-NTMLB-3)
C
      RETURN
      END
C
C*MODULE ELGLIB  *DECK ELGZRD
      SUBROUTINE ELGZRD(BL,ALPHA,BETA)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL LINEAR
C
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
      COMMON /ELGFLS /NEG0,NEG1,NEG2,NEG3,NEG4,NEG5,NEG6,NEG7,NEG8,NEG9
C
      DIMENSION BL(*),ALPHA(*),BETA(*)
C
      IF(NZVAR.EQ.0) RETURN
C
      CALL SEQREW (NEG2)
      READ(NEG2) NN
      READ(NEG2) (BL(I),I=2,NN+1)
      READ(NEG2) (ALPHA(I),I=3,NN+1)
      READ(NEG2) (BETA(I),I=4,NN+1)
C
      RETURN
      END
C
C*MODULE ELGLIB  *DECK NEWSCHMELG
      SUBROUTINE NEWSCHMELG(V,S,WRK,N,M,LDV,MOCC1)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION V(LDV,M),S(*),WRK(N)
      PARAMETER (SMALL=1.0D-15)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     GRAMM-SCHMIDT ORTHOGONALISATION OF -M- VECTORS -V- OF SIZE -N-
C     ON A METRIC -S- (I.E., OVERLAP INTEGRALS).
C
      DO 30 I=MOCC1+1,M
         CALL MTARBR(S,N,V(1,I),1,WRK,N,1)
         DO 10 J=1,I-1
            CALL DAXPY(N,-DDOT(N,WRK,1,V(1,J),1),V(1,J),1,V(1,I),1)
 10      CONTINUE
        CALL MTARBR(S,N,V(1,I),1,WRK,N,1)
         Q=SQRT(ABS(DDOT(N,WRK,1,V(1,I),1)))
         IF(Q.LT.SMALL) THEN
            IF(MASWRK)WRITE(6,*) 'TOO SMALL NORM IN SCHMIDT',Q
            CALL ABRT
         END IF
         CALL DSCAL(N,1.0D+00/Q,V(1,I),1)
 30   CONTINUE
C
      RETURN
      END
C
C*MODULE ELGLIB  *DECK SCHMELG
      SUBROUTINE SCHMELG(V,S,WRK,N,M,LDV,MOCC1,MOCC,MVAC1)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION V(LDV,M),S(*),WRK(N)
      PARAMETER (SMALL=1.0D-15)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     GRAMM-SCHMIDT ORTHOGONALISATION OF -M- VECTORS -V- OF SIZE -N-
C     ON A METRIC -S- (I.E., OVERLAP INTEGRALS).
C
      DO 30 I=MOCC1+1,M
         IF(I.GT.MOCC.AND.I.LE.(MOCC+MVAC1)) GOTO 20
         CALL MTARBR(S,N,V(1,I),1,WRK,N,1)
         DO 10 J=1,I-1
            CALL DAXPY(N,-DDOT(N,WRK,1,V(1,J),1),V(1,J),1,V(1,I),1)
 10      CONTINUE
        CALL MTARBR(S,N,V(1,I),1,WRK,N,1)
         Q=SQRT(ABS(DDOT(N,WRK,1,V(1,I),1)))
         IF(Q.LT.SMALL) THEN
            IF(MASWRK)WRITE(6,*) 'TOO SMALL NORM IN SCHMIDT',Q
            CALL ABRT
         END IF
         CALL DSCAL(N,1.0D+00/Q,V(1,I),1)
 20      CONTINUE
 30   CONTINUE
C
      RETURN
      END
C
C*MODULE ELGLIB  *DECK ELGCRT
      SUBROUTINE ELGCRT(NEG2,NAT,NTMLB,NATCUT,C)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL LINEAR
C
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      DIMENSION C(3,NAT)
C
      DATA OPT/8HOPTIMIZE/
      DATA ANGS/0.52917724924D+00/
C
      IF(RUNTYP.NE.OPT) RETURN
      IF(NZVAR.NE.0) RETURN
      CALL DAREAD(IDAF,IODA,C,3*NAT,1,0)
C
      DO 10 I=1,NAT
         C(1,I)=C(1,I)*ANGS
         C(2,I)=C(2,I)*ANGS
         C(3,I)=C(3,I)*ANGS
 10   CONTINUE
C
      CALL SEQREW (NEG2)
      WRITE(NEG2) NAT-NTMLB-NATCUT
      WRITE(NEG2) (C(1,I),I=NATCUT+1,NAT-NTMLB)
      WRITE(NEG2) (C(2,I),I=NATCUT+1,NAT-NTMLB)
      WRITE(NEG2) (C(3,I),I=NATCUT+1,NAT-NTMLB)
C
      RETURN
      END
C*MODULE ELGLIB  *DECK ELGDOS
      SUBROUTINE ELGDOS(EIGS,C,S,DEN,DENATM,DENTOT,NAO,NMO,NPTS,X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXATM=2000, MXAO=8192)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /ELGLAB/ IATM(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),CC(3,MXATM),IAN(MXATM)
      COMMON /ELGFLS /NEG0,NEG1,NEG2,NEG3,NEG4,NEG5,NEG6,NEG7,NEG8,NEG9
C*******************************************************************
C
C     CALCULATE THE LOCAL DENSITY OF STATES FOR APERIODIC SYSTEMS
C
C             (INTERPOLATION USING GAUSSIAN TYPE FUNCTIONS)
C                    (CONVENTIONAL CALCULATION ONLY)
C
C     ON ENTRY:  EIGS,C: EIGENVALUES AND EIGENVECTORS,
C                     S: OVERLAP MATRIX
C                NAO >= NMO
C
C*******************************************************************
      DIMENSION C(NAO,NMO),EIGS(NAO),S(NAO,NAO),X(NAO,NMO)
      DIMENSION DEN(NAO,NPTS),DENATM(NAT,NPTS),DENTOT(NPTS)
C*******************************************************************
       PI=4.0D+00*ATAN(1.0D+00)
       AUEV = 27.2113957D+00
       GSIGMA = 0.3D+00/AUEV
       GEMIN = EIGS(1) - 0.0D+00*GSIGMA
       GEMAX = EIGS(NMO) + 0.0D+00*GSIGMA
       GDELTA = (GEMAX - GEMIN)/(NPTS - 1)
C
      CALL SEQOPN(NEG0,'ELGDOS','UNKNOWN',.FALSE.,'FORMATTED')
         CALL SEQREW(NEG0)
C
      DO I=1,NAO
      DO J=1,NMO
            X(I,J) = 0.0D+00
         DO L=1,NAO
            X(I,J) = X(I,J) + C(I,J)*S(I,L)*C(L,J)
         ENDDO
      ENDDO
      ENDDO
C
      DO I=1,NAO
        DO IENE=1,NPTS
            DEN(I,IENE) = 0.0D+00
            ENE = GEMIN + (IENE-1)*GDELTA
            DO J=1,NMO
              AAA = (ENE-EIGS(J))**2/(2.0D+00*GSIGMA**2)
              DELT = DEXP(-AAA)/(SQRT(2.0D+00*PI)*GSIGMA)
              DEN(I,IENE) = DEN(I,IENE) + X(I,J)*DELT
            ENDDO
        ENDDO
      ENDDO
C
      DO IENE=1,NPTS
          DO IT=1,NAT
             DENATM(IT,IENE) = 0.0D+00
          ENDDO
        DO I=1,NAO
           IDB=IATM(I)
            DENATM(IDB,IENE) = DENATM(IDB,IENE) + DEN(I,IENE)
        ENDDO
      ENDDO
C
         DO IENE=1,NPTS
           DENTOT(IENE) = 0.0D+00
           DO I=1,NAT
            DENTOT(IENE) = DENTOT(IENE) + DENATM(I,IENE)
           ENDDO
         ENDDO
C
      DO I=1,NAT
         IF(MASWRK) WRITE(NEG0,'('' ATOM: '', I4)') I
         XXX = 0.0D+00
         DO IENE=1,NPTS
           ENE = GEMIN + (IENE-1)*GDELTA
           XXX = XXX + DENATM(I,IENE)*GDELTA
        IF(MASWRK)
     *   WRITE(NEG0,'(F10.4,F14.8)')ENE*AUEV,DENATM(I,IENE)
         ENDDO
           IF(MASWRK) WRITE(NEG0,'(/,10X,F14.8)') XXX
      ENDDO
C
      DENT = 0.0D+00
      IF(MASWRK) WRITE(NEG0,*)'TOTAL DOS'
      DO IENE=1,NPTS
         ENE = GEMIN + (IENE-1)*GDELTA
         DENT = DENT + DENTOT(IENE)*GDELTA
      IF(MASWRK) WRITE(NEG0,'(F10.4,F14.8)')ENE*AUEV,DENTOT(IENE)
      ENDDO
           IF(MASWRK) WRITE(NEG0,'(/,''TOTAL DOS: '',10X,F14.8)') DENT
C
      RETURN
      END
