types_xcpot_inbuild.F90 13.8 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.
!--------------------------------------------------------------------------------
MODULE m_types_xcpot_inbuild
7 8 9 10 11 12 13 14 15 16 17
   !This module contains the xcpot-type used for the in-build xc-implementations
   USE m_types_xcpot_data
   USE m_types_xcpot
   USE m_judft
   IMPLICIT NONE
   PRIVATE
   REAL, PARAMETER, PRIVATE :: hrtr_half = 0.5
   CHARACTER(len=4),PARAMETER:: xc_names(20)=[&
                                'l91 ','x-a ','wign','mjw ','hl  ','bh  ','vwn ','pz  ', &
                                'pw91','pbe ','rpbe','Rpbe','wc  ','PBEs', &
                                'pbe0','hse ','vhse','lhse','exx ','hf  ']
Matthias Redies's avatar
Matthias Redies committed
18 19 20 21 22 23
   
   LOGICAL,PARAMETER:: priv_LDA(20)=[&
                       .FALSE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,&
                       .FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,&
                       .FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.]

24 25 26 27

   LOGICAL,PARAMETER:: priv_gga(20)=[&
                       .TRUE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,&
                       .TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,&
28
                       .TRUE.,.TRUE.,.TRUE.,.TRUE.,.FALSE.,.TRUE.]
29 30 31 32 33 34 35 36 37 38 39 40 41

   LOGICAL,PARAMETER:: priv_hybrid(20)=[&
                       .FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,&
                       .FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,&
                       .TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.]

   REAL, PARAMETER       ::  amix_pbe0 = 0.25
   REAL, PARAMETER       ::  amix_hse  = 0.25
   REAL, PARAMETER       ::  amix_hf   = 1.00

   TYPE, EXTENDS(t_xcpot):: t_xcpot_inbuild
#ifdef CPP_MPI
      INTEGER             :: icorr=0 !not private to allow bcasting it around
42
#else
43
      INTEGER,PRIVATE     :: icorr=0
44 45
#endif

46 47 48 49
      TYPE(t_xcpot_data)   :: DATA

      LOGICAL,ALLOCATABLE :: lda_atom(:)

50
   CONTAINS
51
      !overloading t_xcpot:
Matthias Redies's avatar
Matthias Redies committed
52 53
      PROCEDURE        :: vxc_is_LDA=>xcpot_is_LDA
      PROCEDURE        :: exc_is_LDA=>xcpot_is_LDA
54 55
      PROCEDURE        :: vxc_is_gga=>xcpot_vxc_is_gga
      PROCEDURE        :: exc_is_gga=>xcpot_exc_is_gga
56 57 58 59 60 61 62 63 64 65
      PROCEDURE        :: is_hybrid=>xcpot_is_hybrid
      PROCEDURE        :: get_exchange_weight=>xcpot_get_exchange_weight
      PROCEDURE        :: get_vxc=>xcpot_get_vxc
      PROCEDURE        :: get_exc=>xcpot_get_exc
      !not overloaded
      PROCEDURE        :: get_name=>xcpot_get_name
      PROCEDURE        :: is_name=>xcpot_is_name
      PROCEDURE        :: init=>xcpot_init
   END TYPE t_xcpot_inbuild
   PUBLIC t_xcpot_inbuild
66
CONTAINS
67 68 69 70 71 72 73
   CHARACTER(len=4) FUNCTION xcpot_get_name(xcpot)
      USE m_judft
      IMPLICIT NONE
      CLASS(t_xcpot_inbuild),INTENT(IN)    :: xcpot
      IF (xcpot%icorr==0) CALL judft_error("xc-potential not initialized",calledby="types_xcpot.F90")
      xcpot_get_name=xc_names(xcpot%icorr)
   END FUNCTION xcpot_get_name
74

