types_xcpot_libxc.F90 14.2 KB
Newer Older
1 2 3 4 5 6 7 8
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------

!>This module contains the xcpot-type providing an interface to libxc
MODULE m_types_xcpot_libxc
9 10
#ifdef CPP_LIBXC
   USE xc_f03_lib_m
11
#endif
Matthias Redies's avatar
Matthias Redies committed
12 13 14
   USE m_types_xcpot
   USE m_judft
   IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
15 16 17


#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
18
   PRIVATE :: write_xc_info
Matthias Redies's avatar
Matthias Redies committed
19
#endif
20

Matthias Redies's avatar
Matthias Redies committed
21
   TYPE,EXTENDS(t_xcpot):: t_xcpot_libxc 
22
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
23 24
      TYPE(xc_f03_func_t)      :: vxc_func_x, vxc_func_c
      TYPE(xc_f03_func_t)      :: exc_func_x, exc_func_c
25
#endif
Matthias Redies's avatar
Matthias Redies committed
26 27 28
      INTEGER          :: func_vxc_id_c, func_vxc_id_x !> functionals to be used for potential & density convergence
      INTEGER          :: func_exc_id_c, func_exc_id_x !> functionals to be used in exc- & totale-calculations
      INTEGER          :: jspins
29
   CONTAINS
Matthias Redies's avatar
Matthias Redies committed
30 31
      PROCEDURE        :: vxc_is_LDA          => xcpot_vxc_is_LDA
      PROCEDURE        :: exc_is_LDA          => xcpot_exc_is_LDA
32 33
      PROCEDURE        :: vxc_is_gga          => xcpot_vxc_is_gga
      PROCEDURE        :: exc_is_gga          => xcpot_exc_is_gga
Matthias Redies's avatar
Matthias Redies committed
34
      PROCEDURE        :: exc_is_MetaGGA      => xcpot_exc_is_MetaGGA
Matthias Redies's avatar
Matthias Redies committed
35 36 37 38 39 40 41 42 43
      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
      PROCEDURE,NOPASS :: alloc_gradients     => xcpot_alloc_gradients
      !Not             overloeaded...
      PROCEDURE        :: init                => xcpot_init
   END TYPE t_xcpot_libxc
   PUBLIC t_xcpot_libxc
44 45
CONTAINS

Matthias Redies's avatar
Matthias Redies committed
46 47 48 49 50 51
   SUBROUTINE xcpot_init(xcpot,jspins,vxc_id_x,vxc_id_c,exc_id_x,exc_id_c)
      USE m_judft
      IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(OUT)    :: xcpot
      INTEGER, INTENT(IN)                 :: jspins
      INTEGER, INTENT(IN)                 :: vxc_id_x, vxc_id_c ! potential functional
52
      INTEGER, INTENT(IN)                 :: exc_id_x, exc_id_c ! energy functionals
Matthias Redies's avatar
Matthias Redies committed
53 54
      LOGICAL                             :: same_functionals   ! are vxc and exc equal

55
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
56 57 58 59 60
      INTEGER :: err
      xcpot%jspins = jspins

      xcpot%func_vxc_id_x = vxc_id_x
      xcpot%func_vxc_id_c = vxc_id_c
Matthias Redies's avatar
Matthias Redies committed
61
      
62 63
      xcpot%func_exc_id_x = exc_id_x
      xcpot%func_exc_id_c = exc_id_c
Matthias Redies's avatar
Matthias Redies committed
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80

      if(xcpot%func_vxc_id_x == 0 .or. xcpot%func_exc_id_x == 0 ) then
         CALL judft_error("LibXC exchange- and correlation-function indicies need to be set"&
            ,hint='Try this: ' // ACHAR(10) //&
            '<xcFunctional name="libxc" relativisticCorrections="F">' // ACHAR(10) //& 
            '  <libXC  exchange="1" correlation="1" /> ' // ACHAR(10) //& 
            '</xcFunctional> ')
      endif 

      IF (jspins==1) THEN
         ! potential functionals      
         CALL xc_f03_func_init(xcpot%vxc_func_x, xcpot%func_vxc_id_x, XC_UNPOLARIZED)
         IF (xcpot%func_vxc_id_c>0) CALL xc_f03_func_init(xcpot%vxc_func_c, xcpot%func_vxc_id_c, XC_UNPOLARIZED)

         ! energy functionals
         CALL xc_f03_func_init(xcpot%exc_func_x, xcpot%func_exc_id_x, XC_UNPOLARIZED)
         IF (xcpot%func_exc_id_c>0) CALL xc_f03_func_init(xcpot%exc_func_c, xcpot%func_exc_id_c, XC_UNPOLARIZED)
