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_doswt
       7             :    !
       8             :    !     calculates the weights for each k-point for integrating functions
       9             :    !     of k.  the array w has beeen cleared before entering.
      10             :    !
      11             :    USE m_trisrt
      12             :    USE m_types
      13             : 
      14             :    IMPLICIT NONE
      15             : 
      16             :    CONTAINS
      17             : 
      18           0 :    SUBROUTINE doswt(ei,nemax,jspins,kpts,eig,w)
      19             : 
      20             :       INTEGER,       INTENT(IN) :: jspins
      21             :       TYPE(t_kpts),  INTENT(IN) :: kpts
      22             :       REAL,          INTENT(IN) :: ei
      23             :       INTEGER,       INTENT(IN) :: nemax(:)
      24             :       REAL,          INTENT(IN) :: eig(:,:,:)      !(neig,nkpt,jspins)
      25             :       REAL,          INTENT(OUT):: w(:,:,:)        !(neig,nkpt,jspins)
      26             : 
      27             :       INTEGER :: jspin,iBand,itria
      28             :       INTEGER :: k1,k2,k3
      29             :       INTEGER :: neig
      30             :       REAL    :: e1,e2,e3
      31             :       REAl    :: ee,e32,e31,e21,s
      32           0 :       w=0.0 !init was missing
      33           0 :       DO jspin = 1,jspins
      34           0 :          neig = nemax(jspin)
      35           0 :          DO iBand = 1,neig
      36           0 :             DO itria = 1,kpts%ntet
      37             :                !Get the k-points and eigenvalues
      38             :                !of the current triangle
      39           0 :                k1 = kpts%ntetra(1,itria)
      40           0 :                k2 = kpts%ntetra(2,itria)
      41           0 :                k3 = kpts%ntetra(3,itria)
      42           0 :                e1 = eig(iBand,k1,jspin)
      43           0 :                e2 = eig(iBand,k2,jspin)
      44           0 :                e3 = eig(iBand,k3,jspin)
      45             :                !Sort by ascending eigenvalues
      46           0 :                CALL trisrt(e1,e2,e3,k1,k2,k3)
      47           0 :                IF (e1.LE.-9999.0) CYCLE !Not all eigenvalues available
      48           0 :                IF (ei.LE.e1) CYCLE !triangle not occupied
      49           0 :                IF (ei.GE.e3) THEN
      50             :                   !---> e3<e
      51           0 :                   s = kpts%voltet(itria)/kpts%ntet/3.0
      52           0 :                   w(iBand,k1,jspin) = w(iBand,k1,jspin) + s
      53           0 :                   w(iBand,k2,jspin) = w(iBand,k2,jspin) + s
      54           0 :                   w(iBand,k3,jspin) = w(iBand,k3,jspin) + s
      55           0 :                ELSEIF (ei.GT.e2) THEN
      56             :                   !---> e2<ei<e3
      57           0 :                   ee = e3 - ei
      58           0 :                   e31 = ee/ (e3-e1)
      59           0 :                   e32 = ee/ (e3-e2)
      60           0 :                   s = kpts%voltet(itria)/kpts%ntet/3.0
      61           0 :                   w(iBand,k1,jspin) = w(iBand,k1,jspin) + s* (1.-e31*e31*e32)
      62           0 :                   w(iBand,k2,jspin) = w(iBand,k2,jspin) + s* (1.-e31*e32*e32)
      63           0 :                   w(iBand,k3,jspin) = w(iBand,k3,jspin) + s* (1.-e31*e32*(3.-e31-e32))
      64             :                ELSE
      65             :                   !---> e1<ei<e2
      66           0 :                   ee = ei - e1
      67           0 :                   e31 = ee/ (e3-e1)
      68           0 :                   e21 = ee/ (e2-e1)
      69           0 :                   s = kpts%voltet(itria)/kpts%ntet*e31*e21/3.0
      70           0 :                   w(iBand,k1,jspin) = w(iBand,k1,jspin) + s* (3.0-e21-e31)
      71           0 :                   w(iBand,k2,jspin) = w(iBand,k2,jspin) + s*e21
      72           0 :                   w(iBand,k3,jspin) = w(iBand,k3,jspin) + s*e31
      73             :                ENDIF
      74             : 
      75             :             ENDDO
      76             :          ENDDO
      77             :       ENDDO
      78             : 
      79           0 :    END SUBROUTINE doswt
      80             : END MODULE m_doswt