!
! Copyright (C) 2001-2016 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!----------------------------------------------------------------------------
!
! Usage :
! $ proj_fermi.x -in {input file}
! Then it generates proj.frmsf (for nspin = 1, 4) or
! proj1.frmsf and proj2.frmsf (for nspin = 2)
!
! Input file format (projwfc.x + tail):
! &PROJWFC
! prefix = "..."
! outdir = "..."
! ...
! /
! {Number of target WFCs}
! {Index of WFC1} {Index of WFC2} {Index of WFC3} ...
!
!-----------------------------------------------------------------------
MODULE nesting_proj_routines
!-----------------------------------------------------------------------
  !
  IMPLICIT NONE
  !
CONTAINS
  !
!-----------------------------------------------------------------------
SUBROUTINE read_projwfc_nesting(lbinary)
  !-----------------------------------------------------------------------
  !
  ! ... Read projwfc.x input file and use prefix, outdir, lbinary_data only.
  !
  USE io_files,  ONLY : prefix, tmp_dir
  USE io_global, ONLY : ionode, ionode_id
  USE mp_world,  ONLY : world_comm
  USE spin_orb,  ONLY : lforcet
  USE mp,        ONLY : mp_bcast
  USE kinds,     ONLY : DP
  USE klist,     ONLY : degauss, ngauss
  !
  IMPLICIT NONE
  !
  LOGICAL,INTENT(OUT) :: lbinary
  !
  CHARACTER(LEN=256), EXTERNAL :: trimcheck
  !
  CHARACTER (len=256) :: filpdos, filproj, outdir
  REAL (DP) :: Emin, Emax, DeltaE, ef_0
  INTEGER :: ios
  LOGICAL :: lwrite_overlaps, lbinary_data
  LOGICAL :: lsym, kresolveddos, tdosinboxes, plotboxes, pawproj
  INTEGER, PARAMETER :: N_MAX_BOXES = 999
  INTEGER :: n_proj_boxes, irmin(3,N_MAX_BOXES), irmax(3,N_MAX_BOXES)
  LOGICAL :: lgww  !if .true. use GW QP energies from file bands.dat
  !
  ! Exact the same namelist to that of projwfc.x
  !
  NAMELIST / projwfc / outdir, prefix, ngauss, degauss, lsym, &
             Emin, Emax, DeltaE, filpdos, filproj, lgww, &
             kresolveddos, tdosinboxes, n_proj_boxes, irmin, irmax, plotboxes, &
             lwrite_overlaps, lbinary_data, pawproj, lforcet, ef_0
  !
  !   set default values for variables in namelist
  !
  prefix = 'pwscf'
  CALL get_environment_variable('ESPRESSO_TMPDIR', outdir)
  IF (TRIM(outdir) == ' ') outdir = './'
  lbinary_data = .FALSE.
  !
  ios = 0
  !

  lforcet = .FALSE.

  IF (ionode) THEN
     CALL input_from_file ()
     READ (5, projwfc, iostat = ios)
     tmp_dir = trimcheck (outdir)
  ENDIF
  !
  CALL mp_bcast (ios, ionode_id, world_comm)
  IF (ios /= 0) CALL errore ('do_projwfc', 'reading projwfc namelist', ABS(ios))
  !
  ! ... Broadcast variables (Only used ones)
  !
  CALL mp_bcast(tmp_dir,   ionode_id, world_comm)
  CALL mp_bcast(prefix,    ionode_id, world_comm)
  CALL mp_bcast(lbinary_data,    ionode_id, world_comm)
  CALL mp_bcast(lforcet, ionode_id, world_comm)
  !
  lbinary = lbinary_data
  !
END SUBROUTINE read_projwfc_nesting
!
!-----------------------------------------------------------------------
SUBROUTINE read_atomic_proj_nesting(lbinary_data, wt, ns, nk)
  !-----------------------------------------------------------------------
  !
  ! Read atomic_proj.* generated by projwfc.x
  !
  USE io_files,           ONLY : prefix, tmp_dir, postfix
  USE basis,              ONLY : natomwfc
  USE wvfct,              ONLY : nbnd
  USE fermisurfer_common, ONLY : b_low, b_high
  USE io_global,          ONLY : stdout, ionode, ionode_id
  USE kinds,              ONLY : DP
  USE mp_world,           ONLY : world_comm
  USE mp,                 ONLY : mp_bcast
  USE read_proj,          ONLY : read_xml_proj
  !
  IMPLICIT NONE
  !
  INTEGER,INTENT(IN) :: nk, ns
  LOGICAL,INTENT(IN) :: lbinary_data
  REAL(DP),INTENT(OUT) :: wt(nbnd,nk,ns,2)
  !
  INTEGER :: ik, ik_eff, ibnd, ispin, ierr, nwfc(2), targetwfc(natomwfc,2), ii
  INTEGER :: nb_, nk_, ns_, nwfc_
  REAL(dp) :: nelec, ef
  REAL(dp), ALLOCATABLE :: xk(:,:), wk(:), et(:,:)
  CHARACTER(256) :: tmp
  COMPLEX(DP), ALLOCATABLE :: projs(:,:,:)
  !
  ! Read target wavefunctions from the tail of the input file
  !
  IF (ionode) THEN
     READ(5,*) nwfc(1)
     READ(5,*) targetwfc(1:nwfc(1),1)
     READ(5,*) nwfc(2)
     READ(5,*) targetwfc(1:nwfc(2),2)
  END IF
  CALL mp_bcast(nwfc,      ionode_id, world_comm)
  CALL mp_bcast(targetwfc, ionode_id, world_comm)
  WRITE(stdout,'(5x,a, i6)') "Number of target wavefunction 1 : ",  nwfc(1)
  WRITE(stdout,'(5x,a, 1000i6)') "Target wavefunction 1 : ", targetwfc(1:nwfc(1),1)
  WRITE(stdout,'(5x,a, i6)') "Number of target wavefunction 2 : ",  nwfc(2)
  WRITE(stdout,'(5x,a, 1000i6)') "Target wavefunction 2 : ", targetwfc(1:nwfc(2),2)
  !
  tmp = TRIM(tmp_dir) // TRIM(prefix) // postfix // 'atomic_proj.xml'
  !
  IF (lbinary_data) CALL infomsg('read_atomic_proj','lbinary no longer implemented')
  !
  IF (ionode) THEN
     !
     ! Read projections
     !
     CALL read_xml_proj ( tmp, ierr, nwfc_, nb_, nk_, &
          ns_, nelec, ef, xk, wk, et, projs )
     IF ( ierr /= 0 ) CALL errore('read_atomic_proj','error in reading file',1)
     IF ( ns_ == 2 ) nk_ = nk_ / 2
     IF ( (nwfc_ /= natomwfc) .OR. (nb_ /= nbnd) .OR. (nk_ /= nk) .OR. &
          (ns_ /= ns) ) CALL errore('read_atomic_proj','mismatch',1)
     !
     DO ik = 1, nk
        !
        DO ispin = 1, ns
           !
           ! Store Sum_{target tau l m n} |<n k|tau n m l>|^2 into wt(:,:,:)
           !
           ik_eff = ik + (ispin-1)*nk
           DO ibnd = b_low, b_high
              DO ii = 1, 2
                 wt(ibnd, ik, ispin, ii) = SUM(REAL(CONJG(projs(targetwfc(1:nwfc(ii),ii), ibnd, ik_eff)) &
                 &                                      * projs(targetwfc(1:nwfc(ii),ii), ibnd, ik_eff), DP))
              END DO
           END DO
           !
        END DO
        !
     END DO
     !
  END IF
  !
  CALL mp_bcast(wt, ionode_id, world_comm)
  !