Matthias Redies's avatar
Matthias Redies committed
81
      ELSE
Matthias Redies's avatar
Matthias Redies committed
82 83 84 85 86 87 88 89 90
         ! potential functionals  
         CALL xc_f03_func_init(xcpot%vxc_func_x, xcpot%func_vxc_id_x, XC_POLARIZED)
         IF (xcpot%func_vxc_id_c>0) CALL xc_f03_func_init(xcpot%vxc_func_c, xcpot%func_vxc_id_c, XC_POLARIZED)

         !energy functionals
         CALL xc_f03_func_init(xcpot%exc_func_x, xcpot%func_exc_id_x, XC_POLARIZED)
         IF (xcpot%func_exc_id_c>0) CALL xc_f03_func_init(xcpot%exc_func_c, xcpot%func_exc_id_c, XC_POLARIZED)
      END IF

91
      !check if any potental is a MetaGGA
Matthias Redies's avatar
Matthias Redies committed
92
      IF(ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA] == xc_get_family(xcpot%vxc_func_x))) THEN
93
         call juDFT_error("vxc_x: MetaGGA is not implemented for potentials")
Matthias Redies's avatar
Matthias Redies committed
94 95
      ELSEIF(xcpot%func_vxc_id_c>0) THEN
         IF(ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA] == xc_get_family(xcpot%vxc_func_c))) THEN
96 97 98 99
            call juDFT_error("vxc_x: MetaGGA is not implemented for potentials")
         ENDIF
      ENDIF

Matthias Redies's avatar
Matthias Redies committed
100 101 102 103 104 105 106 107
      CALL write_xc_info(xcpot%vxc_func_x)

      IF (xcpot%func_vxc_id_c>0) THEN
         CALL write_xc_info(xcpot%vxc_func_c)
      ELSE
         WRITE(*,*) "No Correlation functional"
      END IF

108 109
      same_functionals =     (xcpot%func_vxc_id_x == xcpot%func_exc_id_x) &
                       .and. (xcpot%func_vxc_id_c == xcpot%func_exc_id_c)
Matthias Redies's avatar
Matthias Redies committed
110 111 112 113 114 115 116
      IF(.not. same_functionals) THEN
         CALL write_xc_info(xcpot%exc_func_x)
         IF (xcpot%func_exc_id_c>0) THEN
            CALL write_xc_info(xcpot%exc_func_c)
         ELSE
            WRITE(*,*) "No Correlation functional for TotalE"
         ENDIF
Matthias Redies's avatar
Matthias Redies committed
117 118
      ELSE
         write (*,*) "Using same functional for VXC and EXC"
Matthias Redies's avatar
Matthias Redies committed
119
      END IF
120
#else
121
      CALL judft_error("You specified a libxc-exchange correlation potential but FLEUR is not linked against libxc", &
Matthias Redies's avatar
Matthias Redies committed
122
         hint="Please recompile FLEUR with libxc support")
123
#endif
124
   END SUBROUTINE xcpot_init
125

Matthias Redies's avatar
Matthias Redies committed
126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153
   ! LDA
   LOGICAL FUNCTION xcpot_vxc_is_LDA(xcpot)
      IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC
      TYPE(xc_f03_func_info_t)        :: xc_info

      xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
      xcpot_vxc_is_LDA =  XC_FAMILY_LDA == xc_f03_func_info_get_family(xc_info)
#else
      xcpot_vxc_is_LDA = .false.
