subroutine mkdos
      implicit double precision(a-h,o-z)
#include "divcon.dim"
#include "divcon.h"

C     Generates Density of States
C
C     By Arjan van der Vaart

      dimension itmpc(maxsub), itmpb1(maxsub), itmpb2(maxsub)
      logical notb1, notb2

C-RDC       write(iout,'(//" DENSITY OF STATE ANALYSIS:",
C-RDC      $     /" -------------------------")')

C     find the lowest and hightest eigenvalue
      emin = eval(iorbpt(1))
      emax = eval(iorbpt(2)-1)
      do k=2,nsub
         emin = min(emin,eval(iorbpt(k)))
         emax = max(emax,eval(iorbpt(k+1)-1))
      enddo
      if (emin.lt.0.0) then
         emin = int(emin) - 1.0
      else
         emin = int(emin)
      endif
      if (emax.lt.0.0) then
         emax = int(emax)
      else
         emax = int(emax) + 1.0
      endif

      if (stand) then
C-RDC          write(iout,'(/" CONTRACTED ALL DOS PARAMETERS, SINCE ",
C-RDC      $        "A STANDARD CALCULATION",
C-RDC      $        /" WAS PERFORMED. A DENSITY OF STATE ANALYSIS WILL",
C-RDC      $        " NOW BE",
C-RDC      $        /" PERFORMED ON",
C-RDC      $        " ALL ATOMS (1-",I6,") WITH INTERVAL=",F12.3)')
C-RDC      $        natoms, ddos(1)
         ndos = 1
         idos(1) = 1
         idos(2) = natoms
         nc = 1
         itmpc(1) = 1
         nb1 = 0
         nb2 = 0
         call gendos(ndos,nc,itmpc,nb1,itmpb1,nb2,itmpb2,emin,emax)
         return
      endif

      do id=1,ndos
         id2=2*id
         id1=id2-1
         ia1 = idos(id1)
         ia2 = idos(id2)
         nc = 0
         nb1 = 0
         nb2 = 0
         do 100 is=1,nsub
            notb1 = .true.
            notb2 = .true.
            j1 = iatom1(is)
            j2 = iatom1(is+1)-1
            do j=j1,j2
               if (iatoms(j).gt.ia2) goto 100

C     .        Since core atoms have the highest priority,
C     .        we keep on looping over all atoms until a core
C     .        atom is found. We don't have to check for buffer2
C     .        atoms if a buffer1 atom was already found, since
C     .        buffer1 atoms have a higher priority than buffer2
C     .        atoms.
               if ((iatoms(j).ge.ia1).and.(iatoms(j).le.ia2)) then
                  if (iabuff(j).eq.0) then
                     nc = nc + 1
                     itmpc(nc) = is
                     goto 100
                  elseif (notb1.and.(iabuff(j).eq.1)) then
                     nb1 = nb1 + 1
                     itmpb1(nb1) = is
                     notb1 = .false.
                     notb2 = .false.
                  elseif (notb2.and.(iabuff(j).eq.2)) then
                     nb2 = nb2 + 1
                     itmpb2(nb2) = is
                     notb2 = .false.
                  endif
               endif
            enddo
 100     continue

C     .  some subsystems might now be listed twice (or more),
C     .  but for the DOS we only want to access each subsystem
C     .  once (at most).
C     .  so do some sorting now, to get a unique list:
C     .  core has the highest priority, then inner buffer, then
C     .  outer buffer
C     .  if the same subsystem occurs in two or three groups, only
C     .  the entry in the group with highest priority will survive,
C     .  e.g. subsystem 6 is both in core (itmpc) and outer buffer
C     .  (itmpb2), then the entry in the outer buffer will be deleted

C     .  A. Compare outer and inner buffer
         kb1 = 1
         kb2 = 1
 110     if (kb2.le.nb2) then
 120        if (itmpb1(kb1).lt.itmpb2(kb2)) then
               kb1 = kb1 + 1
               if (kb1.gt.nb1) goto 150
               goto 120
            elseif (itmpb1(kb1).gt.itmpb2(kb2)) then
               kb2 = kb2 +1
               goto 110
            else
               do i=kb2,nb2-1
                  itmpb2(i) = itmpb2(i+1)
               enddo
               nb2 = nb2 - 1
               goto 110
            endif
         endif