END SUBROUTINE read_atomic_proj_nesting
!
SUBROUTINE nesting_proj_delta1(nbnd_fs,nspin_lsda,nirr_k,irr_k,eig1,eig2,wfc1,wfc2,nesting_factor)
  !---------------------------------------------------------------------
  !
  ! This routine computed the weight for the double-delta function.
  !
  USE kinds, ONLY : dp
  USE mp, ONLY : mp_sum
  USE mp_world, ONLY : world_comm
  USE klist,  ONLY: nkstot
  USE lsda_mod,   ONLY : nspin
  USE ktetra, ONLY : tetra, ntetra, nntetra, wlsm
  !
  INTEGER,INTENT(IN) :: nbnd_fs, nspin_lsda, nirr_k
  INTEGER,INTENT(IN) :: irr_k(3,nirr_k)
  REAL(dp),INTENT(IN) :: eig1(nbnd_fs, nkstot, nspin_lsda), eig2(nbnd_fs*nirr_k,nkstot, nspin_lsda)
  REAL(dp),INTENT(IN) :: wfc1(nbnd_fs, nkstot, nspin_lsda), wfc2(nbnd_fs*nirr_k,nkstot, nspin_lsda)
  REAL(dp),INTENT(OUT) :: nesting_factor(nirr_k,nspin_lsda)
  !
  INTEGER :: ns, nt, ii, ibnd, itetra(4), start_t, last_t
  REAL(dp) :: e(4), a(4,4), V, tsmall(3,4)
  REAL(dp),ALLOCATABLE :: ei0(:,:), ej0(:,:), ej1(:,:), pi0(:,:), pj0(:,:), pi1(:), pj1(:,:), nesting_factor1(:)
  !
  nesting_factor(1:nirr_k,1:nspin_lsda) = 0.0_dp
  !
  CALL divide(world_comm, 6 * nkstot, start_t, last_t)
  !
  !$OMP PARALLEL DEFAULT(NONE) &
  !$OMP & PRIVATE(ns,nt,ibnd,ii,e,ei0,ej0,pi0,pj0,a,V,tsmall,ej1,pi1,pj1,itetra,nesting_factor1) &
  !$OMP & SHARED(nspin_lsda,nbnd_fs,nirr_k,start_t,last_t,eig1,eig2,wfc1,wfc2,tetra,wlsm,nntetra,irr_k) &
  !$OMP & REDUCTION(+:nesting_factor)
  !
  ALLOCATE(ei0(4,nbnd_fs), ej0(4,nbnd_fs*nirr_k), ej1(3,nbnd_fs*nirr_k), pi1(3), &
  &        pi0(4,nbnd_fs), pj0(4,nbnd_fs*nirr_k), pj1(3,nbnd_fs*nirr_k), nesting_factor1(nirr_k))
  !
  DO ns = 1, nspin_lsda
     !
     !$OMP DO
     DO nt = start_t, last_t
        !
        ei0(1:4, 1:nbnd_fs) = 0.0_dp
        ej0(1:4, 1:nbnd_fs*nirr_k) = 0.0_dp
        pi0(1:4, 1:nbnd_fs) = 0.0_dp
        pj0(1:4, 1:nbnd_fs*nirr_k) = 0.0_dp
        DO ii = 1, nntetra
           !
           DO ibnd = 1, nbnd_fs
              ei0(1:4, ibnd) = ei0(1:4, ibnd) &
              &             + wlsm(1:4,ii) * eig1(ibnd, tetra(ii, nt), ns)
              pi0(1:4, ibnd) = pi0(1:4, ibnd) &
              &             + wlsm(1:4,ii) * wfc1(ibnd, tetra(ii, nt), ns)
           END DO
           !
           DO ibnd = 1, nbnd_fs*nirr_k
              ej0(1:4, ibnd) = ej0(1:4, ibnd) &
              &             + wlsm(1:4,ii) * eig2(ibnd, tetra(ii, nt), ns)
              pj0(1:4, ibnd) = pj0(1:4, ibnd) &
              &             + wlsm(1:4,ii) * wfc2(ibnd, tetra(ii, nt), ns)
           END DO
           !
        END DO
        !
        DO ibnd = 1, nbnd_fs
           !
           itetra(1) = 0
           e(1:4) = ei0(1:4, ibnd)
           call hpsort (4, e, itetra)
           !
           DO ii = 1, 4
              a(1:4,ii) = (0.0_dp - e(ii)) / (e(1:4) - e(ii))
           END DO
           !
           IF(e(1) < 0.0_dp .AND. 0.0_dp <= e(2)) THEN
              !
              ! A
              !
              !V = 3.0_dp * a(2,1) * a(3,1) * a(4,1) / (0.0_dp - e(1))
              V = 3.0_dp * a(2,1) * a(3,1)           / (e(4) - e(1))
              !
              tsmall(1, 1:4) = (/a(1,2), a(2,1), 0.0_dp, 0.0_dp/)
              tsmall(2, 1:4) = (/a(1,3), 0.0_dp, a(3,1), 0.0_dp/)
              tsmall(3, 1:4) = (/a(1,4), 0.0_dp, 0.0_dp, a(4,1)/)
              !
              pi1(1:3                 ) = MATMUL(tsmall(1:3,1:4), pi0(itetra(1:4), ibnd))
              ej1(1:3,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:3,1:4), ej0(itetra(1:4), 1:nbnd_fs*nirr_k))
              pj1(1:3,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:3,1:4), pj0(itetra(1:4), 1:nbnd_fs*nirr_k))
              !
              CALL nesting_proj_delta2(nbnd_fs,nirr_k,irr_k,ej1,pi1,pj1,nesting_factor1)
              !
              nesting_factor(1:nirr_k,ns) = nesting_factor(1:nirr_k,ns) + V * nesting_factor1(1:nirr_k)
              !
           ELSE IF( e(2) < 0.0_dp .AND. 0.0_dp <= e(3)) THEN
              !
              ! B - 1
              !
              !V = 3.0_dp * a(3,1) * a(4,1) * a(2,4) / (0.0_dp - e(1))
              V = 3.0_dp           * a(4,1) * a(2,4) / (e(3) - e(1))
              !
              tsmall(1, 1:4) = (/a(1,3), 0.0_dp, a(3,1), 0.0_dp/)
              tsmall(2, 1:4) = (/a(1,4), 0.0_dp, 0.0_dp, a(4,1)/)
              tsmall(3, 1:4) = (/0.0_dp, a(2,4), 0.0_dp, a(4,2)/)
              !
              pi1(1:3                 ) = MATMUL(tsmall(1:3,1:4), pi0(itetra(1:4), ibnd))
              ej1(1:3,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:3,1:4), ej0(itetra(1:4), 1:nbnd_fs*nirr_k))
              pj1(1:3,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:3,1:4), pj0(itetra(1:4), 1:nbnd_fs*nirr_k))
              !
              CALL nesting_proj_delta2(nbnd_fs,nirr_k,irr_k,ej1,pi1,pj1,nesting_factor1)
              !
              nesting_factor(1:nirr_k,ns) = nesting_factor(1:nirr_k,ns) + V * nesting_factor1(1:nirr_k)
              !
              ! B - 2
              !
              !V = 3.0_dp * a(2,3) * a(3,1) * a(4,2) / (0.0_dp - e(1))
              V = 3.0_dp * a(2,3)           * a(4,2) / (e(3) - e(1))
              !
              tsmall(1, 1:4) = (/a(1,3), 0.0_dp, a(3,1), 0.0_dp/)
              tsmall(2, 1:4) = (/0.0_dp, a(2,3), a(3,2), 0.0_dp/)
              tsmall(3, 1:4) = (/0.0_dp, a(2,4), 0.0_dp, a(4,2)/)
              !
              pi1(1:3                 ) = MATMUL(tsmall(1:3,1:4), pi0(itetra(1:4), ibnd))
              ej1(1:3,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:3,1:4), ej0(itetra(1:4), 1:nbnd_fs*nirr_k))
              pj1(1:3,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:3,1:4), pj0(itetra(1:4), 1:nbnd_fs*nirr_k))
              !
              CALL nesting_proj_delta2(nbnd_fs,nirr_k,irr_k,ej1,pi1,pj1,nesting_factor1)
              !
              nesting_factor(1:nirr_k,ns) = nesting_factor(1:nirr_k,ns) + V * nesting_factor1(1:nirr_k)
              !
           ELSE IF(e(3) < 0.0_dp .AND. 0.0_dp < e(4)) THEN
              !
              ! C
              !
              !V = 3.0_dp * a(1,4) * a(2,4) * a(3,4) / (e(4) - 0.0_dp)
              V = 3.0_dp * a(1,4) * a(2,4)           / (e(4) - e(3))
              !
              tsmall(1, 1:4) = (/a(1,4), 0.0_dp, 0.0_dp, a(4,1)/)
              tsmall(2, 1:4) = (/0.0_dp, a(2,4), 0.0_dp, a(4,2)/)
              tsmall(3, 1:4) = (/0.0_dp, 0.0_dp, a(3,4), a(4,3)/)
              !
              pi1(1:3                 ) = MATMUL(tsmall(1:3,1:4), pi0(itetra(1:4), ibnd       ))
              ej1(1:3,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:3,1:4), ej0(itetra(1:4), 1:nbnd_fs*nirr_k))
              pj1(1:3,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:3,1:4), pj0(itetra(1:4), 1:nbnd_fs*nirr_k))
              !
              CALL nesting_proj_delta2(nbnd_fs,nirr_k,irr_k,ej1,pi1,pj1,nesting_factor1)
              !
              nesting_factor(1:nirr_k,ns) = nesting_factor(1:nirr_k,ns) + V * nesting_factor1(1:nirr_k)
              !
           END IF
           !
        END DO
        !
     END DO ! nt
     !$OMP END DO
     !
  END DO ! ns
  !
  DEALLOCATE(ei0, ej0, ej1, pi0, pj0, pi1, pj1, nesting_factor1)
  !
  !$OMP END PARALLEL
  !
  nesting_factor(1:nirr_k,1:nspin_lsda) = nesting_factor(1:nirr_k,1:nspin_lsda) / REAL(ntetra, dp)
  IF(nspin == 1) nesting_factor(1:nirr_k,1) = 2.0_dp * nesting_factor(1:nirr_k,1)
  !
  CALL mp_sum(nesting_factor, world_comm)
  !
