c23456789
      PROGRAM dosipr
     
      parameter (atoms=1000)
      parameter (orbls=9)
      parameter (maxsta=2000) 
      character orbi,comp,occu,sp
      integer nions, parts, norb
      real nu,wf(1000,9,2000)  !wf(atom,orbital,state)
      real ipr(maxsta),ipryla(maxsta),iprala(maxsta),ener(maxsta)
      integer nks,states,nel,p,no
      real ksen(maxsta),oci(maxsta),dos(maxsta),sum1,const,sig
    
 
      states=1387  !no. of states
      nel=int(real(2454)/2.0)  !no. of electrons
      ef=0.550679 ! grep "CHEMICAL POTENTIAL" cpmd.output
      
      nions=460
      norb= 9  !how many orbitals s+p=5, s+p+d=9
     
      sig=0.04
      const=1.0/(sig*sqrt(2*3.14159265))

      parts=int(real(states)/8.0)
      states=parts*8

      open(10,file='ksener.dat') !from "EIGENVALUES(EV) AND OCCUPATION:" to "CHEMICAL POTENTIAL" in cpmd.output
      open(12,file='projec.dat') !from "WAVEFUNCTIONS IN ATOMIC ORBITAL BASIS" to "WAVEFUNCTIONS IN ATOMIC ORBITAL BASIS" in cpmd.output
      open(14,file='dos-ipr.out')

        !read ksener.dat
      do i=1,states,2
      read(10,*) nks,ksen(i),oci(i),nks,ksen(i+1),oci(i+1)
      ksen(i)=ksen(i)-ef
      ksen(i+1)=ksen(i+1)-ef
      end do

        !read projec.dat
      do p=1,parts
      read(12,*) orbi,nu
      read(12,*) comp,nu
      read(12,*) occu,nu
       do j=1,nions
        read(12,*) no,sp,orbi, (wf(j,1,k),k=8*(p-1)+1,8*p)
         do m=2,norb
          read(12,*) orbi, (wf(j,m,k),k=8*(p-1)+1,8*p)
         enddo
       enddo
      enddo
      write(*,*) "ksener.dat and projec.dat readen"

       ! IPR(i)= [sum(wf(i)**4)]/[sum(wf(i)**2)]**2
      do i=1,1700 !states
       do j=1,nions
        do m=1,norb
        ipryla(i)=ipryla(i)+wf(j,m,i)**4
        iprala(i)=iprala(i)+wf(j,m,i)**2
        enddo
       enddo
       if (iprala(i).eq.0) then
        iprala(i)=1
       endif
       ipr(i)=ipryla(i)/(iprala(i)**2)
 
       !Density Of States (DOS)
       ener(i)= -15.00+i*0.01
       sum1=0
       do j=1,states
        sum1=sum1+const*exp(-0.5*((ksen(j)-ener(i))/sig)**2)
       enddo
       dos(i)=sum1/(2*real(nel))

       write(14,'(4x,35g14.5)') ener(i),dos(i),ksen(i),ipr(i)
      enddo
      !gnuplot
      !p 'dos-ipr.out' u 1:2 w l, '' u 3:4 w i

      write(*,*) "Done!"

      STOP
      END
      !by Julen Larrucea.