C     .  B. Compare outer buffer and core
 150     kc = 1
         kb2 = 1
 210     if (kb2.le.nb2) then
 220        if (itmpc(kc).lt.itmpb2(kb2)) then
               kc = kc + 1
               if (kc.gt.nc) goto 250
               goto 220
            elseif (itmpc(kc).gt.itmpb2(kb2)) then
               kb2 = kb2 +1
               goto 210
            else
               do i=kb2,nb2-1
                  itmpb2(i) = itmpb2(i+1)
               enddo
               nb2 = nb2 - 1
               goto 210
            endif
         endif

C     .  C. Compare inner buffer and core
 250     kc = 1
         kb1 = 1
 310     if (kb1.le.nb1) then
 320        if (itmpc(kc).lt.itmpb1(kb1)) then
               kc = kc + 1
               if (kc.gt.nc) goto 350
               goto 320
            elseif (itmpc(kc).gt.itmpb1(kb1)) then
               kb1 = kb1 +1
               goto 310
            else
               do i=kb1,nb1-1
                  itmpb1(i) = itmpb1(i+1)
               enddo
               nb1 = nb1 - 1
               goto 310
            endif
         endif

 350     call gendos(id,nc,itmpc,nb1,itmpb1,nb2,itmpb2,emin,emax)
      enddo

      end

CC---------------------------------------------------------------CC

      subroutine gendos(id,nc,itmpc,nb1,itmpb1,nb2,itmpb2,emin,emax)
      implicit double precision(a-h,o-z)
#include "divcon.dim"
#include "divcon.h"
      parameter (maxh=1000)
      dimension itmpc(maxsub), itmpb1(maxsub), itmpb2(maxsub),
     $     ihistco(maxh), ihistb1o(maxh), ihistb2o(maxh),
     $     ihistcu(maxh), ihistb1u(maxh), ihistb2u(maxh)

      id2 = 2*id
      id1 = id2-1
      nn = (emax-emin)/ddos(id) + 1
      if (nn.gt.maxh) then
         d = (emax-emin)/(maxh-1) + 1.0E-4
C-RDC          write(iout,'(//" THE DENSITY OF STATE FOR ATOMS ",i6," - ",i6,
C-RDC      $        " COULD NOT",
C-RDC      $        /" BE CALCULATED: INTERVAL IS TOO SMALL.",
C-RDC      $        /" INCREASE INTERVAL TO ",F12.3,
C-RDC      $        " OR INCREASE MAXH TO ",i9,
C-RDC      $        /" IN mkdos.F AND RECOMPILE")')
C-RDC      $        idos(id1),idos(id2), d, nn
         return
      endif

C     first for core
      do i=1,maxh
         ihistco(i) = 0
         ihistcu(i) = 0
      enddo
      nhco = 0
      nhcu = 0
      do ic=1,nc
         is = itmpc(ic)
         if (stand) then
C     .     base occupation on the number of electrons
            iocc = (nelec-1)/2
            i1 = iorbpt(is)
            i2 = i1+iocc
            do i=i1,i2
               j = (eval(i)-emin)/ddos(id) + 1
               ihistco(j) = ihistco(j) + 1
            enddo
            nhco = nhco + i2-i1+1
         else
C     .     base occupation on the Fermi energy (since we don't
C     .     know the charge of this subsystem)
            do icr=1,ncores
               i1=isubend(icr-1)+1
               i2=isubend(icr)
               if ((is.ge.i1).and.(is.le.i2)) then
                  ef = efermi(icr)
                  goto 10
               endif
            enddo
 10         i1 = iorbpt(is)
            i2 = iorbpt(is+1)-1
            do i=i1,i2
               if (eval(i).lt.ef) then
                  j = (eval(i)-emin)/ddos(id) + 1
                  ihistco(j) = ihistco(j) + 1
               else
                  ibeg = i
                  goto 20
               endif
            enddo
 20         i2=ibeg-1
            nhco = nhco + i2-i1+1
         endif

         i1 = i2+1
         i2 = iorbpt(is+1)-1
         do i=i1,i2
            j = (eval(i)-emin)/ddos(id) + 1
            ihistcu(j) = ihistcu(j) + 1
         enddo
         nhcu = nhcu + i2-i1+1

      enddo

      nhc = nhcu + nhco

      if (stand) then
C-RDC C     .  we finished the DOS analysis already, so write stuff

C-RDC          write(iout,'(/" OCCUPIED ORBITALS:",
C-RDC      $        //"     ENERGY (eV)",
C-RDC      $        /"     (START OF",
C-RDC      $        /"      INTERVAL)",
C-RDC      $        /"                            N       %",
C-RDC      $        /" --------------------------------------")')