END SUBROUTINE nesting_proj_delta1
!
!-------------------------------------------------------------------------------
SUBROUTINE nesting_proj_delta2(nbnd_fs,nirr_k,irr_k,ej0,pi0,pj0,nesting_factor)
  !-----------------------------------------------------------------------------
  !
  ! 2nd step of tetrahedra method.
  !
  USE kinds, ONLY : dp
  !
  INTEGER,INTENT(IN) :: nbnd_fs, nirr_k
  INTEGER,INTENT(IN) :: irr_k(3,nirr_k)
  REAL(dp),INTENT(IN) :: ej0(3,nbnd_fs,nirr_k), pi0(3), pj0(3,nbnd_fs,nirr_k)
  REAL(dp),INTENT(OUT) :: nesting_factor(nirr_k)
  !
  INTEGER :: ibnd, itetra(3), ii, ik
  REAL(dp) :: e(3), a(3,3), V, w(3)
  !
  nesting_factor(1:nirr_k) = 0.0_dp
  !
  DO ik = 1, nirr_k
     !
     IF(ALL(irr_k(1:3,ik) == 0)) CYCLE
     !
     DO ibnd = 1, nbnd_fs
        !
        IF(MAXVAL(ABS(ej0(1:3,ibnd,ik))) < 1e-10_dp) &
        & CALL errore("nesting_delta2", "Nesting occurs.", ibnd)
        !
        itetra(1) = 0
        e(1:3) = ej0(1:3,ibnd,ik)
        call hpsort (3, e, itetra)
        !
        DO ii = 1, 3
           a(1:3,ii) = (0.0_dp - e(ii)) / (e(1:3) - e(ii))
        END DO
        !
        IF((e(1) < 0.0_dp .AND. 0.0_dp <= e(2)) .OR. (e(1) <= 0.0_dp .AND. 0.0_dp < e(2))) THEN
           !
           !V = a(2,1) * a(3,1) / (0.0_dp - e(1))
           V = a(2,1)           / (e(3) - e(1))
           !
           w(itetra(1)) = V * (a(1,2) + a(1,3))
           w(itetra(2)) = V * a(2,1)
           w(itetra(3)) = V * a(3,1)
           !
           nesting_factor(ik) = nesting_factor(ik) + SUM(w(1:3)*pi0(1:3)*pj0(1:3,ibnd,ik))
           !
        ELSE IF((e(2) <= 0.0_dp .AND. 0.0_dp < e(3)) .OR. (e(2) < 0.0_dp .AND. 0.0_dp <= e(3))) THEN
           !
           !V = a(1,3) * a(2,3) / (e(3) - 0.0_dp)
           V = a(1,3)           / (e(3) - e(2))
           !
           w(itetra(1)) = V * a(1,3)
           w(itetra(2)) = V * a(2,3)
           w(itetra(3)) = V * (a(3,1) + a(3,2))
           !
           nesting_factor(ik) = nesting_factor(ik) + SUM(w(1:3)*pi0(1:3)*pj0(1:3,ibnd,ik))
           !
        END IF
        !
     END DO ! ib
     !
  END DO
  !
