!! CASTEPDOS - Extract density of states from CASTEP
!! Copyright (C) 2010 G.J.Ackland, University of Edinburgh.
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2 of the License, or
!! (at your option) any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
!!
!! You can contact the author by email on g.j.ackland@ed.ac.uk
!!
!!========================================================================
c  Program to create an interpolated DOS from MP k-point data from CASTEP
c  Inputs required are number of kpoints in the MP grid, the Fermi Energy, 
c  and an array Energy(NKX,NKY,NKZ,NBANDS)
C  The program works by rebuilding the complete k-point set
C  in an array enmat with the appropriate energies (this requires
C  figuring out equivalent point in the reciprocal space and is
C  non-trivial for non-cubic cells.
c  Code by Graeme Ackland, see e.g.#
c# "Embrittlement and the bistable crystal structure of zirconium hydride" G.J.Ackland, Phys.Rev.Letters, 80, 2233 (1998)
c Read PREAMBLE
        PROGRAM CASTEPDOS
        common/symms/ r(3,3,48), d(3,48)
        common/cell/ b0(3,3),g0(3,3),ginv(3,3)
        character*19 ajnk
        character*25 bjnk
        character*30 cjnk
        character(len=80) :: in_file,sym_file,seedname
        character(len=80) :: out_file_1,out_file_2,out_11
        write(*,*)" Enter castep seedname "
        read (*,*)seedname
        in_file=trim(seedname)//".bands"
        open (file=in_file,unit=33,form='formatted',status='old')

        write(*,*)"Getting eigenvalues from ", in_file
        sym_file=trim(seedname)//"-symm.cell"
        open (file=sym_file,unit=35,
     &        form='formatted',status='old',err=99)
        write(*,*)"Symmetry from ", sym_file
        nsymm=1
        goto 101
 99     sym_file=trim(seedname)//".symm"
        open (file=sym_file,unit=34,
     &        form='formatted',status='old',err=100)
        write(*,*)"Symmetry from ",sym_file
        nsymm=2
        goto 101
 100    nsymm=0
        write(*,*)"No symmetry file: assuming all k-points explicit"    
 101    write(*,*)" Enter NKX, NKY, NKZ "

        read(*,*) NKX, NKY, NKZ
        write(*,*) NKX, NKY, NKZ
!  Read in from the .bands file
 111    format(A20,I10)
 112    format(A25,I5)
 113    format(A33,f10.3)
        write(*,*) "You're welcome to castepdos by gjackland@ed.ac.uk"
        read(33,111)ajnk, NKP
c        write(*,*)ajnk, NKP
        read(33,112)bjnk, NSPINS
c        write(*,*)NSPINS
c      if(NSPINS.NE.1)write(*,*)"Sorry magnetic dos not yet supported ", 
c     &   "Workaround by making two .bands files by hand"
        read(33,113)ajnk, ELECT
        NELECT = ELECT
        read(33,112)bjnk, NBANDS
        read(33,113)cjnk, EFERMI 
        read(33,113)cjnk
c        write(*,*)NELECT
         write(*,*)"Nbands " , NBANDS
c        write(*,*)EFERMI 
        if(NSPINS.NE.1) write(*,*)"EFermi=",EFERMI

        out_file_1=trim(seedname)//".dos"
        out_11=trim(seedname)//".11"
        open (file=out_file_1,unit=10,form='formatted',status='unknown')
        open (file=out_11,unit=11,form='formatted',status='unknown')
        if(NSPINS.EQ.2) THEN
         out_file_2=trim(seedname)//".dos2"
         open (file=out_file_2,unit=9,form='formatted',status='unknown')
        endif 

       do j =1,3      
        read(33,*)(B0(I,J),i=1,3) 
       enddo
      call recip(b0,g0,vol)
      call invert(g0,ginv,vol)
        if(nsymm.eq.1)      call readsymm_snap(nsymm)
        if(nsymm.eq.2)      call readsymm(nsymm)
        call dos(NKP,NKX,NKY,NKZ,NBANDS,EFERMI,NSYMM,NSPINS)
        stop
        end program castepdos

        subroutine DOS(NKP,NKX,NKY,NKZ,NBANDS,EFERMI,NSYMM,NSPINS)
        common/symms/ r(3,3,48), d(3,48)
        common/cell/ b0(3,3),g0(3,3),ginv(3,3)
        dimension Enmat(NKX,NKY,NKZ,NBANDS,NSPINS)
        dimension Enmat2(2*NKX+1,2*NKY+1,2*NKZ+1,NBANDS,NSPINS)
        character*10 CHARKP,CHAR

    
c  Read in Energies from Monkhorst Pack
        do  100 KLOOP = 1,NKP
        read(33,*)CHARKP,N,AKX,AKY,AKZ,Weight
        read(33,*)CHAR
c fill both k and -k between 1 and NKx
        IX = NINT( NKX*AKX +(NKX+1)/2.0 )
        IY = NINT( NKY*AKY +(NKY+1)/2.0 )
        IZ = NINT( NKZ*AKZ +(NKZ+1)/2.0 )
        IXB = NINT(-NKX*AKX+(NKX+1)/2.0 )
        IYB = NINT(-NKY*AKY+(NKY+1)/2.0 )
        IZB = NINT(-NKZ*AKZ+(NKZ+1)/2.0 )
c        write(*,*)"Read ",N,IX,IY,IZ,AKX,AKY,AKZ
        do  NB = 1,NBANDS
        read(33,*) Enmat(IX,IY,IZ,NB,1)
        enmat(IXB,IYB,IZB,NB,1)= Enmat(IX,IY,IZ,NB,1)
        enddo
        IF(NSPINS.Eq.2)then
        read(33,*)CHAR
        do  NB = 1,NBANDS
        read(33,*) Enmat(IX,IY,IZ,NB,2)
        enmat(IXB,IYB,IZB,NB,2)= Enmat(IX,IY,IZ,NB,2)
        enddo
        endif
c  apply symmetry elements
        if(nsymm.eq.0)goto 100
        do 50 I = 1,nsymm
c Apply rotation
         BKX =  r(1,1,i)*AKX+r(1,2,i)*AKY+r(1,3,i)*AKZ
         BKY =  r(2,1,i)*AKX+r(2,2,i)*AKY+r(2,3,i)*AKZ
         BKZ =  r(3,1,i)*AKX+r(3,2,i)*AKY+r(3,3,i)*AKZ

c  Only allow elements which conserve k  
          BKX1=g0(1,1)*BKX + g0(1,2)*BKY + g0(1,3)*BKZ
          BKY1=g0(2,1)*BKX + g0(2,2)*BKY + g0(2,3)*BKZ
          BKZ1=g0(3,1)*BKX + g0(3,2)*BKY + g0(3,3)*BKZ

          AKX1=g0(1,1)*AKX + g0(1,2)*AKY + g0(1,3)*AKZ
          AKY1=g0(2,1)*AKX + g0(2,2)*AKY + g0(2,3)*AKZ
          AKZ1=g0(3,1)*AKX + g0(3,2)*AKY + g0(3,3)*AKZ

c  Put it back in first zone
         if(BKX.ge.0.5)BKX=BKX-1.0  
         if(BKY.ge.0.5)BKY=BKY-1.0  
         if(BKZ.ge.0.5)BKZ=BKZ-1.0  
         if(BKX.Lt.-0.5)BKX=BKX+1.0  
         if(BKY.Lt.-0.5)BKY=BKY+1.0  
         if(BKZ.LT.-0.5)BKZ=BKZ+1.0  

        gridx =  NKX*BKX +(NKX+1)/2.0
        gridy =  NKY*BKY +(NKY+1)/2.0 
        gridz =  NKZ*BKZ +(NKZ+1)/2.0 
        IXR = NINT(gridx)
        IYR = NINT(gridy )
        IZR = NINT(gridz )
        IXB = NINT( -NKX*BKX +(NKX+1)/2.0 )
        IYB = NINT( -NKY*BKY +(NKY+1)/2.0 )
        IZB = NINT( -NKZ*BKZ +(NKZ+1)/2.0 )
c        write(*,*)IXR,IYR,IZR,gridx,gridy,gridz,BKX,BKY,BKZ 
c  Exclude anything which doesn't map to a grid point
        if(abs(IXR-gridx).gt.0.05) goto 50
        if(abs(IYR-gridy).gt.0.05) goto 50
        if(abs(IZR-gridz).gt.0.05) goto 50
         do  NB = 1,NBANDS
         do  ISP = 1,NSPINS
          enmat(IXR,IYR,IZR,NB,ISP)= Enmat(IX,IY,IZ,NB,ISP)
          enmat(IXB,IYB,IZB,NB,ISP)= Enmat(IX,IY,IZ,NB,ISP)
         enddo
         enddo
 654    format(10I4,f10.4)
 655    format(I4,6F9.3)

  50    continue 
 100    continue


c        do  101 K = 1,NKZ
c        do  101 J = 1,NKY
c        do  101 I = 1,NKX
c        do  101 ISP = 1,NSPINS
c           write(11,211)i,j,k,(Enmat(I,J,K,M,ISP),M=1,6)
c 211      format(3i3,6f10.4)
c 101    continue  
c           write(11,*)"Expanded array"
        do 600 ISP = 1,NSPINS
c  Create double sized array
        do 102 IX=1,NKX
         JX = 2*IX
         JXP = 2*IX+1
         IXP = IX+1
         IXM = IX-1
         IXP2 = IX+2
          IF(IXM.LT.1)IXM=IXM+NKX
          IF(IXP.GT.NKX)IXP=IXP-NKX
          IF(IXP2.GT.NKX)IXP2=IXP2-NKX
        do 102 IY=1,NKY
         JY = 2*IY
         JYP = 2*IY+1
         IYP = IY+1
         IYM=IY-1
         IYP2=IY+2
          IF(IYM.LT.1)IYM=IYM+NKY
          IF(IYP.GT.NKY)IYP=IYP-NKY
          IF(IYP2.GT.NKY)IYP2=IYP2-NKY
        do 102 IZ=1,NKZ
         JZ = 2*IZ
         JZP = 2*IZ+1
         IZP = IZ + 1
         IZP2 = IZ + 2
         IZM = IZ - 1
          IF(IZM.LT.1)IZM=IZM+NKZ
          IF(IZP.GT.NKZ)IZP=IZP-NKZ
          IF(IZP2.GT.NKZ)IZP2=IZP2-NKZ
       do 102 IB=1,NBANDS
         Enmat2(JX,JY,JZ,IB,ISP)= Enmat(IX,IY,IZ,IB,ISP)
         Enmat2(JXP,JY,JZ,IB,ISP)= 
     &       0.6d0*(Enmat(IX,IY,IZ,IB,ISP)+Enmat(IXP,IY,IZ,IB,ISP))
     &      -0.1d0*(Enmat(IXM,IY,IZ,IB,ISP)+Enmat(IXP2,IY,IZ,IB,ISP))


         Enmat2(JX,JYP,JZ,IB,ISP)= 
     &       0.6d0*(Enmat(IX,IY,IZ,IB,ISP)+Enmat(IX,IYP,IZ,IB,ISP))
     &      -0.1d0*(Enmat(IX,IYM,IZ,IB,ISP)+Enmat(IX,IYP2,IZ,IB,ISP))
 
         Enmat2(JX,JY,JZP,IB,ISP)= 
     &        0.6d0*(Enmat(IX,IY,IZ,IB,ISP)+Enmat(IX,IY,IZP,IB,ISP))
     &      -0.1d0*(Enmat(IX,IY,IZM,IB,ISP)+Enmat(IX,IY,IZP2,IB,ISP))

         Enmat2(JXP,JYP,JZ,IB,ISP)= 
     &       0.30d0*(Enmat(IX,IY,IZ,IB,ISP)+Enmat(IXP,IY,IZ,IB,ISP)+
     &               Enmat(IX,IYP,IZ,IB,ISP)+Enmat(IXP,IYP,IZ,IB,ISP))
     &      -0.025d0*(Enmat(IXM,IY,IZ,IB,ISP)+Enmat(IXP2,IY,IZ,IB,ISP)+
     &               Enmat(IX,IYM,IZ,IB,ISP)+Enmat(IX,IYP2,IZ,IB,ISP)+
     &               Enmat(IXM,IYP,IZ,IB,ISP)+Enmat(IXP2,IYP,IZ,IB,ISP)+
     &               Enmat(IXP,IYM,IZ,IB,ISP)+Enmat(IXP,IYP2,IZ,IB,ISP))

         Enmat2(JX,JYP,JZP,IB,ISP)= 
     &       0.30d0*(Enmat(IX,IY,IZ,IB,ISP)+Enmat(IX,IYP,IZ,IB,ISP) +
     &               Enmat(IX,IY,IZP,IB,ISP)+Enmat(IX,IYP,IZP,IB,ISP))
     &      -0.025d0*(Enmat(IX,IY,IZM,IB,ISP)+Enmat(IX,IY,IZP2,IB,ISP)+
     &               Enmat(IX,IYM,IZ,IB,ISP)+Enmat(IX,IYP2,IZ,IB,ISP)+
     &               Enmat(IX,IYP,IZM,IB,ISP)+Enmat(IX,IYP,IZP2,IB,ISP)+
     &               Enmat(IX,IYM,IZP,IB,ISP)+Enmat(IX,IYP2,IZP,IB,ISP))


         Enmat2(JXP,JY,JZP,IB,ISP)= 
     &        0.30d0*(Enmat(IX,IY,IZ,IB,ISP)+Enmat(IX,IY,IZP,IB,ISP) +
     &                Enmat(IXP,IY,IZ,IB,ISP)+Enmat(IXP,IY,IZP,IB,ISP))
     &      -0.025d0*(Enmat(IX,IY,IZM,IB,ISP)+Enmat(IX,IY,IZP2,IB,ISP)+
     &               Enmat(IXM,IY,IZ,IB,ISP)+Enmat(IXP2,IY,IZ,IB,ISP)+
     &               Enmat(IXP,IY,IZM,IB,ISP)+Enmat(IXP,IY,IZP2,IB,ISP)+
     &               Enmat(IXM,IY,IZP,IB,ISP)+Enmat(IXP2,IY,IZP,IB,ISP))

         Enmat2(JXP,JYP,JZP,IB,ISP)= 0.125*( Enmat(IX,IY,IZ,IB,ISP) +
     &                                Enmat(IXP,IYP,IZ,IB,ISP) +
     &                                Enmat(IXP,IY,IZP,IB,ISP) +
     &                                Enmat(IX,IYP,IZP,IB,ISP) +
     &                                Enmat(IXP,IY,IZ,IB,ISP) +
     &                                Enmat(IX,IYP,IZ,IB,ISP) +
     &                                Enmat(IX,IY,IZP,IB,ISP) +
     &                                Enmat(IXP,IYP,IZP,IB,ISP) )
 102    enddo
  



c  BZ-X Faces assuming quadratic from previous two real points
        do  IY=2,2*NKY
        do  IZ=2,2*NKZ
        do  IB=1,NBANDS
         Enmat2(1,IY,IZ,IB,ISP) = 0.2d0 * (
     &             Enmat2(2,IY,IZ,IB,ISP)*6d0-Enmat2(3,IY,IZ,IB,ISP))
c        Enmat2(2*NKX+1,IY,IZ,IB,ISP) = 0.2d0 * (
c    &              Enmat2(2*NKX,IY,IZ,IB,ISP)*6d0-Enmat2(2*NKX-1,IY,IZ,IB,ISP))
        enddo
        enddo
        enddo


c  BZ-Y Faces assuming quadratic from previous two points
        do IX=2,2*NKX
        do IZ=2,2*NKZ
        do IB=1,NBANDS
         Enmat2(IX,1,IZ,IB,ISP) = 0.2d0 * (
     &              Enmat2(IX,2,IZ,IB,ISP)*6d0-Enmat2(IX,3,IZ,IB,ISP))
c         Enmat2(IX,2*NKY+1,IZ,IB,ISP) = 0.2d0 * (
c     &              Enmat2(IX,2*NKY,IZ,IB,ISP)*6d0-Enmat2(IX,2*NKY-1,IZ,IB,ISP))
        enddo
        enddo
        enddo

      
c  BZ-Z Faces assuming quadratic  quadratic along 110
        do IX=2,2*NKX
        do IY=2,2*NKY
        do IB=1,NBANDS
         Enmat2(IX,IY,1,IB,ISP) = 0.2d0 * (
     &              Enmat2(IX,IY,2,IB,ISP)*6d0-Enmat2(IX,IY,3,IB,ISP))
c         Enmat2(IX,IY,2*NKZ+1,IB,ISP) = 0.2d0 * (
c     &              Enmat2(IX,IY,2*NKZ,IB,ISP)*6d0-Enmat2(IX,IY,2*NKZ-1,IB,ISP))
        enddo
        enddo
        enddo


c  BZ-X EDGE assuming quadratic along 011
        do IX=2,2*NKX
        do  IB=1,NBANDS
         Enmat2(IX,1,1,IB,ISP)             = 0.5d0 * (
     &     3.0* Enmat2(IX,2,2,IB,ISP)- Enmat2(IX,3,3,IB,ISP))
         Enmat2(IX,2*NKY+1,1,IB,ISP)       = 0.5d0 * (
     &     3.0* Enmat2(IX,2*NKY,2,IB,ISP) - Enmat2(IX,2*NKY-1,3,IB,ISP))
         Enmat2(IX,1,2*NKZ+1,IB,ISP)       = 0.5d0 * (
     &     3.0* Enmat2(IX,2,2*NKZ,IB,ISP) - Enmat2(IX,3,2*NKZ-1,IB,ISP))
c         Enmat2(IX,2*NKY+1,2*NKZ+1,IB,ISP) = 0.5d0 * (
c     &  3.0* Enmat2(IX,2*NKY,2*NKZ,IB,ISP) - Enmat2(IX,2*NKY-1,2*NKZ-1,IB,ISP))
        enddo
        enddo


c  BZ-y EDGE assuming  quadratic along 101
        do IY=2,2*NKY
        do  IB=1,NBANDS
         Enmat2(1,IY,1,IB,ISP)             = 0.5d0 * (
     &    3.0*  Enmat2(2,IY,2,IB,ISP) - Enmat2(3,IY,3,IB,ISP))
         Enmat2(1,IY,2*NKZ+1,IB,ISP)       = 0.5d0 * (
     &    3.0*  Enmat2(2,IY,2*NKZ,IB,ISP) - Enmat2(3,IY,2*NKZ-1,IB,ISP))
         Enmat2(2*NKX+1,IY,1,IB,ISP)       = 0.5d0 * (
     &    3.0*  Enmat2(2*NKX,IY,2,IB,ISP) - Enmat2(2*NKX-1,IY,3,IB,ISP))
c         Enmat2(2*NKX+1,IY,2*NKZ+1,IB,ISP) = 0.5d0 * (
c     &  3.0* Enma t2(2*NKX,IY,2*NKZ,IB,ISP) - Enmat2(2*NKX-1,IY,2*NKZ-1,IB,ISP))
        enddo
        enddo

c  BZ-Z EDGE assuming average of adjacent
        do IZ=2,2*NKZ
        do  IB=1,NBANDS
         Enmat2(1,1,IZ,IB,ISP)           =  0.5d0 * (
     &      3.0*Enmat2(2,2,IZ,IB,ISP) - Enmat2(3,3,IZ,IB,ISP))
         Enmat2(1,2*NKY+1,IZ,IB,ISP)       = 0.5d0 * (
     &      3.0*Enmat2(2,2*NKY,IZ,IB,ISP) - Enmat2(3,2*NKY-1,IZ,IB,ISP))
         Enmat2(2*NKX+1,1,IZ,IB,ISP)       = 0.5d0 * (
     &      3.0*Enmat2(2*NKX,2,IZ,IB,ISP) - Enmat2(2*NKX-1,3,IZ,IB,ISP))
c         Enmat2(2*NKX+1,2*NKY+1,IZ,IB,ISP) = 0.5d0 * (
c     &   3.0*Enmat2(2*NKX ,2*NKY,IZ,IB,ISP) - Enmat2(2*NKX-1,2*NKY-1,IZ,IB,ISP))
        enddo
        enddo

c  Corners
        do  IB=1,NBANDS
        Enmat2(1,1,1,IB,ISP) = 0.4* ( Enmat2(2,1,1,IB,ISP) 
     &+ Enmat2(1,2,1,IB,ISP) +  Enmat2(1,1,2,IB,ISP))     
     &  -0.06666666666666*(    Enmat2(2,2,1,IB,ISP) 
     &  + Enmat2(2,1,2,IB,ISP) +  Enmat2(1,2,2,IB,ISP))

        Enmat2(2*NKX+1,1,1,IB,ISP) = 0.4* (
     &        Enmat2(2*NKX,1,1,IB,ISP) + Enmat2(2*NKX+1,2,1,IB,ISP) 
     &                              +  Enmat2(2*NKX+1,1,2,IB,ISP))       
     &  -0.06666666666666*(  Enmat2(2*NKX,2,1,IB,ISP) 
     &     + Enmat2(2*NKX+1,2,2,IB,ISP)  +  Enmat2(2*NKX,1,2,IB,ISP)) 


        Enmat2(1,2*NKY+1,1,IB,ISP) = 0.4* (
     &    Enmat2(2,2*NKY+1,1,IB,ISP) + Enmat2(1,2*NKY,1,IB,ISP)
     &                              +  Enmat2(1,2*NKY+1,2,IB,ISP))      
     &  -0.06666666666666*(  Enmat2(2,2*NKY,1,IB,ISP)
     &+ Enmat2(1,2*NKY,2,IB,ISP)     +  Enmat2(2,2*NKY+1,2,IB,ISP))  


        Enmat2(1,1,2*NKZ+1,IB,ISP) = 0.4 * (
     &    Enmat2(2,1,2*NKZ+1,IB,ISP) + Enmat2(1,1,2*NKZ,IB,ISP)
     &                              +  Enmat2(1,2,2*NKZ+1,IB,ISP))      
     &  -0.06666666666666*(  Enmat2(2,2,2*NKZ,IB,ISP)+
     &  Enmat2(2,1,2*NKZ,IB,ISP) +  Enmat2(1,2,2*NKZ,IB,ISP))


         Enmat2(1,2*NKY+1,2*NKZ+1,IB,ISP) = 
     &  0.4*( Enmat2(2,2*NKY+1,2*NKZ+1,IB,ISP) 
     &  +Enmat2(1,2*NKY,2*NKZ+1,IB,ISP)+Enmat2(1,2*NKY+1,2*NKZ,IB,ISP)) 
     &  -0.06666666666666*(  Enmat2(2,2*NKY,2*NKZ+1,IB,ISP)+
     &  Enmat2(2,2*NKY+1,2*NKZ,IB,ISP) +  Enmat2(1,2*NKY,2*NKZ,IB,ISP))


         Enmat2(2*NKX+1,1,2*NKZ+1,IB,ISP) = 0.4* (
     &    Enmat2(2*NKX,1,2*NKZ+1,IB,ISP) 
     & +Enmat2(2*NKX+1,2,2*NKZ+1,IB,ISP)+Enmat2(2*NKX+1,1,2*NKZ,IB,ISP))     
     &  -0.06666666666666*(  Enmat2(2*NKX,2,2*NKZ+1,IB,ISP) +
     &   Enmat2(2*NKX,1,2*NKZ,IB,ISP) + Enmat2(2*NKX+1,2,2*NKZ,IB,ISP))

         Enmat2(2*NKX+1,2*NKY+1,1,IB,ISP) = 0.4* (
     &   Enmat2(2*NKX,2*NKY+1,1,IB,ISP) + Enmat2(2*NKX+1,2*NKY,1,IB,ISP) 
     &                            +  Enmat2(2*NKX+1,2*NKY+1,2,IB,ISP))     
     &  -0.06666666666666*(    Enmat2(2*NKX,2*NKY,1,IB,ISP)+
     &  Enmat2(2*NKX,2*NKY+1,2,IB,ISP) + Enmat2(2*NKX+1,2*NKY,2,IB,ISP))

         Enmat2(2*NKX+1,2*NKY+1,2*NKZ+1,IB,ISP) = 0.4* (
     &   Enmat2(2*NKX,2*NKY+1,2*NKZ+1,IB,ISP)
     & + Enmat2(2*NKX+1,2*NKY,2*NKZ+1,IB,ISP) 
     & + Enmat2(2*NKX+1,2*NKY+1,2*NKZ,IB,ISP))     
     &  -0.06666666666666*(   Enmat2(2*NKX,2*NKY,2*NKZ+1,IB,ISP) +
     &  Enmat2(2*NKX,2*NKY+1,2*NKZ,IB,ISP) + 
     &  Enmat2(2*NKX+1,2*NKY,2*NKZ,IB,ISP))


c     write(*,*) Enmat2(1,1,1,IB,ISP),Enmat2(2,1,1,IB,ISP), Enmat2(1,2,1,IB,ISP),   Enmat2(1,1,2,IB,ISP)   
        enddo

       call makedos(NKP,NKX,NKY,NKZ,NBANDS,NSPINS,EFERMI,Enmat2,ISP)
 600  enddo
      return
      end subroutine dos


      subroutine makedos
     &            (NKP,NKX,NKY,NKZ,NBANDS,NSPINS,EFERMI,Enmat2,ISP)
c  Set histogram boundaries
        parameter (NHIST=2000)
        dimension Enmat2(2*NKX+1,2*NKY+1,2*NKZ+1,NBANDS,NSPINS)
        dimension XHIST(NHIST)
        dimension GHIST(NHIST)
        dimension UHIST(NHIST)
       
        do I = 1,NHIST
         XHIST(I)=0.0
         GHIST(I)=0.0
         UHIST(I)=0.0
        enddo
        Enmax = Enmat2(1,1,1,1,ISP)
        Enmin = Enmat2(1,1,1,1,ISP)

        do  101 K = 1,2*NKZ+1
        do  101 J = 1,2*NKY+1
        do  101 I = 1,2*NKX+1
        do  N = 1,NBANDS
          if (Enmat2(I,J,K,N,ISP).GT.Enmax) Enmax = Enmat2(I,J,K,N,ISP)
          if (Enmat2(I,J,K,N,ISP).LT.Enmin) Enmin = Enmat2(I,J,K,N,ISP)
        enddo
          !!write(11,211)i,j,k,(Enmat2(I,J,K,M,ISP),M=1,6)
 211      format(3i3,6f10.4)
 101    continue  




c  Add 10%

        Enrange = Enmax-Enmin
        Enmax = Enmax + 0.1*Enrange
        Enmin = Enmin - 0.1*Enrange
        Enrange = Enmax-Enmin
   
c  expand k-point set by interpolating between each calculated point
c  Can make interpolation as fine as you like by changing IFINE
        ifine = 100

c loop over all kpoints and bands  Fill in between extreme k-points
        do 200 IK=1,2*nkx
        do 200 JK=1,2*nky
c        do 200 K=1,(nkz-1)*ifine
        do 200 K=1,(2*nkz) 
        do 200 IKF=1,IFINE
        do 200 JKF=1,IFINE

        I = (IK-1)*IFINE+IKF
        J = (JK-1)*IFINE+JKF
c Find surrounding kpoints
        ilow = ik
        ihigh = ik+1
        jlow = jk
        jhigh = jk+1
        klow = k
        khigh = k+1
c  Relative position
  
        t = float(IKF)/float(IFINE)
        t1 = 1-t
        u = float(JKF)/float(IFINE)
        u1 = 1-u
         v = float(klow-1)/float(NKZ-1)
         v1 = 1-v
        
c BZ boundaries: assume bands go flat at the boundary in a quadratic way.

       if(ilow .eq. 0) ilow = nkx
        if(jlow .eq. 0) jlow = nky
        if(klow .eq. 0) klow = nkz
        if(ihigh.gt.nkx+1)  ihigh = 1
        if(jhigh .gt.nky+1) jhigh = 1
        if(khigh .gt.nkz+1) khigh = 1

c  Evaluate energies by simple interpolation between the 8 real k-points
c  surrounding the interpolated one.
C  Need to incorporate a more sophisticated scheme here to get boundaries right.
        do 201 iban=1,nbands

        y1 = enmat2(ilow,jlow,klow,iban,isp)
        y2 = enmat2(ihigh,jlow,klow,iban,isp)
        y3 = enmat2(ihigh,jhigh,klow,iban,isp)
        y4 = enmat2(ilow,jhigh,klow,iban,isp)
        y5 = enmat2(ilow,jlow,khigh,iban,isp)
        y6 = enmat2(ihigh,jlow,khigh,iban,isp)
        y7 = enmat2(ihigh,jhigh,khigh,iban,isp)
        y8 = enmat2(ilow,jhigh,khigh,iban,isp)

        efine    =  t1*u1*v1*y1 + t*u1*v1*y2 +
     &         t*u*v1*y3 + t1*u*v1*y4  + t1*u1*v*y5 +
     &         t*u1*v*y6 + t*u*v*y7    + t1*u*v*y8 

cc  End points of line along z

        efinev  =  t1*u1*y5 + t*u1*y6 + t*u*y7 + t1*u*y8 
        efinev1 =  t1*u1*y1 + t*u1*y2 + t*u*y3 + t1*u*y4
c  Integrate along z-line
        iev  = NINT(NHIST*(efinev-enmin)/(enmax-enmin))
        iev1 = NINT(NHIST*(efinev1-enmin)/(enmax-enmin))
        if(iev.eq.iev1) then 
            xhist(iev) = xhist(iev) + 1.0
        elseif (iev.lt.iev1) then
            wdos = 1.0/(iev1-iev)
              xhist(iev) = xhist(iev) - wdos/2.0
              xhist(iev1) = xhist(iev1)- wdos/2.0
            do ixh = iev,iev1
             xhist(ixh) = xhist(ixh)+ wdos
            enddo
        elseif (iev.gt.iev1) then
            wdos = 1.0/(iev-iev1)
              xhist(iev) = xhist(iev) - wdos/2.0
              xhist(iev1) = xhist(iev1)- wdos/2.0
            do ixh = iev1,iev
             xhist(ixh) = xhist(ixh)+ wdos
            enddo
        endif
c  Gaussian Smearing thing.
c  Evaluate point in the histogram
        efrac = (efine-enmin)/(enmax-enmin)
        IEFINE = NHIST*efrac
c       IF(IEFINE.GT.0.AND.IEFINE.LT.NHIST) 
c     &        GHIST(IEFINE)=GHIST(IEFINE)+1.0
 201    continue
 200    continue
        write(*,*)"Completed integration"
c  How many electrons?
         ESUM=0.0      
        do i=1,NHIST
         ESUM=ESUM + XHIST(I)
         GHIST(I) = ESUM
        enddo
         Scale = NBANDS/ESUM
        do I =1,2*NKX+1
        do J =1,2*NKY+1
        do K =1,2*NKZ+1
        do IB =1,NBANDS
        y1 = enmat2(i,j,k,ib,isp)
        efrac = (y1-enmin)/(enmax-enmin)
        IEFINE = NHIST*efrac
c  Uncomment for a Smeared Dos
c       IF(IEFINE.GT.0.AND.IEFINE.LT.NHIST) THEN 
c             UHIST(IEFINE-2)=UHIST(IEFINE)+1.0
c             UHIST(IEFINE-1)=UHIST(IEFINE)+4.0
c             UHIST(IEFINE)=UHIST(IEFINE)+6.0
c             UHIST(IEFINE+1)=UHIST(IEFINE)+4.0
c             UHIST(IEFINE+2)=UHIST(IEFINE)+1.0
c       ENDIF
       enddo
       enddo
       enddo
       enddo
 ! Output in eV

        NUMKP = 8*NKX*NKY*NKZ
        ! efermi = 0  for debug.
        ehartree = 27.21138342902473
c        write(*,*) ndos
        do i=1,NHIST
       write((11-ISP),*) ((I*enrange)/NHIST + enmin-efermi)*ehartree,
     &                  XHIST(I)*SCALE*NHIST/(enrange*ehartree)
     &                  ,GHIST(I)*SCALE*NHIST/(enrange*ehartree)
c     &                   UHIST(I)/IFINE**2
        enddo 
c        stop " Created DOS  Ef=0 Units eV"
         return
        end
 
      subroutine readsymm(nsymm)
      common/symms/ r(3,3,48), d(3,48)
        common/cell/ b0(3,3),g0(3,3),ginv(3,3)
      dimension rin(3,3),temp(3,3)
      character*60 jnk
      character*20 bracket
c       write(*,*) "In readsymm"
c       write(*,12)b0
c       write(*,12)g0
c       write(*,12)ginv
 12     format(3f8.4,/,3f8.4,/,3f8.4)
      ns = 48
      nsymm=0
      do i=1,ns  
      read(34,*,end=100)jnk      
      read(34,*,end=100)bracket, rin(1,1),rin(1,2),rin(1,3)
      read(34,*,end=100)bracket, rin(2,1),rin(2,2),rin(2,3)
      read(34,*,end=100)bracket, rin(3,1),rin(3,2),rin(3,3)
      read(34,*,end=100)jnk      
      read(34,*,end=100)bracket, d(1,i),d(2,i),d(3,i)
      nsymm=nsymm+1
      read(34,*,end=100)jnk      
      read(34,*,end=100)jnk      
c      write(*,12) rin
       call matrix_multiply(3,rin,g0,temp)
       call matrix_multiply(3,ginv,temp,rin)
c      write(*,*)i, "becomes"
c      write(*,12) rin
        do j=1,3
        do k=1,3
          r(j,k,i)= rin(j,k) 
        enddo
        enddo
      enddo
 100  continue
      write(*,*)"read ",nsymm," symmetry elements"  
      
      return
      end




c  Reads symmetry elements and KPOINTS from symmetry_snap output
      subroutine readsymm_snap(nsymm)
      common/symms/ r(3,3,48), d(3,48)
        common/cell/ b0(3,3),g0(3,3),ginv(3,3)
      dimension rin(3,3),temp(3,3)
      character*30 jnk,jnk2,jnk3
      character*60 symmblock
      character*20 bracket
c       write(*,*) "In readsymm"
c       write(*,12)b0
c       write(*,12)g0
c       write(*,12)ginv
 12     format(3f8.4,/,3f8.4,/,3f8.4)
      ns = 48
      nsymm=0
      symmblock="%BLOCK symmetry_ops"  
 10      read(35,*,end=100)jnk,jnk2      
      if(jnk2.NE."symmetry_ops") goto 10
      do i=1,ns  
      read(35,*,end=100)jnk,jnk2
c      write(*,*)jnk,jnk2
      if(jnk2.EQ."symmetry_ops") goto 100
      read(35,*,end=100)rin(1,1),rin(1,2),rin(1,3)
      read(35,*,end=100)rin(2,1),rin(2,2),rin(2,3)
      read(35,*,end=100)rin(3,1),rin(3,2),rin(3,3)
      read(35,*,end=100)d(1,i),d(2,i),d(3,i)
      nsymm=nsymm+1
       call matrix_multiply(3,rin,g0,temp)
       call matrix_multiply(3,ginv,temp,rin)
c      write(*,*)i, "becomes"
c      write(*,12) rin
        do j=1,3
        do k=1,3
          r(j,k,i)= rin(j,k) 
        enddo
        enddo
      enddo
 100  continue
      write(*,*)"read ",nsymm," symmetry elements"  
      
      return
      end








      subroutine recip(a,b,vol)
c Compute unit cell volume and reciprocal lattice primitive translations,
c the latter  WITHOUT the factor of 2*Pi which is sometimes used.
c i.e. (b1 b2 b3) = inverse(transpose (a1 a2 a3)  )  )

      implicit none
      real a,b
      dimension  a(3,3),b(3,3)
      real t1,t2,t3,det,den,vol