C-RDC          do i=nn,1,-1
C-RDC             if (ihistco(i).ne.0) goto 30
C-RDC          enddo
C-RDC  30      iend=i
C-RDC          y = emin-ddos(id)
C-RDC          tcp = 0.0
C-RDC          do i=1,iend
C-RDC             y = y + ddos(id)
C-RDC             cp = dble(100*ihistco(i))/dble(nhc)
C-RDC             tcp = tcp + cp
C-RDC             write(iout,'(1x,F15.3,4x,i9,2x,f6.2)')
C-RDC      $           y,ihistco(i), cp
C-RDC          enddo
C-RDC          write(iout,'(/"                    ---------  ------  ",
C-RDC      $     /" TOTAL          ",4x,i9,2x,f6.2)')
C-RDC      $        nhco,tcp

C-RDC          write(iout,'(/" UNOCCUPIED ORBITALS:",
C-RDC      $        //"     ENERGY (eV)",
C-RDC      $        /"     (START OF",
C-RDC      $        /"      INTERVAL)",
C-RDC      $        /"                            N       %",
C-RDC      $        /" --------------------------------------")')

         do i=1,nn
            if (ihistcu(i).ne.0) goto 40
         enddo
 40      ibeg=i
         y = emin+dble(ibeg-2)*ddos(id)
         tcp = 0.0
         do i=ibeg,nn
            y = y + ddos(id)
            cp = dble(100*ihistcu(i))/dble(nhc)
            tcp = tcp + cp
C-RDC             write(iout,'(1x,F15.3,4x,i9,2x,f6.2)')
C-RDC      $           y,ihistcu(i), cp
         enddo
C-RDC          write(iout,'(/"                    ---------  ------  ",
C-RDC      $     /" TOTAL          ",4x,i9,2x,f6.2)')
C-RDC      $        nhcu,tcp

C-RDC          write(iout,'(/" ALL ORBITALS:",
C-RDC      $        //"     ENERGY (eV)",
C-RDC      $        /"     (START OF",
C-RDC      $        /"      INTERVAL)",
C-RDC      $        /"                            N       %",
C-RDC      $        /" --------------------------------------")')

         y = emin-ddos(id)
         tcp = 0.0
         do i=1,nn
            y = y + ddos(id)
            cp = dble(100*(ihistcu(i)+ihistco(i)))/dble(nhc)
            tcp = tcp + cp
C-RDC             write(iout,'(1x,F15.3,4x,i9,2x,f6.2)')
C-RDC      $           y,ihistcu(i)+ihistco(i), cp
         enddo
C-RDC          write(iout,'(/"                    ---------  ------  ",
C-RDC      $     /" TOTAL          ",4x,i9,2x,f6.2)')
C-RDC      $        nhc,tcp

      else

C     .  still have to do the DOS of the buffers

         do i=1,maxh
            ihistb1o(i) = ihistco(i)
            ihistb1u(i) = ihistcu(i)
         enddo

C     .  do the inner buffer
         nhb1o = nhco
         nhb1u = nhcu
         do ic=1,nb1
            is = itmpb1(ic)
C     .     base occupation on the Fermi energy (since we don't
C     .     know the charge of this subsystem)
            do icr=1,ncores
               i1=isubend(icr-1)+1
               i2=isubend(icr)
               if ((is.ge.i1).and.(is.le.i2)) then
                  ef = efermi(icr)
                  goto 50
               endif
            enddo
 50         i1 = iorbpt(is)
            i2 = iorbpt(is+1)-1
            do i=i1,i2
               if (eval(i).lt.ef) then
                  j = (eval(i)-emin)/ddos(id) + 1
                  ihistb1o(j) = ihistb1o(j) + 1
               else
                  ibeg = i
                  goto 60
               endif
            enddo
 60         i2=ibeg-1
            nhb1o = nhb1o + i2-i1+1

            i1 = i2+1
            i2 = iorbpt(is+1)-1
            do i=i1,i2
               j = (eval(i)-emin)/ddos(id) + 1
               ihistb1u(j) = ihistb1u(j) + 1
            enddo
            nhb1u = nhb1u + i2-i1+1
         enddo
         nhb1 = nhb1o + nhb1u
         do i=1,maxh
            ihistb2o(i) = ihistb1o(i)
            ihistb2u(i) = ihistb1u(i)
         enddo

C     .  then the outer buffer
         nhb2o = nhb1o
         nhb2u = nhb1u
         do ic=1,nb2
            is = itmpb2(ic)