END SUBROUTINE nesting_proj_delta2
!
SUBROUTINE nesting_proj_theta1(nbnd_fs,nspin_lsda,nirr_k,eig1,eig2,wfc1,wfc2,nesting_factor)
  !---------------------------------------------------------------------
  !
  ! This routine computed the weight for the double-delta function.
  !
  USE kinds, ONLY : dp
  USE mp, ONLY : mp_sum
  USE mp_world, ONLY : world_comm
  USE klist,  ONLY: nkstot
  USE lsda_mod,   ONLY : nspin
  USE ktetra, ONLY : tetra, ntetra, nntetra, wlsm
  !
  INTEGER,INTENT(IN) :: nbnd_fs, nspin_lsda, nirr_k
  REAL(dp),INTENT(IN) :: eig1(nbnd_fs, nkstot, nspin_lsda), eig2(nbnd_fs*nirr_k,nkstot, nspin_lsda), &
  &                      wfc1(nbnd_fs, nkstot, nspin_lsda), wfc2(nbnd_fs*nirr_k,nkstot, nspin_lsda)
  REAL(dp),INTENT(OUT) :: nesting_factor(nirr_k,nspin_lsda)
  !
  INTEGER :: ns, nt, ii, ibnd, itetra(4), start_t, last_t
  REAL(dp) :: e(4), a(4,4), thr = 1.0e-8_dp, V, tsmall(4,4), ei1(4), pi1(4)
  REAL(dp),ALLOCATABLE :: ei0(:,:), ej0(:,:), ej1(:,:), pi0(:,:), pj0(:,:), pj1(:,:), nesting_factor1(:)
  !
  nesting_factor(1:nirr_k,1:nspin_lsda) = 0.0_dp
  !
  CALL divide(world_comm, 6 * nkstot, start_t, last_t)
  !
  !$OMP PARALLEL DEFAULT(NONE) &
  !$OMP & PRIVATE(ns,nt,ibnd,ii,e,ei0,ej0,pi0,pj0,a,V,tsmall,ei1,ej1,pi1,pj1,itetra,nesting_factor1) &
  !$OMP & SHARED(nspin_lsda,nbnd_fs,nirr_k,start_t,last_t,eig1,eig2,wfc1,wfc2,tetra,wlsm,nntetra,thr) &
  !$OMP & REDUCTION(+:nesting_factor)
  !
  ALLOCATE(ei0(4,nbnd_fs), ej0(4,nbnd_fs*nirr_k), ej1(4,nbnd_fs*nirr_k), nesting_factor1(nirr_k), &
  &        pi0(4,nbnd_fs), pj0(4,nbnd_fs*nirr_k), pj1(4,nbnd_fs*nirr_k))
  !
  DO ns = 1, nspin_lsda
     !
     !$OMP DO
     DO nt = start_t, last_t
        !
        ei0(1:4, 1:nbnd_fs) = 0.0_dp
        ej0(1:4, 1:nbnd_fs*nirr_k) = 0.0_dp
        pi0(1:4, 1:nbnd_fs) = 0.0_dp
        pj0(1:4, 1:nbnd_fs*nirr_k) = 0.0_dp
        DO ii = 1, nntetra
           !
           DO ibnd = 1, nbnd_fs
              ei0(1:4, ibnd) = ei0(1:4, ibnd) &
              &             + wlsm(1:4,ii) * eig1(ibnd, tetra(ii, nt), ns)
              pi0(1:4, ibnd) = pi0(1:4, ibnd) &
              &             + wlsm(1:4,ii) * wfc1(ibnd, tetra(ii, nt), ns)
           END DO
           !
           DO ibnd = 1, nbnd_fs*nirr_k
              ej0(1:4, ibnd) = ej0(1:4, ibnd) &
              &             + wlsm(1:4,ii) * eig2(ibnd, tetra(ii, nt), ns)
              pj0(1:4, ibnd) = pj0(1:4, ibnd) &
              &             + wlsm(1:4,ii) * wfc2(ibnd, tetra(ii, nt), ns)
           END DO
           !
        END DO
        !
        DO ibnd = 1, nbnd_fs
           !
           itetra(1) = 0
           e(1:4) = ei0(1:4, ibnd)
           call hpsort (4, e, itetra)
           !
           DO ii = 1, 4
              a(1:4,ii) = (0.0_dp - e(ii)) / (e(1:4) - e(ii))
           END DO
           !
           IF(e(1) <= 0.0_dp .AND. 0.0_dp < e(2)) THEN
              !
              ! A - 1
              !
              V = a(2,1) * a(3,1) * a(4,1)
              !
              IF(V > thr) THEN
                 !
                 tsmall(1, 1:4) = (/1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/)
                 tsmall(2, 1:4) = (/a(1,2), a(2,1), 0.0_dp, 0.0_dp/)
                 tsmall(3, 1:4) = (/a(1,3), 0.0_dp, a(3,1), 0.0_dp/)
                 tsmall(4, 1:4) = (/a(1,4), 0.0_dp, 0.0_dp, a(4,1)/)
                 !
                 ei1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
                 pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4), ibnd))
                 ej1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), ej0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 pj1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 !
                 CALL nesting_proj_theta2(nbnd_fs,nirr_k,ei1,ej1,pi1,pj1,nesting_factor1)
                 !
                 nesting_factor(1:nirr_k,ns) = nesting_factor(1:nirr_k,ns) + V * nesting_factor1(1:nirr_k)
                 !
              END IF
              !
           ELSE IF( e(2) <= 0.0_dp .AND. 0.0_dp < e(3)) THEN
              !
              ! B - 1
              !
              V = a(3,1) * a(4,1) * a(2,4)
              !
              IF(V > thr) THEN
                 !
                 tsmall(1, 1:4) = (/1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/)
                 tsmall(2, 1:4) = (/a(1,3), 0.0_dp, a(3,1), 0.0_dp/)
                 tsmall(3, 1:4) = (/a(1,4), 0.0_dp, 0.0_dp, a(4,1)/)
                 tsmall(4, 1:4) = (/0.0_dp, a(2,4), 0.0_dp, a(4,2)/)
                 !
                 ei1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
                 pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4), ibnd))
                 ej1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), ej0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 pj1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 !
                 CALL nesting_proj_theta2(nbnd_fs,nirr_k,ei1,ej1,pi1,pj1,nesting_factor1)
                 !
                 nesting_factor(1:nirr_k,ns) = nesting_factor(1:nirr_k,ns) + V * nesting_factor1(1:nirr_k)
                 !
              END IF
              !
              ! B - 2
              !
              V = a(3,2) * a(4,2)
              !
              IF(V > thr) THEN
                 !
                 tsmall(1, 1:4) = (/1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/)
                 tsmall(2, 1:4) = (/0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp/)
                 tsmall(3, 1:4) = (/0.0_dp, a(2,3), a(3,2), 0.0_dp/)
                 tsmall(4, 1:4) = (/0.0_dp, a(2,4), 0.0_dp, a(4,2)/)
                 !
                 ei1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
                 pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4), ibnd))
                 ej1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), ej0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 pj1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 !
                 CALL nesting_proj_theta2(nbnd_fs,nirr_k,ei1,ej1,pi1,pj1,nesting_factor1)
                 !
                 nesting_factor(1:nirr_k,ns) = nesting_factor(1:nirr_k,ns) + V * nesting_factor1(1:nirr_k)
                 !
              END IF
              !
              ! B - 3
              !
              V = a(2,3) * a(3,1) * a(4,2)
              !
              IF(V > thr) THEN
                 !
                 tsmall(1, 1:4) = (/1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/)
                 tsmall(2, 1:4) = (/a(1,3), 0.0_dp, a(3,1), 0.0_dp/)
                 tsmall(3, 1:4) = (/0.0_dp, a(2,3), a(3,2), 0.0_dp/)
                 tsmall(4, 1:4) = (/0.0_dp, a(2,4), 0.0_dp, a(4,2)/)
                 !
                 ei1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
                 pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4), ibnd))
                 ej1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), ej0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 pj1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 !
                 CALL nesting_proj_theta2(nbnd_fs,nirr_k,ei1,ej1,pi1,pj1,nesting_factor1)
                 !
                 nesting_factor(1:nirr_k,ns) = nesting_factor(1:nirr_k,ns) + V * nesting_factor1(1:nirr_k)
                 !
              END IF
              !
           ELSE IF( e(3) <= 0.0_dp .AND. 0.0_dp < e(4)) THEN
              !
              ! C - 1
              !
              V = a(4,3)
              !
              IF(V > thr) THEN
                 !
                 tsmall(1, 1:4) = (/1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/)
                 tsmall(2, 1:4) = (/0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp/)
                 tsmall(3, 1:4) = (/0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp/)
                 tsmall(4, 1:4) = (/0.0_dp, 0.0_dp, a(3,4), a(4,3)/)
                 !
                 ei1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
                 pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4), ibnd))
                 ej1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), ej0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 pj1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 !
                 CALL nesting_proj_theta2(nbnd_fs,nirr_k,ei1,ej1,pi1,pj1,nesting_factor1)
                 !
                 nesting_factor(1:nirr_k,ns) = nesting_factor(1:nirr_k,ns) + V * nesting_factor1(1:nirr_k)
                 !
              END IF
              !
              ! C - 2
              !
              V = a(3,4) * a(4,2)
              !
              IF(V > thr) THEN
                 !
                 tsmall(1, 1:4) = (/1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/)
                 tsmall(2, 1:4) = (/0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp/)
                 tsmall(3, 1:4) = (/0.0_dp, a(2,4), 0.0_dp, a(4,2)/)
                 tsmall(4, 1:4) = (/0.0_dp, 0.0_dp, a(3,4), a(4,3)/)
                 !
                 ei1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
                 pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4), ibnd))
                 ej1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), ej0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 pj1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 !
                 CALL nesting_proj_theta2(nbnd_fs,nirr_k,ei1,ej1,pi1,pj1,nesting_factor1)
                 !
                 nesting_factor(1:nirr_k,ns) = nesting_factor(1:nirr_k,ns) + V * nesting_factor1(1:nirr_k)
                 !
              END IF
              !
              ! C - 3
              !
              V = a(3,4) * a(2,4) * a(4,1)
              !
              IF(V > thr) THEN
                 !
                 tsmall(1, 1:4) = (/1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/)
                 tsmall(2, 1:4) = (/a(1,4), 0.0_dp, 0.0_dp, a(4,1)/)
                 tsmall(3, 1:4) = (/0.0_dp, a(2,4), 0.0_dp, a(4,2)/)
                 tsmall(4, 1:4) = (/0.0_dp, 0.0_dp, a(3,4), a(4,3)/)
                 !
                 ei1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
                 pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4), ibnd))
                 ej1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), ej0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 pj1(1:4,1:nbnd_fs*nirr_k) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4), 1:nbnd_fs*nirr_k))
                 !
                 CALL nesting_proj_theta2(nbnd_fs,nirr_k,ei1,ej1,pi1,pj1,nesting_factor1)
                 !
                 nesting_factor(1:nirr_k,ns) = nesting_factor(1:nirr_k,ns) + V * nesting_factor1(1:nirr_k)
                 !
              END IF
              !
           ELSE IF( e(4) <= 0.0_dp ) THEN
              !
              ! D - 1
              !
              ei1(1:4) = e(1:4)
              pi1(1:4) = pi0(itetra(1:4),ibnd)
              ej1(1:4,1:nbnd_fs*nirr_k) = ej0(itetra(1:4), 1:nbnd_fs*nirr_k)
              pj1(1:4,1:nbnd_fs*nirr_k) = pj0(itetra(1:4), 1:nbnd_fs*nirr_k)
              !
              CALL nesting_proj_theta2(nbnd_fs,nirr_k,ei1,ej1,pi1,pj1,nesting_factor1)
              !
              nesting_factor(1:nirr_k,ns) = nesting_factor(1:nirr_k,ns) + nesting_factor1(1:nirr_k)
              !
           END IF
           !
        END DO
        !
     END DO ! nt
     !$OMP END DO
     !
  END DO ! ns
  !
  DEALLOCATE(ei0, ej0, ej1, nesting_factor1, pi0, pj0, pj1)
  !
  !$OMP END PARALLEL
  !
  nesting_factor(1:nirr_k,1:nspin_lsda) = nesting_factor(1:nirr_k,1:nspin_lsda) / REAL(ntetra, dp)
  IF(nspin == 1) nesting_factor(1:nirr_k,1) = 2.0_dp * nesting_factor(1:nirr_k,1)
  !
  CALL mp_sum(nesting_factor, world_comm)
  !
