Line data    Source code

       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : MODULE m_ptdos
       7             :    !-------------------------------------------------------------------------
       8             :    ! Density of states calculated by linear triangular method
       9             :    !-------------------------------------------------------------------------
      10             :    USE m_types_kpts
      11             :    USE m_tetsrt
      12             : 
      13             :    IMPLICIT NONE
      14             : 
      15             :    CONTAINS
      16             : 
      17           0 :    SUBROUTINE ptdos(jspins,kpts,eGrid,eig,qal,g,energyShift)
      18             : 
      19             :       INTEGER,       INTENT(IN)  :: jspins
      20             :       TYPE(t_kpts),  INTENT(IN)  :: kpts
      21             :       REAL,          INTENT(IN)  :: qal(:,:,:)  !(nbands,nkpt,jspins)
      22             :       REAL,          INTENT(IN)  :: eGrid(:)
      23             :       REAL,          INTENT(IN)  :: eig(:,:,:)   !(nbands,nkpt,jspins)
      24             :       REAL,          INTENT(OUT) :: g(:,:)      !(ne,jspins)
      25             :       REAL, OPTIONAL, INTENT(IN) :: energyShift
      26             : 
      27             :       INTEGER :: iBand,itria,iGrid,ispin
      28             :       INTEGER :: ind(3),k(3)
      29             :       REAL    :: sfac,fa,shift
      30             :       REAL    :: ei(3)
      31             : 
      32           0 :       shift = 0.0
      33           0 :       IF(PRESENT(energyShift)) shift = energyShift
      34             : 
      35             :       !Spin-degeneracy factor
      36           0 :       sfac = 2.0/jspins
      37             : 
      38           0 :       g = 0.0
      39           0 :       DO ispin = 1, SIZE(qal,3)
      40           0 :          DO itria = 1 , kpts%ntet
      41           0 :             fa = sfac*kpts%voltet(itria)/kpts%ntet
      42           0 :             k = kpts%ntetra(:,itria)
      43           0 :             DO iBand = 1 , SIZE(eig,1)
      44             :                !---------------------------
      45             :                !eigenvalues at the corners
      46             :                !of the current triangle
      47             :                !---------------------------
      48           0 :                ei = eig(iBand,k,ispin) - shift
      49             : 
      50             :                !sort in ascending order
      51           0 :                ind = tetsrt(ei)
      52             : 
      53             :                !
      54             :                !calculate partial densities of states
      55             :                !
      56             : 
      57           0 :                DO iGrid = 1 , SIZE(eGrid)
      58             :                   g(iGrid,ispin) = g(iGrid,ispin) + fa &
      59           0 :                                   * dostet( eGrid(iGrid),ei(ind),qal(iBand,k(ind),ispin) )
      60             :                ENDDO
      61             : 
      62             :             ENDDO
      63             :          ENDDO
      64             :       ENDDO
      65             : 
      66           0 :    END SUBROUTINE ptdos
      67             : 
      68           0 :    REAL FUNCTION dostet(e,ei,q)
      69             :       !
      70             :       !     partial density of states for one tetrahedron
      71             :       !     note that e1.le.e2.le.e3 is assumed
      72             :       !
      73             :       !     ei : eigenvalues at the corners of the triangle (sorted)
      74             :       !     q : partial charges at the corners
      75             : 
      76             :       REAL, INTENT(IN) :: e      !Energy point where the DOS is calculated
      77             :       REAL, INTENT(IN) :: ei(:)  !(3)
      78             :       REAL, INTENT(IN) :: q(:)   !(3)
      79             : 
      80             :       REAL :: e21 , e31 , e32 , ee
      81             :       REAL, PARAMETER :: tol = 1e-6
      82             : 
      83           0 :       dostet = 0.0
      84           0 :       IF ( e.LT.ei(1) ) RETURN
      85           0 :       IF ( e.LE.ei(2) ) THEN
      86             :          !
      87             :          !     case 1: e between e1 and e2
      88             :          !
      89           0 :          ee = e - ei(1)
      90           0 :          e21 = ei(2) - ei(1)
      91           0 :          IF ( e21.LT.tol ) RETURN
      92           0 :          e31 = ei(3) - ei(1)
      93           0 :          IF ( e31.LT.tol ) RETURN
      94             :          dostet = ee/(e21*e31)*(2.0*q(1) &
      95             :                                 + (ee/e21*(q(2)-q(1)) &
      96           0 :                                 + ee/e31*(q(3)-q(1))))
      97             :       ELSE
      98           0 :          IF ( e.GT.ei(3) ) RETURN
      99             :          !
     100             :          !     case 2: e between e2 and e3
     101             :          !
     102           0 :          e31 = ei(3) - ei(1)
     103           0 :          IF ( e31.LT.tol ) RETURN
     104           0 :          e32 = ei(3) - ei(2)
     105           0 :          IF ( e32.LT.tol ) RETURN
     106             :          dostet = (ei(3)-e)/(e31*e32)*(q(1)+q(2) &
     107             :                                         +(e-ei(1))/e31*(q(3)-q(1)) &
     108           0 :                                         +(e-ei(2))/e32*(q(3)-q(2)))
     109           0 :          RETURN
     110             :       ENDIF
     111           0 :    END FUNCTION dostet
     112             : 
     113             : END MODULE m_ptdos