C     .     base occupation on the Fermi energy (since we don't
C     .     know the charge of this subsystem)
            do icr=1,ncores
               i1=isubend(icr-1)+1
               i2=isubend(icr)
               if ((is.ge.i1).and.(is.le.i2)) then
                  ef = efermi(icr)
                  goto 70
               endif
            enddo
 70         i1 = iorbpt(is)
            i2 = iorbpt(is+1)-1
            do i=i1,i2
               if (eval(i).lt.ef) then
                  j = (eval(i)-emin)/ddos(id) + 1
                  ihistb2o(j) = ihistb2o(j) + 1
               else
                  ibeg = i
                  goto 80
               endif
            enddo
 80         i2=ibeg-1
            nhb2o = nhb2o + i2-i1+1

            i1 = i2+1
            i2 = iorbpt(is+1)-1
            do i=i1,i2
               j = (eval(i)-emin)/ddos(id) + 1
               ihistb2u(j) = ihistb2u(j) + 1
            enddo
            nhb2u = nhb2u + i2-i1+1
         enddo
         nhb2 = nhb2o + nhb2u

C-RDC          write(iout,'(//" THE DENSITY OF STATE FOR ATOMS ",I6," - ",
C-RDC      $        I6," AS CORE ATOMS ",
C-RDC      $        /" ARE OBTAINED FROM THE EIGENVALUES OF SUBSYSTEMS")')
C-RDC      $        idos(id1),idos(id2)
C-RDC          write(iout,'(8(1x,i6))') (itmpc(i),i=1,nc)
C-RDC          write(iout,'(" THE DOS FOR ATOMS ",I6," - ",
C-RDC      $        I6," AS CORE OR INNER BUFFER ATOMS ",
C-RDC      $        /" ARE OBTAINED FROM THE EIGENVALUES OF THE ",
C-RDC      $        "SUBSYSTEMS ABOVE AND SUBSYSTEMS")')
C-RDC      $        idos(id1),idos(id2)
         if (nb1.eq.0) then
C-RDC             write(iout,'("  ---")')
         else
C-RDC             write(iout,'(8(1x,i6))') (itmpb1(i),i=1,nb1)
         endif
C-RDC          write(iout,'(" THE DOS FOR ATOMS ",I6," - ",
C-RDC      $        I6," AS CORE, INNER OR OUTER BUFFER ATOMS ",
C-RDC      $        /" ARE OBTAINED FROM THE EIGENVALUES OF ALL ",
C-RDC      $        "THESE SUBSYSTEM, PLUS",/" SUBSYSTEMS")')
C-RDC      $        idos(id1),idos(id2)
         if (nb2.eq.0) then
C-RDC             write(iout,'("  ---")')
         else
C-RDC             write(iout,'(8(1x,i6))') (itmpb2(i),i=1,nb2)
         endif

C-RDC C     .  and write the Occupied DOS
C-RDC          write(iout,'(/" OCCUPIED ORBITALS:",
C-RDC      $        //"     ENERGY (eV)             CORE     ",
C-RDC      $        "          CORE               CORE",
C-RDC      $        /"     (START OF                                ",
C-RDC      $        "+ BUFF1            + BUFF1",
C-RDC      $        /"      INTERVAL)                           ",
C-RDC      $        "                       + BUFF2",
C-RDC      $        //"                            N       %",
C-RDC      $        "          N       %          N       %",
C-RDC      $        /" -------------------------------------------",
C-RDC      $        "-------------------------------")')

         do i=nn,1,-1
            if (ihistco(i).ne.0) goto 130
         enddo
 130     iendc=i
         do i=nn,1,-1
            if (ihistb1o(i).ne.0) goto 140
         enddo
 140     iendb1=i
         do i=nn,1,-1
            if (ihistb2o(i).ne.0) goto 150
         enddo
 150     iendb2=i
         iend = max(iendc,iendb1,iendb2)

         y = emin-ddos(id)
         tcp = 0.0
         tb1 = 0.0
         tb2 = 0.0
         do i=1,iend
            y = y + ddos(id)
            cp = dble(100*ihistco(i))/dble(nhc)
            tcp = tcp + cp
            b1p = dble(100*ihistb1o(i))/dble(nhb1)
            tb1 = tb1 + b1p
            b2p = dble(100*ihistb2o(i))/dble(nhb2)
            tb2 = tb2 + b2p