75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116
   SUBROUTINE xcpot_init(xcpot,namex,relcor,ntype)
      USE m_judft
      IMPLICIT NONE
      CLASS(t_xcpot_inbuild),INTENT(INOUT)    :: xcpot
      CHARACTER(len=*),INTENT(IN)  :: namex
      LOGICAL,INTENT(IN)           :: relcor
      INTEGER,INTENT(IN)           :: ntype
      INTEGER:: n
      !Determine icorr from name

      ALLOCATE(xcpot%lda_atom(ntype))
      xcpot%lda_atom=.FALSE.
      xcpot%icorr=0
      DO n=1,SIZE(xc_names)
         IF (TRIM(ADJUSTL(namex))==TRIM(xc_names(n))) THEN
            xcpot%icorr=n
         ENDIF
      ENDDO
      if (xcpot%icorr==0) CALL judft_error("Unkown xc-potential:"//namex,calledby="types_xcpot.F90")
      xcpot%data%krla=MERGE(1,0,relcor)

      !Code from exchpbe to speed up determination of constants
      IF (xcpot%is_name("rpbe")) THEN
         xcpot%data%uk=1.2450
      ELSE
         xcpot%data%uk=0.8040
      ENDIF
      IF (xcpot%is_name("PBEs")) THEN     ! pbe_sol
         xcpot%data%um=0.123456790123456d0
      ELSE
         xcpot%data%um=0.2195149727645171e0
      ENDIF
      xcpot%data%is_hse=xcpot%is_name("hse").OR.xcpot%is_name("lhse").OR.xcpot%is_name("vhse")
      xcpot%data%is_rpbe=xcpot%is_name("Rpbe") !Rpbe
      xcpot%data%is_wc=xcpot%is_name("wc")
      xcpot%data%is_pbes=xcpot%is_name("PBEs")
      xcpot%data%is_pbe0=xcpot%is_name("pbe0")
      xcpot%data%is_mjw=xcpot%is_name("mjw")
      xcpot%data%is_bh=xcpot%is_name("bh")
      xcpot%DATA%exchange_weight=xcpot%get_exchange_weight()

   END SUBROUTINE xcpot_init
Matthias Redies's avatar
Matthias Redies committed
117 118
  
   !! LDA
Matthias Redies's avatar
Matthias Redies committed
119 120
   
   LOGICAL FUNCTION xcpot_is_lda(xcpot)
Matthias Redies's avatar
Matthias Redies committed
121 122
      IMPLICIT NONE
      CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot
Matthias Redies's avatar
Matthias Redies committed
123
      xcpot_is_lda = (.not. xcpot%vxc_is_gga()) .and. (.not. xcpot%is_hybrid())
Matthias Redies's avatar
Matthias Redies committed
124
   END FUNCTION xcpot_is_lda
Matthias Redies's avatar
Matthias Redies committed
125
   
Matthias Redies's avatar
Matthias Redies committed
126
   !! GGA
127 128 129 130 131 132 133 134

   LOGICAL FUNCTION xcpot_vxc_is_gga(xcpot)
      IMPLICIT NONE
      CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot
      xcpot_vxc_is_gga=priv_gga(xcpot%icorr)
   END FUNCTION xcpot_vxc_is_gga

   LOGICAL FUNCTION xcpot_exc_is_gga(xcpot)
135 136
      IMPLICIT NONE
      CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot
137 138
      xcpot_exc_is_gga = xcpot%vxc_is_gga()
   END FUNCTION xcpot_exc_is_gga
139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160

   LOGICAL FUNCTION xcpot_is_hybrid(xcpot)
      IMPLICIT NONE
      CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot
      xcpot_is_hybrid=priv_hybrid(xcpot%icorr)
   END FUNCTION xcpot_is_hybrid

   FUNCTION xcpot_get_exchange_weight(xcpot) RESULT(a_ex)
      USE m_judft
      IMPLICIT NONE
      CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot

      REAL:: a_ex

      a_ex=-1
      IF (xcpot%is_name("pbe0")) a_ex=amix_pbe0
      IF (xcpot%is_name("hf")) a_ex=amix_hf
      IF (xcpot%is_name("hse")) a_ex=amix_hse
      IF (xcpot%is_name("vhse")) a_ex=amix_hse
   END FUNCTION xcpot_get_exchange_weight

   SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad)