END SUBROUTINE nesting_proj_theta1
!
!-------------------------------------------------------------------------------
SUBROUTINE nesting_proj_theta2(nbnd_fs,nirr_k,ei0,ej0,pi0,pj0,nesting_factor)
  !-----------------------------------------------------------------------------
  !
  ! 2nd step of tetrahedra method.
  !
  USE kinds, ONLY : dp
  !
  INTEGER,INTENT(IN) :: nbnd_fs, nirr_k
  REAL(dp),INTENT(IN) :: ei0(4), ej0(4,nbnd_fs,nirr_k), pi0(4), pj0(4,nbnd_fs,nirr_k)
  REAL(dp),INTENT(OUT) :: nesting_factor(nirr_k)
  !
  INTEGER :: ibnd, itetra(4), ii, ik
  REAL(dp) :: e(4), a(4,4), V, ei1(4), ej1(4), pi1(4), pj1(4), thr = 1.0e-8_dp, &
  &           tsmall(4,4), nesting_factor1
  !
  nesting_factor(1:nirr_k) = 0.0_dp
  !
  DO ik = 1, nirr_k
     DO ibnd = 1, nbnd_fs
        !
        itetra(1) = 0
        e(1:4) = ej0(1:4,ibnd,ik)
        call hpsort (4, e, itetra)
        !
        DO ii = 1, 4
           a(1:4,ii) = (0.0_dp - e(ii)) / (e(1:4) - e(ii))
        END DO
        !
        IF(0_dp <= e(1)) THEN
             !
             ! A - 1
             !
             ej1(1:4) = e(1:4)
             ei1(1:4) = ei0(itetra(1:4))
             pj1(1:4) = pj0(itetra(1:4),ibnd,ik)
             pi1(1:4) = pi0(itetra(1:4))
             !
             CALL nesting_proj_lindhard(ei1,ej1,pi1,pj1,nesting_factor1)
             !
             nesting_factor(ik) = nesting_factor(ik) + nesting_factor1
             !
        ELSE IF((e(1) < 0.0_dp .AND. 0.0_dp <= e(2)) .OR. (e(1) <= 0.0_dp .AND. 0.0_dp < e(2))) THEN
           !
           ! B - 1
           !
           V = a(1,2)
           !
           IF(V > thr) THEN
              !
              tsmall(1, 1:4) = (/a(1,2), a(2,1), 0.0_dp, 0.0_dp/)
              tsmall(2, 1:4) = (/0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp/)
              tsmall(3, 1:4) = (/0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp/)
              tsmall(4, 1:4) = (/0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp/)
              !
              ej1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
              ei1(1:4) = MATMUL(tsmall(1:4,1:4), ei0(itetra(1:4)))
              pj1(1:4) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4),ibnd,ik))
              pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4)))
              !
              CALL nesting_proj_lindhard(ei1,ej1,pi1,pj1,nesting_factor1)
              !
              nesting_factor(ik) = nesting_factor(ik) + V * nesting_factor1
              !       
           END IF
           !
           ! B - 2
           !
           V = a(1,3) * a(2,1)
           !
           IF(V > thr) THEN
              !
              tsmall(1, 1:4) = (/a(1,2), a(2,1), 0.0_dp, 0.0_dp/)
              tsmall(2, 1:4) = (/a(1,3), 0.0_dp, a(3,1), 0.0_dp/)
              tsmall(3, 1:4) = (/0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp/)
              tsmall(4, 1:4) = (/0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp/)
              !
              ej1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
              ei1(1:4) = MATMUL(tsmall(1:4,1:4), ei0(itetra(1:4)))
              pj1(1:4) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4),ibnd,ik))
              pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4)))
              !
              CALL nesting_proj_lindhard(ei1,ej1,pi1,pj1,nesting_factor1)
              !
              nesting_factor(ik) = nesting_factor(ik) + nesting_factor1
              !
           END IF
           !
           ! B - 3
           !
           V = a(1,4) * a(2,1) * a(3,1)
           !
           IF(V > thr) THEN
              !
              tsmall(1, 1:4) = (/a(1,2), a(2,1), 0.0_dp, 0.0_dp/)
              tsmall(2, 1:4) = (/a(1,3), 0.0_dp, a(3,1), 0.0_dp/)
              tsmall(3, 1:4) = (/a(1,4), 0.0_dp, 0.0_dp, a(4,1)/)
              tsmall(4, 1:4) = (/0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp/)
              !
              ej1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
              ei1(1:4) = MATMUL(tsmall(1:4,1:4), ei0(itetra(1:4)))
              pj1(1:4) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4),ibnd,ik))
              pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4)))
              !
              CALL nesting_proj_lindhard(ei1,ej1,pi1,pj1,nesting_factor1)
              !
              nesting_factor(ik) = nesting_factor(ik) + nesting_factor1
              !       
           END IF
           !          
        ELSE IF((e(2) < 0.0_dp .AND. 0.0_dp <= e(3)) .OR. (e(2) <= 0.0_dp .AND. 0.0_dp < e(3))) THEN
           !          
           ! C - 1
           !
           V = a(2,4) * a(1,4) * a(3,1)
           !
           IF(V > thr) THEN
              !
              tsmall(1, 1:4) = (/a(1,3), 0.0_dp, a(3,1), 0.0_dp/)
              tsmall(2, 1:4) = (/a(1,4), 0.0_dp, 0.0_dp, a(4,1)/)
              tsmall(3, 1:4) = (/0.0_dp, a(2,4), 0.0_dp, a(4,2)/)
              tsmall(4, 1:4) = (/0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp/)
              !
              ej1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
              ei1(1:4) = MATMUL(tsmall(1:4,1:4), ei0(itetra(1:4)))
              pj1(1:4) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4),ibnd,ik))
              pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4)))
              !
              CALL nesting_proj_lindhard(ei1,ej1,pi1,pj1,nesting_factor1)
              !
              nesting_factor(ik) = nesting_factor(ik) + nesting_factor1
              !      
           END IF
           !
           ! C - 2
           !
           V = a(1,3) * a(2,3)
           !
           IF(V > thr) THEN
              !
              tsmall(1, 1:4) = (/a(1,3), 0.0_dp, a(3,1), 0.0_dp/)
              tsmall(2, 1:4) = (/0.0_dp, a(2,3), a(3,2), 0.0_dp/)
              tsmall(3, 1:4) = (/0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp/)
              tsmall(4, 1:4) = (/0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp/)
              !
              ej1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
              ei1(1:4) = MATMUL(tsmall(1:4,1:4), ei0(itetra(1:4)))
              pj1(1:4) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4),ibnd,ik))
              pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4)))
              !
              CALL nesting_proj_lindhard(ei1,ej1,pi1,pj1,nesting_factor1)
              !
              nesting_factor(ik) = nesting_factor(ik) + nesting_factor1
              !
           END IF
           !
           ! C - 3
           ! 
           V = a(1,3) * a(2,4) * a(3,2)
           !
           IF(V > thr) THEN
              !
              tsmall(1, 1:4) = (/a(1,3), 0.0_dp, a(3,1), 0.0_dp/)
              tsmall(2, 1:4) = (/0.0_dp, a(2,3), a(3,2), 0.0_dp/)
              tsmall(3, 1:4) = (/0.0_dp, a(2,4), 0.0_dp, a(4,2)/)
              tsmall(4, 1:4) = (/0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp/)
              !
              ej1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
              ei1(1:4) = MATMUL(tsmall(1:4,1:4), ei0(itetra(1:4)))
              pj1(1:4) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4),ibnd,ik))
              pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4)))
              !
              CALL nesting_proj_lindhard(ei1,ej1,pi1,pj1,nesting_factor1)
              !
              nesting_factor(ik) = nesting_factor(ik) + nesting_factor1
              !
           END IF
           !          
        ELSE IF((e(3) < 0.0_dp .AND. 0.0_dp <= e(4)) .OR. (e(3) <= 0.0_dp .AND. 0.0_dp < e(4))) THEN
           !
           ! D - 1
           !
           V = a(3,4) * a(2,4) * a(1,4) 
           !          
           IF(V > thr) THEN
              !
              tsmall(1, 1:4) = (/a(1,4), 0.0_dp, 0.0_dp, a(4,1)/)
              tsmall(2, 1:4) = (/0.0_dp, a(2,4), 0.0_dp, a(4,2)/)
              tsmall(3, 1:4) = (/0.0_dp, 0.0_dp, a(3,4), a(4,3)/)
              tsmall(4, 1:4) = (/0.0_dp, 0.0_dp, 0.0_dp, 1.0_dp/)
              !
              ej1(1:4) = MATMUL(tsmall(1:4,1:4), e(1:4))
              ei1(1:4) = MATMUL(tsmall(1:4,1:4), ei0(itetra(1:4)))
              pj1(1:4) = MATMUL(tsmall(1:4,1:4), pj0(itetra(1:4),ibnd,ik))
              pi1(1:4) = MATMUL(tsmall(1:4,1:4), pi0(itetra(1:4)))
              !
              CALL nesting_proj_lindhard(ei1,ej1,pi1,pj1,nesting_factor1)
              !
              nesting_factor(ik) = nesting_factor(ik) + nesting_factor1
              !        
           END IF
           !
        END IF
        !
     END DO ! ib
     !
  END DO
  !