C-RDC             write(iout,'(1x,F15.3,2x,3(2x,i9,2x,f6.2))')
C-RDC      $           y,ihistco(i), cp,
C-RDC      $           ihistb1o(i), b1p, ihistb2o(i), b2p
         enddo
C-RDC          write(iout,'(/"                    ---------  ------  ",
C-RDC      $        "---------  ------  ---------  ------",
C-RDC      $     /" TOTAL          ",2x,3(2x,i9,2x,f6.2))')
C-RDC      $        nhco,tcp,nhb1o,tb1,nhb2o,tb2

C-RDC C     .  and write the Unoccupied DOS
C-RDC          write(iout,'(/" UNOCCUPIED ORBITALS:",
C-RDC      $        //"     ENERGY (eV)             CORE     ",
C-RDC      $        "          CORE               CORE",
C-RDC      $        /"     (START OF                                ",
C-RDC      $        "+ BUFF1            + BUFF1",
C-RDC      $        /"      INTERVAL)                           ",
C-RDC      $        "                       + BUFF2",
C-RDC      $        //"                            N       %",
C-RDC      $        "          N       %          N       %",
C-RDC      $        /" -------------------------------------------",
C-RDC      $        "-------------------------------")')

         do i=1,nn
            if (ihistcu(i).ne.0) goto 160
         enddo
 160     ibegc=i
         do i=1,nn
            if (ihistb1u(i).ne.0) goto 170
         enddo
 170     ibegb1=i
         do i=1,nn
            if (ihistb2u(i).ne.0) goto 180
         enddo
 180     ibegb2=i
         ibeg = min(ibegc,ibegb1,ibegb2)

         y = emin+dble(ibeg-2)*ddos(id)
         tcp = 0.0
         tb1 = 0.0
         tb2 = 0.0
         do i=ibeg,nn
            y = y + ddos(id)
            cp = dble(100*ihistcu(i))/dble(nhc)
            tcp = tcp + cp
            b1p = dble(100*ihistb1u(i))/dble(nhb1)
            tb1 = tb1 + b1p
            b2p = dble(100*ihistb2u(i))/dble(nhb2)
            tb2 = tb2 + b2p
C-RDC             write(iout,'(1x,F15.3,2x,3(2x,i9,2x,f6.2))')
C-RDC      $           y,ihistcu(i), cp,
C-RDC      $           ihistb1u(i), b1p, ihistb2u(i), b2p
         enddo
C-RDC          write(iout,'(/"                    ---------  ------  ",
C-RDC      $        "---------  ------  ---------  ------",
C-RDC      $     /" TOTAL          ",2x,3(2x,i9,2x,f6.2))')
C-RDC      $        nhcu,tcp,nhb1u,tb1,nhb2u,tb2

C-RDC C     .  and write the total DOS
C-RDC          write(iout,'(/" ALL ORBITALS:",
C-RDC      $        //"     ENERGY (eV)             CORE     ",
C-RDC      $        "          CORE               CORE",
C-RDC      $        /"     (START OF                                ",
C-RDC      $        "+ BUFF1            + BUFF1",
C-RDC      $        /"      INTERVAL)                           ",
C-RDC      $        "                       + BUFF2",
C-RDC      $        //"                            N       %",
C-RDC      $        "          N       %          N       %",
C-RDC      $        /" -------------------------------------------",
C-RDC      $        "-------------------------------")')

         y = emin-ddos(id)
         tcp = 0.0
         tb1 = 0.0
         tb2 = 0.0
         do i=1,nn
            y = y + ddos(id)
            cp = dble(100*(ihistco(i)+ihistcu(i)))/dble(nhc)
            tcp = tcp + cp
            b1p = dble(100*(ihistb1o(i)+ihistb1u(i)))/dble(nhb1)
            tb1 = tb1 + b1p
            b2p = dble(100*(ihistb2o(i)+ihistb2u(i)))/dble(nhb2)
            tb2 = tb2 + b2p
C-RDC             write(iout,'(1x,F15.3,2x,3(2x,i9,2x,f6.2))')
C-RDC      $           y,ihistco(i)+ihistcu(i), cp,
C-RDC      $           ihistb1o(i)+ihistb1u(i), b1p,
C-RDC      $           ihistb2o(i)+ihistb2u(i), b2p
         enddo
C-RDC          write(iout,'(/"                    ---------  ------  ",
C-RDC      $        "---------  ------  ---------  ------",
C-RDC      $     /" TOTAL          ",2x,3(2x,i9,2x,f6.2))')
C-RDC      $        nhc,tcp,nhb1,tb1,nhb2,tb2

      endif

      end