161
!
162 163 164 165 166 167 168 169 170
      USE m_xcxal, ONLY : vxcxal
      USE m_xcwgn, ONLY : vxcwgn
      USE m_xcbh,  ONLY : vxcbh
      USE m_xcvwn, ONLY : vxcvwn
      USE m_xcpz,  ONLY : vxcpz
      USE m_vxcl91
      USE m_vxcwb91
      USE m_vxcpw91
      USE m_vxcepbe
171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187
      IMPLICIT NONE
!c
!c---> running mode parameters
!c
      CLASS(t_xcpot_inbuild),INTENT(IN) :: xcpot
      INTEGER, INTENT (IN)     :: jspins
!c
!c---> charge density
!c
      REAL,INTENT (IN) :: rh(:,:)
!c
!c---> xc potential
!c
      REAL, INTENT (OUT) :: vx (:,:)
      REAL, INTENT (OUT) :: vxc(:,:)

      ! optional arguments for GGA
188
      TYPE(t_gradients),INTENT(INOUT),OPTIONAL::grad
189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
!c
!c ---> local scalars
      INTEGER :: ngrid
      REAL, PARAMETER :: hrtr_half = 0.5

      !used to be dummy arguments for testing
      INTEGER,PARAMETER   :: idsprs=0,isprsv=0,iofile=6
      REAL,PARAMETER      :: sprsv=0.0
      LOGICAL,PARAMETER   :: lwbc=.false. ! l-white-bird-current (ta)
!c
!c.....------------------------------------------------------------------
!c
!c-----> determine exchange correlation potential
!c
      vx (:,:) = 0.0
      vxc(:,:) = 0.0
      ngrid=SIZE(rh,1)
206

207
      IF (xcpot%needs_grad()) THEN
208 209
         IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives")
         IF (xcpot%is_name("l91")) THEN    ! local pw91
210 211 212
            CALL vxcl91(jspins,ngrid,ngrid,rh,grad%agrt(:ngrid),grad%agru(:ngrid),grad%agrd(:ngrid), grad%g2rt(:ngrid),&
                 grad%g2ru(:ngrid),grad%g2rd(:ngrid),grad%gggrt(:ngrid),grad%gggru(:ngrid),grad%gggrd(:ngrid),&
                 grad%gzgr(:ngrid), vx(:ngrid,:),vxc(:ngrid,:), isprsv,sprsv)
213 214
         ELSEIF (xcpot%is_name("pw91")) THEN  ! pw91
            IF (lwbc) THEN
215 216 217
               CALL vxcwb91(jspins,ngrid,ngrid,rh(:ngrid,:),grad%agrt(:ngrid),grad%agru(:ngrid),grad%agrd(:ngrid),&
                 grad%g2rt(:ngrid),grad%g2ru(:ngrid),grad%g2rd(:ngrid),grad%gggrt(:ngrid),grad%gggru(:ngrid),&
                 grad%gggrd(:ngrid),grad%gzgr(:ngrid), vx(:ngrid,:),vxc(:ngrid,:), idsprs,isprsv,sprsv)
218 219
            ELSE

220 221 222
               CALL vxcpw91(jspins,ngrid,ngrid,rh(:ngrid,:),grad%agrt(:ngrid),grad%agru(:ngrid),grad%agrd(:ngrid),&
                 grad%g2rt(:ngrid),grad%g2ru(:ngrid),grad%g2rd(:ngrid),grad%gggrt(:ngrid),grad%gggru(:ngrid),&
                 grad%gggrd,grad%gzgr, vx(:ngrid,:),vxc(:ngrid,:), idsprs,isprsv,sprsv)
223 224 225

            ENDIF
         ELSE  ! pbe or similar