c     Compute determinant of (a1 a2 a3)
      t1 = a(2,2) * a(3,3) - a(2,3) * a(3,2)
      t2 = a(2,3) * a(3,1) - a(2,1) * a(3,3)
      t3 = a(2,1) * a(3,2) - a(2,2) * a(3,1)
      det = a(1,1) * t1 + a(1,2) * t2 + a(1,3) * t3
      vol = abs(det)
      den = 1.d0/det

      b(1,1) = t1 * den
      b(1,2) = t2 * den
      b(1,3) = t3 * den
      b(2,1) = (a(1,3)*a(3,2)-a(1,2)*a(3,3))*den
      b(2,2) = (a(1,1)*a(3,3)-a(1,3)*a(3,1))*den
      b(2,3) = (a(1,2)*a(3,1)-a(1,1)*a(3,2))*den
      b(3,1) = (a(1,2)*a(2,3)-a(1,3)*a(2,2))*den
      b(3,2) = (a(1,3)*a(2,1)-a(1,1)*a(2,3))*den
      b(3,3) = (a(1,1)*a(2,2)-a(1,2)*a(2,1))*den

      return
      end

  !--------------------------------------------------------------
  !
  !  subroutine matrix_invert
  !
  !  inverts a martix
  !
  !--------------------------------------------------------------
      subroutine invert(a,b,deta)

      integer i, j
      real a,b
      dimension a(3,3),b(3,3)
      real  deta

    !! determinant - in three stages
      deta=a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2))
      deta=deta + a(2,1)*(a(3,2)*a(1,3)-a(1,2)*a(3,3)) 
      deta=deta + a(3,1)*(a(1,2)*a(2,3)-a(1,3)*a(2,2))
      b(1,1)=a(2,2)*a(3,3)-a(2,3)*a(3,2)
      b(2,1)=a(3,1)*a(2,3)-a(2,1)*a(3,3)
      b(3,1)=a(2,1)*a(3,2)-a(3,1)*a(2,2)
      b(1,2)=a(3,2)*a(1,3)-a(1,2)*a(3,3)
      b(2,2)=a(1,1)*a(3,3)-a(1,3)*a(3,1)
      b(3,2)=a(3,1)*a(1,2)-a(1,1)*a(3,2)
      b(1,3)=a(1,2)*a(2,3)-a(1,3)*a(2,2)
      b(2,3)=a(2,1)*a(1,3)-a(2,3)*a(1,1)
      b(3,3)=a(1,1)*a(2,2)-a(1,2)*a(2,1)

    !! check for determinant of zero
      if(abs(deta).eq.0.d0)then
       stop "Error in matrix_invert: STOP."
      end if
    
    !! normalise
      do i=1,3
      do j=1,3
       b(i,j)=b(i,j)/deta
      enddo
      enddo
      return
      end subroutine invert

       subroutine matrix_multiply(n,a,b,c)
      real a,b,c
      integer  i, j, k
      integer  n                                 !< extent of matrix
      dimension a(n,n), b(n,n) !< input matrices
      dimension c(n,n)         !< resultant matrix
    
      do i=1,n
       do j=1,n
          c(i,j)=0.0d0
       end do
       do k=1,n
          do j=1,n
             c(i,j)=c(i,j)+a(i,k)*b(k,j)
          end do
       end do
      end do
      return

      end subroutine matrix_multiply