#endif
   END FUNCTION xcpot_vxc_is_LDA

   LOGICAL FUNCTION xcpot_exc_is_LDA(xcpot)
      IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC
      TYPE(xc_f03_func_info_t)        :: xc_info

      xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
      xcpot_exc_is_LDA = XC_FAMILY_LDA == xc_f03_func_info_get_family(xc_info)
#else
      xcpot_exc_is_LDA = .false.
#endif
   END FUNCTION xcpot_exc_is_LDA

   ! GGA
154
   LOGICAL FUNCTION xcpot_vxc_is_gga(xcpot)
155
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
156
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
157
#ifdef CPP_LIBXC
158 159
      TYPE(xc_f03_func_info_t)        :: xc_info

160
      xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
Matthias Redies's avatar
Matthias Redies committed
161
      xcpot_vxc_is_gga =  ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f03_func_info_get_family(xc_info))
162 163 164 165 166 167 168 169 170 171 172 173
#else
      xcpot_vxc_is_gga=.false.
#endif
   END FUNCTION xcpot_vxc_is_gga

   LOGICAL FUNCTION xcpot_exc_is_gga(xcpot)
      IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC
      TYPE(xc_f03_func_info_t)        :: xc_info

      xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
Matthias Redies's avatar
Matthias Redies committed
174
      xcpot_exc_is_gga =  ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f03_func_info_get_family(xc_info))
175
#else
176
      xcpot_exc_is_gga=.false.
177
#endif
178
   END FUNCTION xcpot_exc_is_gga
179

Matthias Redies's avatar
Matthias Redies committed
180
   LOGICAL FUNCTION xcpot_exc_is_MetaGGA(xcpot)
Matthias Redies's avatar
Matthias Redies committed
181 182
      IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
Matthias Redies's avatar
Matthias Redies committed
183
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
184
      TYPE(xc_f03_func_info_t)        :: xc_info
Matthias Redies's avatar
Matthias Redies committed
185

Matthias Redies's avatar
Matthias Redies committed
186
      xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
Matthias Redies's avatar
Matthias Redies committed
187
      xcpot_exc_is_MetaGGA=ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA]==xc_f03_func_info_get_family(xc_info))
188
#else
Matthias Redies's avatar
Matthias Redies committed
189
      xcpot_exc_is_MetaGGA =  .False.
Matthias Redies's avatar
Matthias Redies committed
190
#endif
Matthias Redies's avatar
Matthias Redies committed
191
   END FUNCTION xcpot_exc_is_MetaGGA
Matthias Redies's avatar
Matthias Redies committed
192

193 194
   LOGICAL FUNCTION xcpot_is_hybrid(xcpot)
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
195
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
196
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
197
      TYPE(xc_f03_func_info_t)        :: xc_info
Matthias Redies's avatar
Matthias Redies committed
198

Matthias Redies's avatar
Matthias Redies committed
199 200
      xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
      xcpot_is_hybrid=ANY([XC_FAMILY_HYB_MGGA, XC_FAMILY_HYB_GGA]==xc_f03_func_info_get_family(xc_info))
201
#else
Matthias Redies's avatar
Matthias Redies committed
202
      xcpot_is_hybrid = .False.
203
#endif
204
   END FUNCTION xcpot_is_hybrid
205

206 207 208
   FUNCTION xcpot_get_exchange_weight(xcpot) RESULT(a_ex)
      USE m_judft
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
209
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
210

Matthias Redies's avatar
Matthias Redies committed
211
      REAL:: a_ex
212
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
213
      a_ex=xc_f03_hyb_exx_coef(xcpot%vxc_func_x)
214
#endif
Matthias Redies's avatar
Matthias Redies committed
215 216 217 218 219 220 221 222 223 224 225 226 227
   END FUNCTION xcpot_get_exchange_weight

   !***********************************************************************
   SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad)
      !***********************************************************************
      IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
      INTEGER, INTENT (IN)     :: jspins
      REAL,INTENT (IN)         :: rh(:,:)   !Dimensions here
      REAL, INTENT (OUT)       :: vx (:,:)  !points,spin
      REAL, INTENT (OUT  )     :: vxc(:,:)  !
      ! optional arguments for GGA
      TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