226
            CALL vxcepbe(xcpot%DATA,jspins,ngrid,ngrid,rh(:ngrid,:), grad%agrt,grad%agru,grad%agrd,grad%g2ru,grad%g2rd,grad%gggrt,grad%gggru,grad%gggrd, vx(:ngrid,:),vxc(:ngrid,:))
227 228 229
         ENDIF
      ELSE  !LDA potentials
         IF (xcpot%is_name("x-a"))  THEN   ! X-alpha method
230
            CALL vxcxal(xcpot%data%krla,jspins, ngrid,ngrid,rh(:ngrid,:), vx(:ngrid,:),vxc(:ngrid,:))
231
         ELSEIF (xcpot%is_name("wign")) THEN    ! Wigner interpolation formula
232
            CALL vxcwgn(xcpot%data%krla,jspins, ngrid,ngrid,rh(:ngrid,:), vx(:ngrid,:),vxc(:ngrid,:))
233
         ELSEIF (xcpot%is_name("mjw").OR.xcpot%is_name("bh")) THEN ! von Barth,Hedin correlation
234
            CALL vxcbh(iofile,xcpot%data,jspins, ngrid,ngrid,rh(:ngrid,:), vx(:ngrid,:),vxc(:ngrid,:))
235

236
         ELSEIF (xcpot%is_name("vwn")) THEN     ! Vosko,Wilk,Nusair correlation
237
            CALL vxcvwn(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh(:ngrid,:), vx(:ngrid,:),vxc(:ngrid,:))
238
         ELSEIF (xcpot%is_name("pz")) THEN     ! Perdew,Zunger correlation
239
            CALL vxcpz(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh(:ngrid,:), vx(:ngrid,:),vxc(:ngrid,:))
240 241 242
         ELSEIF (xcpot%is_name("hf")) THEN
            ! Hartree-Fock  calculation: X-alpha potential is added to generate a rational local potential,
            !                            later it is subtracted again
243
            CALL juDFT_error('HF should now be treated as a GGA functional', calledby='xcpot_get_vxc')
244
            CALL vxcxal(xcpot%data%krla,jspins, ngrid,ngrid,rh(:ngrid,:), vx(:ngrid,:),vxc(:ngrid,:))
245 246 247 248 249 250 251 252 253 254 255
            !         vxc=0
         ELSEIF (xcpot%is_name("exx")) THEN
            ! if exact exchange calculation do nothing
            vxc = 0
         ELSE
            CALL juDFT_error("Unkown LDA potential",calledby="type xcpot")
         ENDIF
      ENDIF
!
!-----> hartree units
!
256 257
      vx  = hrtr_half*vx
      vxc = hrtr_half*vxc
258 259

   END SUBROUTINE xcpot_get_vxc
260 261

!***********************************************************************
Matthias Redies's avatar
Matthias Redies committed
262
   SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad,kinEnergyDen_KS)
263 264 265 266 267 268 269 270 271 272 273
!***********************************************************************
      USE m_xcxal, ONLY : excxal
      USE m_xcwgn, ONLY : excwgn
      USE m_xcbh,  ONLY : excbh
      USE m_xcvwn, ONLY : excvwn
      USE m_xcpz,  ONLY : excpz
      USE m_excl91
      USE m_excwb91
      USE m_excpw91
      USE m_excepbe
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
274
      
Matthias Redies's avatar
Matthias Redies committed
275 276 277 278 279
      CLASS(t_xcpot_inbuild),INTENT(IN)     :: xcpot
      INTEGER, INTENT (IN)                  :: jspins
      REAL,INTENT (IN)                      :: rh(:,:)
      REAL, INTENT (OUT)                    :: exc(:)
      TYPE(t_gradients),OPTIONAL,INTENT(IN) ::grad
Matthias Redies's avatar
Matthias Redies committed
280
      REAL, INTENT(IN), OPTIONAL            :: kinEnergyDen_KS(:,:)
281 282 283 284 285 286 287 288 289 290 291 292 293 294 295

