xcallg.f 7.11 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8 9 10 11 12 13 14 15 16 17 18
      MODULE m_xcallg
      use m_juDFT

c  vxcallg: driver subroutine for the exchange-correlation potentials
c  excallg: driver subroutine for the exchange-correlation energy density.
c
c  inputs a total bare electron density rh(r,1) for jspins=1
c         a bare electron density rh(r,1),rh(r,2) of majority and
c          minority for jspins=2.
c  outputs the effective exch-corr energy density due to this electron
c  density in hartree units.
c
19
c  depending on the variable xcpot
20 21 22 23 24 25 26 27 28 29 30 31 32 33
c  the contribution of exchange and correlation is determined in
c  different subroutines.
c
c
c     krla=1 means : relativistic correction of exchange energy
c                    density related to dirac kinetic energy, according
c                    to: a.h. macdonald and s.h. vosko, j. phys. c12,
c                    2977(1979)
c
c     based on a subroutine from s.bluegel
c     r.pentcheva, 22.01.96

      CONTAINS
c***********************************************************************
34
      SUBROUTINE vxcallg(xcpot,lwbc,jspins,mfftwk,nfftwk,rh,
35 36 37 38 39 40 41 42 43
     +                  agr,agru,agrd,g2r,g2ru,g2rd,gggr,gggru,gggrd,
     +                  gzgr,
     +                  vx,vxc)
c***********************************************************************
c
      USE m_vxcl91
      USE m_vxcwb91
      USE m_vxcpw91
      USE m_vxcepbe
44
      USE m_types
45 46 47 48
      IMPLICIT NONE
c
c---> running mode parameters
c
49
      TYPE(t_xcpot),INTENT(IN) :: xcpot
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
      INTEGER, INTENT (IN) :: mfftwk,nfftwk ! radial mesh,number of mesh points
      INTEGER, INTENT (IN) :: jspins
      LOGICAL, INTENT (IN) :: lwbc          ! l-white-bird-current (ta)
c
c---> charge density
c
      REAL,INTENT (IN) :: agr(mfftwk)
      REAL,INTENT (IN) :: rh(mfftwk,jspins)
      REAL,INTENT (IN) :: agru(mfftwk),agrd(mfftwk)
      REAL,INTENT (IN) :: g2r(mfftwk),g2ru(mfftwk),g2rd(mfftwk)
      REAL,INTENT (IN) :: gggr(mfftwk),gggru(mfftwk)
      REAL, INTENT (IN) :: gggrd(mfftwk),gzgr(mfftwk)
c
c---> xc potential
c
      REAL, INTENT (OUT) :: vx (mfftwk,jspins)
      REAL, INTENT (OUT) :: vxc(mfftwk,jspins)
c
c ---> local scalars
      INTEGER ir,js
      REAL, PARAMETER :: hrtr_half = 0.5e0

      !used to be dummy arguments for testing
      INTEGER,PARAMETER   :: idsprs=0,isprsv=0
      REAL,PARAMETER      :: sprsv=0.0
c
c.....------------------------------------------------------------------
c
c-----> determine exchange correlation potential
c
      vx (:,:) = 0.0
      vxc(:,:) = 0.0

83
      IF (xcpot%is_name("l91")) THEN    ! local pw91
84 85 86 87 88 89 90

       CALL vxcl91(
     >             jspins,mfftwk,nfftwk,rh,agr,agru,agrd,
     >             g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr,
     <             vx,vxc,
     >             isprsv,sprsv)

91
      ELSEIF (xcpot%is_name("pw91")) THEN  ! pw91
92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108

       IF (lwbc) THEN
          CALL vxcwb91(
     >                 jspins,mfftwk,nfftwk,rh,agr,agru,agrd,
     >                 g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr,
     <                 vx,vxc,
     >                 idsprs,isprsv,sprsv)
       ELSE

          CALL vxcpw91(
     >                 jspins,mfftwk,nfftwk,rh,agr,agru,agrd,
     >                 g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr,
     <                 vx,vxc,
     >                 idsprs,isprsv,sprsv)

       ENDIF

109
      ELSE  ! pbe or similar
110 111

        CALL vxcepbe(
112
     >               xcpot,jspins,mfftwk,nfftwk,rh,
113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
     >               agr,agru,agrd,g2ru,g2rd,gggr,gggru,gggrd,
     <               vx,vxc)


      ENDIF

