! ======= Set of subroutines to be used for smearing methods ===========
!
! These routines are derived from those in VASP with the generalization
! that the integration results for each k-point and band can be scaled
! by an individual weight function (e.g. as needed for optical spectra)
! if no extra weights are required supply an array WEIGHT containg "1."
! furthermore CELEN was renamed to EIGEN and its data type been changed
! from COMPLEX(q) to REAL(q), the structure WDES of VASP was removed and
! all necessary members of WDES have been listed as separate arguments,
! and Fermi function smearing replaced by Lorentzian line shape smearing.
! Anything related to calculation of Fermi energy, entropy etc. appearing
! in the original VASP routines was taken out completely (useless here!).
!

!***********************SUBROUTINE DENMP *******************************
!
! if ISMEAR=0
! subroutine DENMP calculates a continuous density of states in the
! interval (EMIN,EMAX) by applying a gaussian broadening to the discrete
! eigenvalue spectrum contained in EIGEN(NBANDS,NKPTS). The width of the
! gaussians is SIGMA.
! if ISMEAR<0 alternatively a Lorentzian smearing (width SIGMA) is used
! if ISMEAR>0 the generalized form of Methfessel and Paxton of order
!        N=ISMEAR will be used instead of Gaussians to get the dos ...
! initially it took 2 seconds to calculate occupancies (gK)
!
!***********************************************************************

     SUBROUTINE DENMP(NBANDD,NKDIM,NKPTS,NBANDS,WTKPT,ISPIN,EMIN, &
                      EMAX,ISMEAR,SIGMA,EIGEN,WEIGHT,NEDOS,DOS,DOSI)

      USE prec
      USE constant
      
      IMPLICIT COMPLEX(q) (C)
      IMPLICIT REAL(q) (A-B,D-H,O-Z)

      REAL(q) DOS(NEDOS,ISPIN),DOSI(NEDOS,ISPIN)
      REAL(q) EIGEN(NBANDD,NKDIM,ISPIN)
      REAL(q) WTKPT(NKPTS),WEIGHT(NBANDD,NKDIM,ISPIN)
! local variables
      LOGICAL LOWB,HIGHB
      INTEGER ISP

      PLUSMINUS = 8._q
      IF (ISMEAR==-1) PLUSMINUS = 8000._q
      SIGMA_=ABS(SIGMA)
      IF (SIGMA_==0) RETURN
      RSPIN=3-ISPIN

      DELTAE=(EMAX-EMIN)/(NEDOS-1)
!=======================================================================
! initialize arrays for dos and integr. dos
!=======================================================================
      DOS =0
      DOSI=0
!=======================================================================
! accumulate dos and integrated dos
!=======================================================================
      DO ISP=1,ISPIN
      DO K=1,NKPTS
      DO N=1,NBANDS
        EPS=EIGEN(N,K,ISP)
        WGT=RSPIN*WTKPT(K)*WEIGHT(N,K,ISP)

        NELOW=(EPS-PLUSMINUS*SIGMA_-EMIN)/DELTAE+1
        NEHIG=(EPS+PLUSMINUS*SIGMA_-EMIN)/DELTAE+1
        IF (NELOW<1)     NELOW=1
        IF (NELOW>NEDOS) NELOW=NEDOS
        IF (NEHIG<1)     NEHIG=1
        IF (NEHIG>NEDOS) NEHIG=NEDOS

        DO I=NELOW,NEHIG
          E=EMIN+DELTAE*(I-1)-EPS
          CALL DELSTP(ISMEAR,(E/SIGMA_),DFUN,SFUN)
          EPSDOS=DFUN/SIGMA_
          DOS(I,ISP) =DOS(I,ISP) +(WGT*EPSDOS)
          DOSI(I,ISP)=DOSI(I,ISP)+WGT*SFUN
        ENDDO
        DO I=NEHIG+1,NEDOS
          DOSI(I,ISP)=DOSI(I,ISP)+WGT
        ENDDO
      ENDDO
      ENDDO
      ENDDO

      RETURN
      END SUBROUTINE


!******************** DELSTP *******************************************
!
! Returns generalised delta and step functions (Methfessel & Paxton)
! or (not contained in original code) a Lorentzian and its integral
!
!  Input:
!      n > -1 : order of approximant; x : argument
!      n < 0  : return a Lorentzian line shape
!  Output:
!      D_n (x) ,  S_n (x)
!  Remarks:
!      D_n (x) = exp -x^2 * sum_i=0^n A_i H_2i(x)
!      S_n (x) = (1 - erf x)/2 + exp -x^2 * sum_i=1^n A_i H_{2i-1}(x)
!      where H is a Hermite polynomial and
!      A_i = (-1)^i / ( i! 4^i sqrt(pi) )
!
!***********************************************************************

      SUBROUTINE DELSTP(N,X,D,S)
      USE prec
      USE constant
      IMPLICIT REAL(q) (A-H,O-Z)

      IF (X<-1.E5_q) THEN
         D=0._q
         S=0._q
         RETURN
      END IF
      IF (X>1.E5_q) THEN
         D=0._q
         S=1._q
         RETURN
      END IF
!=======================================================================
!  If n = 0 : assume Gaussian type smearing
!=======================================================================
      IF (N==0) THEN
         D=EXP(-(X*X))/SQRT(PI)
         S=0.5_q+0.5_q*ERRF(X)
         RETURN
      END IF
!=======================================================================
!  If n < 0 : assume Lorentzian type smearing
!=======================================================================
      IF (N<0) THEN
         D=1._q/(X*X+1._q)/PI
         S=0.5_q+ATAN(X)/PI
         RETURN
      END IF
!=======================================================================
! Methfessel & Paxton
!=======================================================================
      EX2=EXP(-(X*X))
      S0=0.5_q*ERRF(X)
      A=1._q/SQRT(PI)
      K=0
      H1=1._q
      H2=2._q*X
      S=0._q
      D=A
      DO I=1,N
         A=A/((-4._q)*I)
         K=K+1
         H3=H1
         H1=H2
         H2=2._q*X*H2-2*K*H3
         S=S+A*H1
         K=K+1
         H3=H1
         H1=H2
         H2=2._q*X*H2-2*K*H3
         D=D+A*H1
      ENDDO
      D=D*EX2
      S=0.5_q+S0-S*EX2
      RETURN
      END