!! Calculates electronic free energy and entropy from a VASP EIGENVAL file
!! cat EIGENVAL | Felect (output T, F, U, -TS in units of eV/atom)
!! Option to set EF by hand
!! Note: requires dense k-point mesh for accurate result - check convergence!
!! Note: mu=EF for T<100K

module MyFermi
contains
function Fermi(kT,mu,E)
    real*16 kT,mu,E(:),Fermi(size(E))
    Fermi = 1./(1.+exp((E-mu)/kT))
end function Fermi
function FermiP(kT,mu,E)
    real*16 kT,mu,E(:),FermiP(size(E))            
    FermiP = exp((E-mu)/kT)/kT/(1.+exp((E-mu)/kT))**2
end function FermiP
end module MyFermi

use MyFermi
implicit none
real*16 T,T0,T1,dT,Stot,kB,Etot,U0,V0,V1
real*16 mu,Eb1,Eb2,wk,eps,rjunk
integer NE,i,Nelect,Nk,Nb,k,b,Ntot,Natoms,ijunk,ispin,FixedMu
character ajunk
character*35 LINE
real*16, allocatable :: E1(:)
real*16, allocatable :: E2(:)
real*16, allocatable :: F1(:)
real*16, allocatable :: F2(:)
real*16, allocatable :: S(:)
real*16, allocatable :: W(:)
real*16 GetEF, GetMu

call parse(T0,T1,dT,mu,FixedMu)
read(*,*) Natoms,ijunk,ijunk,ispin
do i=1,4
   read(*,*)
enddo
read(*,*) Nelect,Nk,Nb
Ntot=Nk*Nb
!! W is weight, E1 and E2 are spin up and spin down eigenvalues
allocate(E1(Ntot),E2(Ntot),F1(Ntot),F2(Ntot),S(Ntot),W(Ntot))
kB=1./11604.5 !! eV/K
eps=1.E-10 !! Fudge factor

NE=0
do k=1,Nk
   read(*,*)
   read(*,*) rjunk,rjunk,rjunk,wk
   do b=1,Nb
      read(*,'(A)') LINE
      NE=NE+1
      if(ispin.eq.2) then
         read(LINE,*) ajunk,Eb1,Eb2
         E1(NE)=Eb1
         E2(NE)=Eb2
      elseif(ispin.eq.1) then
         read(LINE,*) ajunk,Eb1
         E1(NE)=Eb1
         E2(NE)=Eb1
      else
         write(*,*) "Error ispin ",ispin
         write(*,*) LINE
         stop
      endif
      W(NE)=wk
   enddo
enddo

if (NE /= Nk*Nb) then
   write(*,*) "Counting error NE /= Nk*Nb ",NE,Nk,Nb
   stop
endif

if(mu.eq.123456789.) mu=getEF()  !! signal EF not set by hand
U0=(SUM(W*E1,mask=E1.le.mu)+SUM(W*E2,mask=E2.le.mu))/float(Natoms)
!! Correct possible miscount
U0=U0-mu*(SUM(W,mask=E1.le.mu)+SUM(W,mask=E2.le.mu)-Nelect)/float(Natoms)

do T=T0,T1,dT
   if((FixedMu.ne.1).and.(T.gt.100.)) mu=GetMu(Nelect,Ntot,W,kB*T,E1,E2,mu)
   F1=Fermi(kB*T,mu,E1)
   F2=Fermi(kB*T,mu,E2)
   Etot=(SUM(W*E1*F1)+SUM(W*E2*F2))/float(Natoms)-U0
!! Fudge factor eps to avoid log of negative number
   F1=eps+(1.-2.*eps)*F1
   F2=eps+(1.-2.*eps)*F2
   S=F1*log(F1)+(1.-F1)*log(1.-F1)+F2*log(F2)+(1.-F2)*log(1.-F2)
   Stot=-kB*SUM(W*S)/float(Natoms)
   write(*,'(6f12.6)') T,Etot-T*Stot,Etot,-T*Stot,mu
enddo
deallocate(W,F1,F2,E1,E2)
stop
end

real*16 function GetEF()
implicit none
real*16 rjunk
integer i,ijunk
LOGICAL :: doscar_exists
INQUIRE(FILE="DOSCAR", EXIST=doscar_exists)
if(doscar_exists) then
   open(1,file="DOSCAR")
   do i=1,5
      read(1,*)
   enddo
   read(1,*) rjunk,rjunk,ijunk,GetEF,rjunk
   close(1)
else
   write(0,*) "DOSCAR not found. Use option -EF EF to set initial mu"
   stop
end if
return
end

real*16 function GetMu(Nelect,Ntot,W,kT,E1,E2,mu)
  use MyFermi
implicit none
real*16 kT,mu,E1(Ntot),E2(Ntot),W(Ntot)
integer Nelect,Ntot,iter
real*16 I,Iprime,tol
tol=1.e-5
  do iter=1,10
     I=sum(W*Fermi(kT,mu,E1))+sum(W*Fermi(kT,mu,E2))
     if(abs(I-Nelect).lt.tol) then
        GetMu = mu
        return
     endif
     Iprime = sum(W*FermiP(kT,mu,E1))+sum(W*FermiP(kT,mu,E2))
     mu=mu+0.9*(Nelect-I)/Iprime  !! Newton's method
     GetMu=mu
  enddo
  write(*,'(A,f12.4,A,f12.6,A,i6,A,f12.4)') "GetMu failed at T=",kT*11604.5,",  mu=",mu,",  target N=",Nelect,",  obtained ",I
  write(*,'(A)') "Consider higher T0 or higher k-point density."
  stop
end function GetMu

subroutine parse(T0,T1,dT,mu,FixedMu)
implicit none
real*16 T0,T1,dT,mu
character*20 option,value
integer iarg,narg,ia,nargs,FixedMu
iarg=iargc()
nargs=iarg/2
T0 = 100.
T1 = 2000.
dT = 10.
mu = 123456789.  !! signal that EF is not set
if(iarg.ne.2*nargs) call syntax()
do ia=1,nargs
   call getarg(2*ia-1,option)
   call getarg(2*ia,value)
   if(option == "-T0") then
      read(value,*) T0
   else if(option == "-T1") then
      read(value,*) T1
   else if(option == "-dT") then
      read(value,*) dT
   else if(option == "-mu") then
      read(value,*) mu
   else if(option == "-FixedMu") then
      read(value,*) FixedMu
   else if(option == "-h") then
      call syntax()
   endif
enddo
return
end

subroutine syntax()
  write(*,*) "syntax: cat EIGENVAL | Felect -h -T0 (100) -T1 (2000) -dT (10) -mu EF -FixedMu 0"
  write(*,*) "output: T, F, U, -TS  [eV/atom]"
stop
end