!c
!c ---> local scalars
      INTEGER :: ngrid
      REAL, PARAMETER :: hrtr_half = 0.5

      !used to be dummy arguments for testing
      INTEGER,PARAMETER   :: idsprs=0,isprsv=0,iofile=6
      REAL,PARAMETER      :: sprsv=0.0
      LOGICAL,PARAMETER   :: lwbc=.false. ! l-white-bird-current (ta)
!c
!c-----> determine exchange correlation energy density
!c
      exc(:) = 0.0
      ngrid=SIZE(rh,1)
296
      IF (xcpot%exc_is_gga()) THEN
297 298
         IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_exc for a GGA potential without providing derivatives")
         IF (xcpot%is_name("l91")) THEN  ! local pw91
299
            CALL excl91(jspins,ngrid,ngrid,rh(:ngrid,:),grad%agrt,grad%agru,grad%agrd,grad%g2rt,grad%g2ru,grad%g2rd,grad%gggrt,grad%gggru,grad%gggrd,grad%gzgr, exc, isprsv,sprsv)
300 301
         ELSEIF (xcpot%is_name("pw91")) THEN     ! pw91
            IF (lwbc) THEN
302
               CALL excwb91(ngrid,ngrid,rh(:ngrid,1),rh(:ngrid,2),grad%agrt,grad%agru,grad%agrd, grad%g2rt,grad%g2ru,grad%g2rd,grad%gggrt,grad%gggru,grad%gggrd,grad%gzgr, exc, idsprs,isprsv,sprsv)
303
            ELSE
304
               CALL excpw91(jspins,ngrid,ngrid,rh(:ngrid,:),grad%agrt,grad%agru,grad%agrd, grad%g2rt,grad%g2ru,grad%g2rd,grad%gggrt,grad%gggru,grad%gggrd,grad%gzgr, exc, idsprs,isprsv,sprsv)
305 306
            ENDIF
         ELSE
307
            CALL excepbe(xcpot%data,jspins,ngrid,ngrid, rh(:ngrid,:),grad%agrt,grad%agru,grad%agrd,grad%g2ru,grad%g2rd,grad%gggrt,grad%gggru,grad%gggrd, exc)
308 309 310 311 312 313 314 315 316 317 318 319 320
         ENDIF
      ELSE !LDA
         IF (xcpot%is_name("x-a"))  THEN   ! X-alpha method
            CALL excxal(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh, exc)
         ELSEIF (xcpot%is_name("wign")) THEN    ! Wigner interpolation formula
            CALL excwgn(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh, exc)
         ELSEIF (xcpot%is_name("mjw").OR.xcpot%is_name("bh")) THEN ! von Barth,Hedin correlation
            CALL excbh(iofile,xcpot%data,jspins, ngrid,ngrid,rh, exc)
         ELSEIF (xcpot%is_name("vwn")) THEN     ! Vosko,Wilk,Nusair correlation
            CALL excvwn(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh, exc)
         ELSEIF (xcpot%is_name("pz")) THEN     ! Perdew,Zunger correlation
            CALL excpz(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh, exc)
         ELSEIF (xcpot%is_name("hf") .OR. xcpot%is_name("exx")) THEN
321
            CALL juDFT_error('HF should now be treated as a GGA functional', calledby='xcpot_get_exc')
322 323 324 325 326 327 328 329
            exc=0
         ELSE
            CALL juDFT_error("Unkown LDA potential",calledby="type xcpot")
         ENDIF
      ENDIF
!c-----> hartree units
      exc= hrtr_half*exc

330 331 332
   END SUBROUTINE xcpot_get_exc

   LOGICAL FUNCTION xcpot_is_name(xcpot,name)
333 334 335
      CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot
      CHARACTER(len=*),INTENT(IN)  :: name
      xcpot_is_name=(TRIM(xc_names(xcpot%icorr))==TRIM((name)))
336
   END FUNCTION xcpot_is_name
337

338
END MODULE m_types_xcpot_inbuild