Line data    Source code

       1             : MODULE m_dfpt_fermie
       2             :    USE m_juDFT
       3             : #ifdef CPP_MPI
       4             :    USE mpi
       5             : #endif
       6             : 
       7             : CONTAINS
       8           0 :    SUBROUTINE dfpt_fermie(eig_id,dfpt_eig_id,fmpi,kpts,input,noco,results,results1)
       9             :       !! Calculate the perturbed occupation numbers from the unperturbed ones and the
      10             :       !! perturbed eigenenergies.
      11             :       !! This is only done for metals, i.e. systems where the smearing is not set
      12             :       !! to 0.
      13             :       !! Fermi-Dirac smearing is assumed.
      14             :       USE m_types
      15             :       USE m_constants
      16             :       USE m_eig66_io, ONLY : read_eig, write_eig
      17             : 
      18             :       IMPLICIT NONE
      19             : 
      20             :       TYPE(t_results), INTENT(INOUT) :: results, results1
      21             :       TYPE(t_mpi),     INTENT(IN)    :: fmpi
      22             :       TYPE(t_input),   INTENT(IN)    :: input
      23             :       TYPE(t_noco),    INTENT(IN)    :: noco
      24             :       TYPE(t_kpts),    INTENT(IN)    :: kpts
      25             : 
      26             :       INTEGER, INTENT(IN) :: eig_id, dfpt_eig_id
      27             : 
      28             :       REAL    :: efermi, ef_num, ef_den, x
      29             :       INTEGER :: j, jsp, k, nspins, noccbd
      30             : 
      31           0 :       REAL, ALLOCATABLE :: sxm(:,:,:)
      32             : 
      33             : #ifdef CPP_MPI
      34             :       INTEGER, PARAMETER :: comm = MPI_COMM_SELF
      35             :       INTEGER ierr
      36             : #endif
      37             : 
      38           0 :       IF (noco%l_noco) THEN
      39             :          nspins = 1
      40             :       ELSE
      41           0 :          nspins = input%jspins
      42             :       END IF
      43             : 
      44           0 :       ALLOCATE(sxm(MAXVAL(results%neig),kpts%nkpt,nspins))
      45             : 
      46           0 :       IF (fmpi%irank == 0) THEN
      47           0 :          efermi = results%ef
      48           0 :          results1%ef = 0.0
      49           0 :          ef_num = 0.0
      50           0 :          ef_den = 0.0
      51             : 
      52           0 :          DO jsp = 1, nspins
      53           0 :             DO k = 1, kpts%nkpt
      54           0 :                noccbd  = COUNT(results%w_iks(:,k,jsp)*2.0/input%jspins>1.e-8)
      55           0 :                DO j = 1, noccbd
      56           0 :                   x = (results%eig(j,k,jsp)-efermi)/input%tkb
      57           0 :                   sxm(j,k,jsp) = sfermi(-x)
      58           0 :                   ef_num = ef_num + results%w_iks(j,k,jsp) * sxm(j,k,jsp) * results1%eig(j,k,jsp)
      59           0 :                   ef_den = ef_den + results%w_iks(j,k,jsp) * sxm(j,k,jsp)
      60             :                END DO
      61             :             END DO
      62             :          END DO
      63             : 
      64           0 :          IF (ABS(ef_den)>1e-12) THEN
      65           0 :             results1%ef = ef_num/ef_den
      66             :          ELSE
      67             :             results1%ef = 0.0
      68             :          END IF
      69             : 
      70             :          results1%w_iks(:noccbd,:,1:nspins) = -results%w_iks(:noccbd,:,1:nspins) &
      71             :                                             * sxm(:noccbd,:,1:nspins) &
      72           0 :                                             * (results1%eig(:noccbd,:,1:nspins)-results1%ef)/input%tkb
      73             :       END IF
      74             : 
      75           0 :       RETURN
      76           0 :    END SUBROUTINE dfpt_fermie
      77             : 
      78           0 :    REAL FUNCTION sfermi(x)
      79             :       !! Returns the Fermi-Dirac function
      80             :       !! $$s(x)=(e^{x}+1)^{-1}$$
      81             :       !! for \(x=(\epsilon_{\nu\boldsymbol{k}}-E_{F})/(k_{B}T)\).
      82             : 
      83             :       REAL, INTENT(IN) :: x
      84             : 
      85             :       REAL :: expo
      86             : 
      87           0 :       expo = EXP(x)
      88             : 
      89           0 :       sfermi = 1.0/(expo+1.0)
      90             : 
      91             :       RETURN
      92             : 
      93             :    END FUNCTION sfermi
      94             : 
      95             : END MODULE m_dfpt_fermie