types_xcpot_inbuild.F90 13.9 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  ']
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_vxc_is_LDA
      PROCEDURE        :: exc_is_LDA=>xcpot_exc_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 119 120 121 122 123 124 125 126 127 128 129 130 131 132
  
   !! LDA

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

   LOGICAL FUNCTION xcpot_exc_is_LDA(xcpot)
      IMPLICIT NONE
      CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot
      xcpot_exc_is_LDA = xcpot%vxc_is_LDA()
   END FUNCTION xcpot_exc_is_LDA

   !! GGA
133 134 135 136 137 138 139 140

   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)
141 142
      IMPLICIT NONE
      CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot
143 144
      xcpot_exc_is_gga = xcpot%vxc_is_gga()
   END FUNCTION xcpot_exc_is_gga
145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166

   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)
167
!
168 169 170 171 172 173 174 175 176
      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
177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
      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
194
      TYPE(t_gradients),INTENT(INOUT),OPTIONAL::grad
195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211
!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)
212

213
      IF (xcpot%needs_grad()) THEN
214 215
         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
216 217 218
            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)
219 220
         ELSEIF (xcpot%is_name("pw91")) THEN  ! pw91
            IF (lwbc) THEN
221 222 223
               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)
224 225
            ELSE

226 227 228
               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)
229 230 231

            ENDIF
         ELSE  ! pbe or similar
232
            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,:))
233 234 235
         ENDIF
      ELSE  !LDA potentials
         IF (xcpot%is_name("x-a"))  THEN   ! X-alpha method
236
            CALL vxcxal(xcpot%data%krla,jspins, ngrid,ngrid,rh(:ngrid,:), vx(:ngrid,:),vxc(:ngrid,:))
237
         ELSEIF (xcpot%is_name("wign")) THEN    ! Wigner interpolation formula
238
            CALL vxcwgn(xcpot%data%krla,jspins, ngrid,ngrid,rh(:ngrid,:), vx(:ngrid,:),vxc(:ngrid,:))
239
         ELSEIF (xcpot%is_name("mjw").OR.xcpot%is_name("bh")) THEN ! von Barth,Hedin correlation
240
            CALL vxcbh(iofile,xcpot%data,jspins, ngrid,ngrid,rh(:ngrid,:), vx(:ngrid,:),vxc(:ngrid,:))
241

242
         ELSEIF (xcpot%is_name("vwn")) THEN     ! Vosko,Wilk,Nusair correlation
243
            CALL vxcvwn(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh(:ngrid,:), vx(:ngrid,:),vxc(:ngrid,:))
244
         ELSEIF (xcpot%is_name("pz")) THEN     ! Perdew,Zunger correlation
245
            CALL vxcpz(iofile,xcpot%data%krla,jspins, ngrid,ngrid,rh(:ngrid,:), vx(:ngrid,:),vxc(:ngrid,:))
246 247 248
         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
249
            CALL juDFT_error('HF should now be treated as a GGA functional', calledby='xcpot_get_vxc')
250
            CALL vxcxal(xcpot%data%krla,jspins, ngrid,ngrid,rh(:ngrid,:), vx(:ngrid,:),vxc(:ngrid,:))
251 252 253 254 255 256 257 258 259 260 261
            !         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
!
262 263
      vx  = hrtr_half*vx
      vxc = hrtr_half*vxc
264 265

   END SUBROUTINE xcpot_get_vxc
266 267

!***********************************************************************
268
   SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad)
269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
!***********************************************************************
      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
!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(:,:)
289
      
290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310
!c
!c---> xc energy density
!c
      REAL, INTENT (OUT) :: exc(:)

      ! optional arguments for GGA
      TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad
!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)
311
      IF (xcpot%exc_is_gga()) THEN
312 313
         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
314
            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)
315 316
         ELSEIF (xcpot%is_name("pw91")) THEN     ! pw91
            IF (lwbc) THEN
317
               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)
318
            ELSE
319
               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)
320 321
            ENDIF
         ELSE
322
            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)
323 324 325 326 327 328 329 330 331 332 333 334 335
         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
336
            CALL juDFT_error('HF should now be treated as a GGA functional', calledby='xcpot_get_exc')
337 338 339 340 341 342 343 344
            exc=0
         ELSE
            CALL juDFT_error("Unkown LDA potential",calledby="type xcpot")
         ENDIF
      ENDIF
!c-----> hartree units
      exc= hrtr_half*exc

345 346 347
   END SUBROUTINE xcpot_get_exc

   LOGICAL FUNCTION xcpot_is_name(xcpot,name)
348 349 350
      CLASS(t_xcpot_inbuild),INTENT(IN):: xcpot
      CHARACTER(len=*),INTENT(IN)  :: name
      xcpot_is_name=(TRIM(xc_names(xcpot%icorr))==TRIM((name)))
351
   END FUNCTION xcpot_is_name
352

353
END MODULE m_types_xcpot_inbuild