228
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
229 230 231 232
      REAL,ALLOCATABLE::vxc_tmp(:,:),vx_tmp(:,:),vsigma(:,:)
      !libxc uses the spin as a first index, hence we have to transpose....
      ALLOCATE(vxc_tmp(SIZE(vxc,2),SIZE(vxc,1)));vxc_tmp=0.0
      ALLOCATE(vx_tmp(SIZE(vx,2),SIZE(vx,1)));vx_tmp=0.0
233
      IF (xcpot%needs_grad()) THEN
Matthias Redies's avatar
Matthias Redies committed
234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
         IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives")
         ALLOCATE(vsigma,mold=grad%vsigma)
         !where(abs(grad%sigma)<1E-9) grad%sigma=1E-9
         CALL xc_f03_gga_vxc(xcpot%vxc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vx_tmp,vsigma)
         IF (xcpot%func_vxc_id_c>0) THEN
            CALL xc_f03_gga_vxc(xcpot%vxc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vxc_tmp,grad%vsigma)
            grad%vsigma=grad%vsigma+vsigma
            vxc_tmp=vxc_tmp+vx_tmp
         ELSE     
            vxc_tmp=vx_tmp
         ENDIF
      ELSE  !LDA potentials
         CALL xc_f03_lda_vxc(xcpot%vxc_func_x, SIZE(rh,1), TRANSPOSE(rh), vx_tmp)
         IF (xcpot%func_vxc_id_c>0) THEN
            CALL xc_f03_lda_vxc(xcpot%vxc_func_c, SIZE(rh,1), TRANSPOSE(rh), vxc_tmp)
            vxc_tmp=vxc_tmp+vx_tmp
         ENDIF
      ENDIF
      vx=TRANSPOSE(vx_tmp) 
      vxc=TRANSPOSE(vxc_tmp) 

255
#endif
Matthias Redies's avatar
Matthias Redies committed
256 257 258 259 260 261 262 263 264
   END SUBROUTINE xcpot_get_vxc


   SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad)
      IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
      INTEGER, INTENT (IN)     :: jspins
      REAL,INTENT (IN) :: rh(:,:)  !points,spin
      REAL, INTENT (OUT) :: exc(:) !points
265
      ! optional arguments for GGA
Matthias Redies's avatar
Matthias Redies committed
266 267 268
      TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad

      REAL  :: excc(SIZE(exc))
269
#ifdef CPP_LIBXC
270
      IF (xcpot%exc_is_gga()) THEN
Matthias Redies's avatar
Matthias Redies committed
271 272
         IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives")
         CALL xc_f03_gga_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,exc)
273
         IF (xcpot%func_exc_id_c>0) THEN
Matthias Redies's avatar
Matthias Redies committed
274 275 276
            CALL xc_f03_gga_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,excc)
            exc=exc+excc
         END IF
277
      ELSEIF(xcpot%exc_is_LDA()) THEN  !LDA potentials
Matthias Redies's avatar
Matthias Redies committed
278
         CALL xc_f03_lda_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh), exc)
279
         IF (xcpot%func_exc_id_c>0) THEN
Matthias Redies's avatar
Matthias Redies committed
280 281 282
            CALL xc_f03_lda_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh), excc)
            exc=exc+excc
         END IF
283 284 285 286 287
      ELSEIF(xcpot%exc_is_MetaGGA()) THEN
         !call xc_f03_mgga_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh), grad%sigma, transpose(grad%laplace)
         write (*,*) "have to implement that"
      ELSE
         call juDFT_error("exc is part of a known Family", calledby="xcpot_get_exc@libxc")
Matthias Redies's avatar
Matthias Redies committed
288 289
      ENDIF

290 291
#endif
   END SUBROUTINE xcpot_get_exc
292