END SUBROUTINE nesting_proj_theta2
!
SUBROUTINE nesting_proj_lindhard(ei,ej,pi,pj,nesting_factor)
  !---------------------------------------------------------------
  !
  ! This routine compute 1 / (e_{k+q} - e_{k})
  !
  USE kinds, ONLY : dp
  USE dfpt_tetra_mod, ONLY : dfpt_tetra_lindhard_1234, dfpt_tetra_lindhard_1231, &
  &                          dfpt_tetra_lindhard_1233, dfpt_tetra_lindhard_1221, &
  &                          dfpt_tetra_lindhard_1222, dfpt_tetra_lindhard_1211
  !
  IMPLICIT NONE
  !
  REAL(dp),INTENT(IN) :: ei(4), ej(4), pi(4), pj(4)
  REAL(dp),INTENT(OUT) :: nesting_factor
  !
  INTEGER :: ii, itetra(4)
  REAL(dp) :: e(4), le(4), thr, thr2, w(4)
  !
  w = 0.0_dp
  !
  itetra(1) = 0
  e(1:4) = ej(1:4) - ei(1:4)
  call hpsort (4, e, itetra)
  !
  thr = MAXVAL(e(1:4)) * 1e-3_dp
  thr2 = 1e-9_dp
  !
  DO ii = 1, 4
     IF(e(ii) < thr2) THEN
        IF(ii == 3) THEN
           WRITE(*,*) e(1:4)
           CALL errore("nesting_lindhard", "Nesting occurs.", 1)
        END IF
        le(ii) = 0.0_dp
        e(ii) = 0.0_dp
     ELSE
        le(ii) = LOG(e(ii))
     END IF
  END DO
  !
  IF(ABS(e(4) - e(3)) < thr ) THEN
     IF(ABS(e(4) - e(2)) < thr ) THEN
        IF(ABS(e(4) - e(1)) < thr ) THEN
           !
           ! e(4) = e(3) = e(2) = e(1)
           !
           w(itetra(4)) = 0.25_dp / e(4)
           w(itetra(3)) = w(itetra(4))
           w(itetra(2)) = w(itetra(4))
           w(itetra(1)) = w(itetra(4))
           !
        ELSE
           !
           ! e(4) = e(3) = e(2)
           !
           w(itetra(4)) = dfpt_tetra_lindhard_1211(e(4),e(1),le(4),le(1))
           w(itetra(3)) = w(itetra(4))
           w(itetra(2)) = w(itetra(4))
           w(itetra(1)) = dfpt_tetra_lindhard_1222(e(1),e(4),le(1),le(4))
           !
           IF(ANY(w(itetra(1:4)) < 0.0_dp)) THEN
              WRITE(*,'(100e15.5)') e(1:4)
              WRITE(*,'(100e15.5)') w(itetra(1:4))
              CALL errore("dfpt_tetra_lindhard", "Case : 4=3=2", 1)
           END IF
           !
        END IF
     ELSE IF(ABS(e(2) - e(1)) < thr ) THEN
        !
        ! e(4) = e(3), e(2) = e(1)
        !
        w(itetra(4)) = dfpt_tetra_lindhard_1221(e(4),e(2), le(4),le(2))
        w(itetra(3)) = w(itetra(4))
        w(itetra(2)) = dfpt_tetra_lindhard_1221(e(2),e(4), le(2),le(4))
        w(itetra(1)) = w(itetra(2))
        !
        IF(ANY(w(itetra(1:4)) < 0.0_dp)) THEN
           WRITE(*,'(100e15.5)') e(1:4)
           WRITE(*,'(100e15.5)') w(itetra(1:4))
           CALL errore("dfpt_tetra_lindhard", "Case : 4=3, 2=1", 1)
        END IF
        !
     ELSE
        !
        ! e(4) = e(3)
        !
        w(itetra(4)) = dfpt_tetra_lindhard_1231(e(4),e(1),e(2),le(4),le(1),le(2))
        w(itetra(3)) = w(itetra(4))
        w(itetra(2)) = dfpt_tetra_lindhard_1233(e(2),e(1),e(4),le(2),le(1),le(4))
        w(itetra(1)) = dfpt_tetra_lindhard_1233(e(1),e(2),e(4),le(1),le(2),le(4))
        !
        IF(ANY(w(itetra(1:4)) < 0.0_dp)) THEN
           WRITE(*,'(100e15.5)') e(1:4)
           WRITE(*,'(100e15.5)') w(itetra(1:4))
           CALL errore("dfpt_tetra_lindhard", "Case : 4=3", 1)
        END IF
        !
     END IF
  ELSE IF(ABS(e(3) - e(2)) < thr) THEN
     IF(ABS(e(3) - e(1)) < thr) THEN
        !
        ! e(3) = e(2) = e(1)
        !
        w(itetra(4)) = dfpt_tetra_lindhard_1222(e(4),e(3), le(4),le(3))
        w(itetra(3)) = dfpt_tetra_lindhard_1211(e(3),e(4), le(3),le(4))
        w(itetra(2)) = w(itetra(3))
        w(itetra(1)) = w(itetra(3))
        !
        IF(ANY(w(itetra(1:4)) < 0.0_dp)) THEN
           WRITE(*,'(100e15.5)') e(1:4)
           WRITE(*,'(100e15.5)') w(itetra(1:4))
           CALL errore("dfpt_tetra_lindhard", "Case : 3=2=1", 1)
        END IF
        !
     ELSE
        !
        ! e(3) = e(2)
        !
        w(itetra(4)) = dfpt_tetra_lindhard_1233(e(4),e(1),e(3),le(4),le(1),le(3))
        w(itetra(3)) = dfpt_tetra_lindhard_1231(e(3),e(1),e(4),le(3),le(1),le(4))
        w(itetra(2)) = w(itetra(3))
        w(itetra(1)) = dfpt_tetra_lindhard_1233(e(1),e(4),e(3),le(1),le(4),le(3))
        !
        IF(ANY(w(itetra(1:4)) < 0.0_dp)) THEN
           WRITE(*,'(100e15.5)') e(1:4)
           WRITE(*,'(100e15.5)') w(itetra(1:4))
           CALL errore("dfpt_tetra_lindhard", "Case : 3=2", 1)
        END IF
        !
     END IF
  ELSE IF(ABS(e(2) - e(1)) < thr) THEN
     !
     ! e(2) = e(1)
     !
     w(itetra(4)) = dfpt_tetra_lindhard_1233(e(4),e(3),e(2),le(4),le(3),le(2))
     w(itetra(3)) = dfpt_tetra_lindhard_1233(e(3),e(4),e(2),le(3),le(4),le(2))
     w(itetra(2)) = dfpt_tetra_lindhard_1231(e(2),e(3),e(4),le(2),le(3),le(4))
     w(itetra(1)) = w(itetra(2))
     !
     IF(ANY(w(itetra(1:4)) < 0.0_dp)) THEN
        WRITE(*,'(100e15.5)') e(1:4)
        WRITE(*,'(100e15.5)') w(itetra(1:4))
        CALL errore("dfpt_tetra_lindhard", "Case : 2=1", 1)
     END IF
     !
  ELSE
     !
     ! Different each other.
     !
     w(itetra(4)) = dfpt_tetra_lindhard_1234(e(4),e(1),e(2),e(3),le(4),le(1),le(2),le(3))
     w(itetra(3)) = dfpt_tetra_lindhard_1234(e(3),e(1),e(2),e(4),le(3),le(1),le(2),le(4))
     w(itetra(2)) = dfpt_tetra_lindhard_1234(e(2),e(1),e(3),e(4),le(2),le(1),le(3),le(4))
     w(itetra(1)) = dfpt_tetra_lindhard_1234(e(1),e(2),e(3),e(4),le(1),le(2),le(3),le(4))
     !      
     IF(ANY(w(itetra(1:4)) < 0.0_dp)) THEN
        WRITE(*,'(100e15.5)') e(1:4)
        WRITE(*,'(100e15.5)') w(itetra(1:4))
        CALL errore("dfpt_tetra_lindhard", "Case : Different each other", 1)
     END IF
     !
  END IF
  !
  nesting_factor = SUM(w(1:4)*pi(1:4)*pj(1:4))
  !