c
c-----> hartree units
c
      IF (jspins.EQ.2) THEN
          DO ir = 1,nfftwk
              vx(ir,1)      = hrtr_half*vx(ir,1)
              vx(ir,jspins) = hrtr_half*vx(ir,jspins)        
          
              vxc(ir,1)      = hrtr_half*vxc(ir,1)
              vxc(ir,jspins) = hrtr_half*vxc(ir,jspins)
          ENDDO
      ELSEIF (jspins.eq.1) THEN
          DO ir = 1,nfftwk
              vx(ir,1)  = hrtr_half*vx(ir,1)
              
              vxc(ir,1) = hrtr_half*vxc(ir,1)
          ENDDO
      ELSE
          WRITE (6,fmt='('' error in jspins, jspins ='',i2)')
     +      jspins
          CALL juDFT_error("error in jspins",calledby ="xcallg")
      ENDIF

      END SUBROUTINE vxcallg
!*********************************************************************
      SUBROUTINE excallg(
145
     >                  xcpot,lwbc,jspins,nfftwk,
146 147 148 149 150 151 152 153 154
     >                  rh,agr,agru,agrd,
     >                  g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr,
     <                  exc)
c***********************************************************************
c
      USE m_excl91
      USE m_excwb91
      USE m_excpw91
      USE m_excepbe
155 156
      USE m_types
   
157 158 159 160
      IMPLICIT NONE
c
c---> running mode parameters
c
161
      TYPE(t_xcpot),INTENT(IN) :: xcpot
162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
      LOGICAL, INTENT (IN) :: lwbc          ! l-white-bird-current (ta)
      INTEGER, INTENT (IN) :: nfftwk ! radial mesh,number of mesh points
      INTEGER, INTENT (IN) :: jspins
c
c---> charge density  & gradients
c
      REAL,INTENT (IN) :: agr(:)
      REAL,INTENT (IN) :: rh(:,:)
      REAL,INTENT (IN) :: agru(size(agr)),agrd(size(agr))
      REAL,INTENT (IN) :: g2r(size(agr)),g2ru(size(agr)),g2rd(size(agr))
      REAL,INTENT (IN) :: gggr(size(agr)),gggru(size(agr))
      REAL,INTENT (IN) :: gggrd(size(agr)),gzgr(size(agr))
c
c---> xc energy density
c
      REAL, INTENT (OUT) :: exc(size(agr))
c
c ---> local scalars
      INTEGER ir
      REAL, PARAMETER :: hrtr_half = 0.5e0

      !used to be dummy arguments for testing
      INTEGER,PARAMETER   :: idsprs=0,isprsv=0
      REAL,PARAMETER      :: sprsv=0.0
c
c-----> determine exchange correlation energy density
c
c     write(6,'(/'' icorr,krla,igrd,jspins,lwbc,='',5i5,l2)')
c    &  icorr,krla,igrd,jspins,lwbc

192
      IF (xcpot%is_name("l91")) THEN  ! local pw91
193 194 195 196 197 198 199

        CALL excl91(
     >              jspins,size(rh,1),nfftwk,rh,agr,agru,agrd,
     >              g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr,
     <              exc,
     >              isprsv,sprsv)

200
      ELSEIF (xcpot%is_name("pw91")) THEN     ! pw91
201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218

        IF (lwbc) THEN
          CALL excwb91(
     >                 size(agr),nfftwk,
     >                 rh(:,1),rh(:,2),agr,agru,agrd,
     >                 g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr,
     <                 exc,
     >                 idsprs,isprsv,sprsv)
        ELSE

          CALL excpw91(
     >                 jspins,size(agr),nfftwk,rh,agr,agru,agrd,
     >                 g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr,
     <                 exc,
     >                 idsprs,isprsv,sprsv)

         ENDIF

219
      ELSE
220
        CALL excepbe(
221
     >               xcpot,jspins,size(rh,1),nfftwk,
222 223 224 225 226 227 228 229 230 231 232 233 234 235
     >               rh,agr,agru,agrd,g2ru,g2rd,gggr,gggru,gggrd,
     <               exc)

      ENDIF
c
c-----> hartree units
c
      DO ir = 1,nfftwk
        exc(ir) = hrtr_half*exc(ir)
      ENDDO


      END SUBROUTINE excallg
      END MODULE m_xcallg