293 294 295 296 297 298 299 300 301
   SUBROUTINE xcpot_alloc_gradients(ngrid,jspins,grad)
      INTEGER, INTENT (IN)         :: jspins,ngrid
      TYPE(t_gradients),INTENT(INOUT):: grad
      !For libxc we only need the sigma array...
      IF (ALLOCATED(grad%sigma)) DEALLOCATE(grad%sigma,grad%gr,grad%laplace,grad%vsigma)
      ALLOCATE(grad%sigma(MERGE(1,3,jspins==1),ngrid))
      ALLOCATE(grad%gr(3,ngrid,jspins))
      ALLOCATE(grad%laplace(ngrid,jspins))
      ALLOCATE(grad%vsigma(MERGE(1,3,jspins==1),ngrid))
302

303
   END SUBROUTINE xcpot_alloc_gradients
304

Matthias Redies's avatar
Matthias Redies committed
305
#ifdef CPP_LIBXC  
Matthias Redies's avatar
Matthias Redies committed
306 307 308 309 310 311 312
   SUBROUTINE write_xc_info(xc_func, is_E_func)
      IMPLICIT NONE
      LOGICAL,INTENT(IN),OPTIONAL         :: is_E_func
      INTEGER                             :: i
      CHARACTER(len=120)                  :: kind, family
      LOGICAL                             :: is_energy_func

313 314 315
      TYPE(xc_f03_func_t),INTENT(IN)      :: xc_func
      TYPE(xc_f03_func_info_t)            :: xc_info

Matthias Redies's avatar
Matthias Redies committed
316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357
      xc_info = xc_f03_func_get_info(xc_func)
      is_energy_func = merge(is_E_func, .False., PRESENT(is_E_func))

      SELECT CASE(xc_f03_func_info_get_kind(xc_info))
      CASE (XC_EXCHANGE)
         WRITE(kind, '(a)') 'an exchange functional'
      CASE (XC_CORRELATION)
         WRITE(kind, '(a)') 'a correlation functional'
      CASE (XC_EXCHANGE_CORRELATION)
         WRITE(kind, '(a)') 'an exchange-correlation functional'
      CASE (XC_KINETIC)
         WRITE(kind, '(a)') 'a kinetic energy functional'
      CASE default
         WRITE(kind, '(a)') 'of unknown kind'
      END SELECT
      SELECT CASE (xc_f03_func_info_get_family(xc_info))
      CASE (XC_FAMILY_LDA);
         WRITE(family,'(a)') "LDA"
      CASE (XC_FAMILY_GGA);
         WRITE(family,'(a)') "GGA"
      CASE (XC_FAMILY_HYB_GGA);
         WRITE(family,'(a)') "Hybrid GGA"
      CASE (XC_FAMILY_MGGA);
         WRITE(family,'(a)') "MGGA"
      CASE (XC_FAMILY_HYB_MGGA);
         WRITE(family,'(a)') "Hybrid MGGA"
      CASE default;
         WRITE(family,'(a)') "unknown"
      END SELECT

      IF(.not. is_energy_func) THEN
         WRITE(*,'("The functional ''", a, "'' is ", a, ", it belongs to the ''", a, "'' family and is defined in the reference(s):")') &
            TRIM(xc_f03_func_info_get_name(xc_info)), TRIM(kind), TRIM(family)
      ELSE
         WRITE(*,'("The functional used for TotalE ''", a, "'' is ", a, ", it belongs to the ''", a, "'' family and is defined in the reference(s):")') &
            TRIM(xc_f03_func_info_get_name(xc_info)), TRIM(kind), TRIM(family)
      ENDIF

      i = 0
      DO WHILE(i >= 0)
         WRITE(*, '(a,i1,2a)') '[', i+1, '] ', TRIM(xc_f03_func_reference_get_ref(xc_f03_func_info_get_references(xc_info, i)))
      END DO
358
   END SUBROUTINE write_xc_info
359 360 361 362 363 364 365 366

   FUNCTION xc_get_family(xc_func) result(family)
      IMPLICIT NONE
      TYPE(xc_f03_func_t)  :: xc_func  
      integer              :: family
      family = xc_f03_func_info_get_family(xc_f03_func_get_info(xc_func))
   END FUNCTION xc_get_family

Matthias Redies's avatar
Matthias Redies committed
367
#endif
368

369
END MODULE m_types_xcpot_libxc