END SUBROUTINE nesting_proj_lindhard
END MODULE nesting_proj_routines
!
!----------------------------------------------------------------------------
PROGRAM nesting_proj
  !----------------------------------------------------------------------------
  !
  ! Usage :
  ! $ nesting_proj.x -in {input file}
  ! Then it generates proj.frmsf (for nspin = 1, 4) or
  ! proj1.frmsf and proj2.frmsf (for nspin = 2)
  !
  ! Input file format (projwfc.x + tail):
  ! &PROJWFC
  ! prefix = "..."
  ! outdir = "..."
  ! ...
  ! /
  ! {Number of target WFCs}
  ! {Index of WFC1} {Index of WFC2} {Index of WFC3} ...
  !
  USE parameters,           ONLY : npk
  USE cell_base,            ONLY : at, bg
  USE io_files,             ONLY : prefix, tmp_dir
  USE mp_global,            ONLY : mp_startup, mp_global_end
  USE environment,          ONLY : environment_start, environment_end
  USE kinds,                ONLY : DP
  USE wvfct,                ONLY : nbnd, et
  USE start_k,              ONLY : nk1, nk2, nk3
  USE lsda_mod,             ONLY : nspin
  USE ener,                 ONLY : ef, ef_up, ef_dw
  USE klist,                ONLY : nks, two_fermi_energies, xk, nkstot, wk
  USE fermisurfer_common,   ONLY : b_low, b_high, rotate_k_fs, write_fermisurfer
  USE nesting_proj_routines,  ONLY : read_projwfc_nesting, read_atomic_proj_nesting, nesting_proj_delta1, nesting_proj_theta1
  USE ktetra,               ONLY : opt_tetra_init, tetra_type, opt_tetra_partialdos
  !
  IMPLICIT NONE
  !
  INTEGER :: ik, i1, i2, i3, ns, nk, nbnd_fs, nirr_k, ikqv(3), s_dummy(3,3,48), t_rev_dummy(48)
  REAL(DP) :: ef1, ef2, dosf(2), pdosf(2), edos
  INTEGER,ALLOCATABLE :: equiv(:,:,:), irr_k(:,:)
  REAL(DP),ALLOCATABLE :: eig1(:,:,:,:,:), wfc1(:,:,:,:,:,:), eig2(:,:,:,:,:,:), wfc2(:,:,:,:,:,:), wt(:,:,:), &
  &                       nesting_irr(:,:), nesting_full(:,:,:,:)
  LOGICAL :: lbinary_data, needwf = .FALSE.
  !
  CHARACTER(LEN=256), EXTERNAL :: trimcheck
  !
  CALL mp_startup ()
  CALL environment_start ('NESTING_PROJ')
  !
  ! ... Read projwfc.x input file and get prefix, outdir, lbinary_data
  !
  CALL read_projwfc_nesting(lbinary_data)
  !
  ! ... Read XML file generated by pw.x
  !
  CALL read_file_new ( needwf )
  !
  ! ... Find equivalent k point in irr-BZ for whole BZ
  !
  ALLOCATE(equiv(nk1, nk2, nk3))
  CALL rotate_k_fs(equiv)
  !
  ef1 = 0.0_dp
  ef2 = 0.0_dp
  IF (nspin == 2) THEN
     ns = 2
     IF(two_fermi_energies) THEN
        ef1 = ef_up
        ef2 = ef_dw
     ELSE
        ef1 = ef
        ef2 = ef
     END IF
  ELSE
     ns = 1
  END IF
  nk = nks / ns
  !
  ! ... Read {prefix}.save/atomic_proj.* generated by projwfc.x
  !
  ALLOCATE(wt(nbnd,nks,2))
  CALL read_atomic_proj_nesting(lbinary_data, wt, ns, nk)
  !
  ! k in irreducible BZ
  !
  nirr_k = nk
  nks = nk1 * nk2 * nk3
  nkstot = nks
  ALLOCATE(nesting_irr(nirr_k,ns), nesting_full(nk1,nk2,nk3,ns), irr_k(3,nirr_k))
  !
  DO ik = 1, nirr_k
     irr_k(1:3,ik) = NINT(MATMUL(xk(1:3,ik), at(1:3,1:3)) * REAL((/nk1, nk2, nk3/), kind=dp))
  END DO
  !
  ik = 0
  DO i1 = 1, nk1
     DO i2 = 1, nk2
        DO i3 = 1, nk3
           ik = ik + 1
           xk(1:3,ik) = REAL((/i1, i2, i3/), DP) / REAL((/nk1, nk2, nk3/), DP)
           WHERE((/i1, i2, i3/)*2 >= (/nk1, nk2, nk3/)) xk(1:3,ik) = xk(1:3,ik) - 1.0_dp
           xk(1:3,ik) = MATMUL(bg(1:3,1:3), xk(1:3,ik))
           wk(ik) = 1.0_dp / REAL(nks, DP)
        END DO
     END DO
  END DO
  !
  ! Use tetrahedron method without symmetry
  !
  tetra_type = 2
  s_dummy(1:3,1:3,1:48) = 0
  t_rev_dummy(1:48) = 0
  DO i1 = 1, 3
     s_dummy(i1,i1,1:48) = 1
  END DO
  CALL opt_tetra_init(1, s_dummy, .False., t_rev_dummy, at, bg, npk, &
  &                   0, 0, 0, nk1, nk2, nk3, nks, xk, 1)
  !
  ! Nesting function function delta(e_k)*delta(e_{k+q})
  !
  nbnd_fs = b_high - b_low + 1
  ALLOCATE(eig1(nbnd_fs, nk3, nk2, nk1, ns), &
  &        eig2(nbnd_fs, nirr_k, nk3, nk2, nk1, ns))
  ALLOCATE(wfc1(nbnd_fs, nk3, nk2, nk1, ns, 2), &
  &        wfc2(nbnd_fs, nirr_k, nk3, nk2, nk1, ns))
  !
  ! ... Map e_k(Measured from E_F) and projected WFCs into whole BZ 
  !
  DO i1 = 1, nk1
     DO i2 = 1, nk2
        DO i3 = 1, nk3
           !
           IF(nspin == 2) THEN
              eig1(1:nbnd_fs,i3,i2,i1,1     ) = et(b_low:b_high, equiv(i1,i2,i3)     ) - ef1
              eig1(1:nbnd_fs,i3,i2,i1,2     ) = et(b_low:b_high, equiv(i1,i2,i3) + nk) - ef2
              wfc1(1:nbnd_fs,i3,i2,i1,1, 1:2) = wt(b_low:b_high, equiv(i1,i2,i3)     , 1:2)
              wfc1(1:nbnd_fs,i3,i2,i1,2, 1:2) = wt(b_low:b_high, equiv(i1,i2,i3) + nk, 1:2)
           ELSE
              eig1(1:nbnd_fs,i3,i2,i1,1     ) = et(b_low:b_high, equiv(i1,i2,i3)     ) - ef
              wfc1(1:nbnd_fs,i3,i2,i1,1, 1:2) = wt(b_low:b_high, equiv(i1,i2,i3)     , 1:2)
           END IF
           !
        END DO ! i3 = 1, nk3
     END DO ! i2 = 1, nk2
  END DO ! i1 = 1, nk1
  !
  ! Compute E(k+q) by shifting grid
  !
  DO ik = 1, nirr_k
     DO i1 = 1, nk1
        DO i2 = 1, nk2
           DO i3 = 1, nk3
              ikqv(1:3) = (/i1, i2, i3/) - 1 + irr_k(1:3,ik)
              ikqv(1:3) = MODULO(ikqv(1:3), (/nk1, nk2, nk3/)) + 1
              eig2(1:nbnd_fs,ik,i3,i2,i1,1:ns) = eig1(1:nbnd_fs,ikqv(3),ikqv(2),ikqv(1),1:ns)
              wfc2(1:nbnd_fs,ik,i3,i2,i1,1:ns) = wfc1(1:nbnd_fs,ikqv(3),ikqv(2),ikqv(1),1:ns,2)
           END DO
        END DO
     END DO
  END DO
  !
  nesting_irr(1:nirr_k,1:ns) = 0.0_dp
  !
  CALL nesting_proj_delta1(nbnd_fs,ns,nirr_k,irr_k,eig1,eig2,wfc1(:,:,:,:,:,1),wfc2,nesting_irr)
  !
  ! Extend nesting factor to the full BZ
  !
  DO i3 = 1, nk3
     DO i2 = 1, nk2
        DO i1 = 1, nk1
           nesting_full(i1,i2,i3,1:ns) = nesting_irr(equiv(i1,i2,i3), 1:ns)
        END DO
     END DO
  END DO
  !
  ! ... Output in the FermiSurfer format
  !
  b_low = 1
  b_high = 1
  IF (nspin == 2) THEN
     !
     CALL write_fermisurfer(nesting_full(1:nk1, 1:nk2, 1:nk3, 1), &
     &                      nesting_full(1:nk1, 1:nk2, 1:nk3, 1), TRIM(prefix) // "_nesting_proj_delta1.frmsf")
     CALL write_fermisurfer(nesting_full(1:nk1, 1:nk2, 1:nk3, 2), &
     &                      nesting_full(1:nk1, 1:nk2, 1:nk3, 2), TRIM(prefix) // "nesting_proj_delta2.frmsf")
     !
  ELSE
     !
     CALL write_fermisurfer(nesting_full(1:nk1, 1:nk2, 1:nk3, 1), &
     &                      nesting_full(1:nk1, 1:nk2, 1:nk3, 1), TRIM(prefix) // "_nesting_proj_delta.frmsf")
     !
  END IF
  !
  ! Nesting function with (f_k - f_{k+q})/(e_{k+q} - e_k)
  !
  nesting_irr(1:nirr_k,1:ns) = 0.0_dp 
  CALL nesting_proj_theta1(nbnd_fs,ns,nirr_k,eig1,eig2,wfc1(:,:,:,:,:,1),wfc2,nesting_irr)
  nesting_irr(1:nirr_k, 1:ns) = 2.0_dp * nesting_irr(1:nirr_k, 1:ns)
  IF (nspin == 1) nesting_irr(1:nirr_k, 1:ns) = 2.0_dp * nesting_irr(1:nirr_k, 1:ns)
  !
  ! Drude term
  !
  
  DEALLOCATE(et)
  ALLOCATE(et(nbnd_fs,nk1*nk2*nk3*ns))
  et(1:nbnd_fs,1:nk1*nk2*nk3*ns) = RESHAPE(eig1, (/nbnd_fs, nk1*nk2*nk3*ns/))
  nbnd = nbnd_fs
  CALL opt_tetra_partialdos( ns, .FALSE., 0, 1, 1, &
  &                          edos, 1.0_dp, wfc1(:,:,:,:,:,1), pdosf, dosf, ns )
  !
  DO ik = 1, nirr_k
    IF(ALL(irr_k(1:3,ik) == 0)) THEN
      nesting_irr(ik, 1:ns) = nesting_irr(ik, 1:ns) + pdosf(1:ns)
    END IF
  END DO
  !
  ! Extend nesting factor to the full BZ
  !
  DO i3 = 1, nk3
     DO i2 = 1, nk2
        DO i1 = 1, nk1
           nesting_full(  i1,i2,i3,1:ns) = nesting_irr(equiv(i1,i2,i3), 1:ns)
        END DO
     END DO
  END DO
  !
  ! ... Output in the FermiSurfer format
  !
  b_low = 1
  b_high = 1
  IF (nspin == 2) THEN
     !
     CALL write_fermisurfer(nesting_full(1:nk1, 1:nk2, 1:nk3, 1), &
     &                      nesting_full(1:nk1, 1:nk2, 1:nk3, 1), &
     &                      TRIM(tmp_dir) // TRIM(prefix) // "_nesting_proj_chi1.frmsf")
     CALL write_fermisurfer(nesting_full(1:nk1, 1:nk2, 1:nk3, 2), &
     &                      nesting_full(1:nk1, 1:nk2, 1:nk3, 2), &
     &                      TRIM(tmp_dir) // TRIM(prefix) // "_nesting_proj_chi2.frmsf")
     !
  ELSE
     !
     CALL write_fermisurfer(nesting_full(1:nk1, 1:nk2, 1:nk3, 1), &
     &                      nesting_full(1:nk1, 1:nk2, 1:nk3, 1), &
     &                      TRIM(tmp_dir) // TRIM(prefix) // "_nesting_proj_chi.frmsf")
     !
  END IF
  !
  DEALLOCATE(eig1, eig2, wfc1, wfc2, nesting_full, nesting_irr, equiv)
  !

  CALL environment_end ('NESTING_PROJ')
#if defined(__MPI)
  CALL mp_global_end ( )
#endif
  !
END PROGRAM nesting_proj
