types_xcpot_libxc.F90 21 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
   USE m_types_xcpot
   USE m_judft
Matthias Redies's avatar
Matthias Redies committed
14
   use m_types_misc
Matthias Redies's avatar
Matthias Redies committed
15
   IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
16 17 18


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

Matthias Redies's avatar
Matthias Redies committed
22
   TYPE,EXTENDS(t_xcpot):: t_xcpot_libxc 
23
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
24 25
      TYPE(xc_f03_func_t)      :: vxc_func_x, vxc_func_c
      TYPE(xc_f03_func_t)      :: exc_func_x, exc_func_c
26
#endif
Matthias Redies's avatar
Matthias Redies committed
27 28 29
      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
30
   CONTAINS
Matthias Redies's avatar
Matthias Redies committed
31 32 33
      PROCEDURE        :: vx_is_LDA     => xcpot_vx_is_LDA
      PROCEDURE        :: vx_is_GGA     => xcpot_vx_is_GGA
      PROCEDURE        :: vx_is_MetaGGA => xcpot_vx_is_MetaGGA
34 35 36 37 38 39 40 41
      
      PROCEDURE        :: vc_is_LDA => xcpot_vc_is_LDA
      PROCEDURE        :: vc_is_GGA => xcpot_vc_is_GGA

      PROCEDURE        :: exc_is_LDA     => xcpot_exc_is_LDA
      PROCEDURE        :: exc_is_gga     => xcpot_exc_is_gga
      PROCEDURE        :: exc_is_MetaGGA => xcpot_exc_is_MetaGGA

Matthias Redies's avatar
Matthias Redies committed
42 43 44 45 46 47 48 49 50
      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
51 52
CONTAINS

Matthias Redies's avatar
Matthias Redies committed
53 54 55 56 57 58
   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
59
      INTEGER, INTENT(IN)                 :: exc_id_x, exc_id_c ! energy functionals
Matthias Redies's avatar
Matthias Redies committed
60 61
      LOGICAL                             :: same_functionals   ! are vxc and exc equal

62
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
63 64 65 66 67
      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
68
      
69 70
      xcpot%func_exc_id_x = exc_id_x
      xcpot%func_exc_id_c = exc_id_c
Matthias Redies's avatar
Matthias Redies committed
71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87

      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
88
         
Matthias Redies's avatar
Matthias Redies committed
89
      ELSE
Matthias Redies's avatar
Matthias Redies committed
90 91 92 93 94 95 96 97 98
         ! 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

99
      !check if any potental is a MetaGGA
Matthias Redies's avatar
Matthias Redies committed
100 101
      IF(ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA] == xc_get_family(xcpot%vxc_func_c))) THEN
         call juDFT_error("vxc_x: MetaGGA is not implemented for correclation potentials")
102 103
      ENDIF

Matthias Redies's avatar
Matthias Redies committed
104 105 106 107 108 109 110 111
      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

112 113
      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
114 115 116 117 118 119 120
      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
121 122
      ELSE
         write (*,*) "Using same functional for VXC and EXC"
Matthias Redies's avatar
Matthias Redies committed
123
      END IF
Matthias Redies's avatar
Matthias Redies committed
124 125

      xcpot%kinED%set = .False.
126
#else
127
      CALL judft_error("You specified a libxc-exchange correlation potential but FLEUR is not linked against libxc", &
Matthias Redies's avatar
Matthias Redies committed
128
         hint="Please recompile FLEUR with libxc support")
129
#endif
130
   END SUBROUTINE xcpot_init
131

Matthias Redies's avatar
Matthias Redies committed
132
   ! LDA
133
   LOGICAL FUNCTION xcpot_vx_is_LDA(xcpot)
Matthias Redies's avatar
Matthias Redies committed
134 135 136 137 138 139
      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)
140
      xcpot_vx_is_LDA =  XC_FAMILY_LDA == xc_f03_func_info_get_family(xc_info)
Matthias Redies's avatar
Matthias Redies committed
141
#else
142
      xcpot_vx_is_LDA = .false.
Matthias Redies's avatar
Matthias Redies committed
143
#endif
144 145 146 147 148 149 150
   END FUNCTION xcpot_vx_is_LDA
   
   LOGICAL FUNCTION xcpot_vc_is_LDA(xcpot)
      IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC
      TYPE(xc_f03_func_info_t)        :: xc_info
Matthias Redies's avatar
Matthias Redies committed
151

152 153 154 155 156 157 158
      xc_info = xc_f03_func_get_info(xcpot%vxc_func_c)
      xcpot_vc_is_LDA =  XC_FAMILY_LDA == xc_f03_func_info_get_family(xc_info)
#else
      xcpot_vc_is_LDA = .false.
#endif
   END FUNCTION xcpot_vc_is_LDA
   
Matthias Redies's avatar
Matthias Redies committed
159 160 161 162 163 164
   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

Matthias Redies's avatar
Matthias Redies committed
165 166
      xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
      xcpot_exc_is_LDA = (XC_FAMILY_LDA == xc_f03_func_info_get_family(xc_info))
Matthias Redies's avatar
Matthias Redies committed
167 168 169 170 171 172
#else
      xcpot_exc_is_LDA = .false.
#endif
   END FUNCTION xcpot_exc_is_LDA

   ! GGA
173
   LOGICAL FUNCTION xcpot_vc_is_gga(xcpot)
174
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
175
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
176
#ifdef CPP_LIBXC
177 178
      TYPE(xc_f03_func_info_t)        :: xc_info

179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
      xc_info = xc_f03_func_get_info(xcpot%vxc_func_c)
      xcpot_vc_is_gga =  ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f03_func_info_get_family(xc_info))
#else
      xcpot_vc_is_gga=.false.
#endif
   END FUNCTION xcpot_vc_is_gga

   LOGICAL FUNCTION xcpot_vx_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_c)
      xcpot_vx_is_gga =  ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f03_func_info_get_family(xc_info))
194
#else
195
      xcpot_vx_is_gga=.false.
196
#endif
197
   END FUNCTION xcpot_vx_is_gga
198

Matthias Redies's avatar
Matthias Redies committed
199 200 201 202 203 204
   LOGICAL FUNCTION xcpot_vx_is_MetaGGA(xcpot)
      IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC
      TYPE(xc_f03_func_info_t)        :: xc_info

Matthias Redies's avatar
Matthias Redies committed
205
      xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
Matthias Redies's avatar
Matthias Redies committed
206 207 208 209 210 211
      xcpot_vx_is_MetaGGA =  ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA]==xc_f03_func_info_get_family(xc_info))
#else
      xcpot_vx_is_MetaGGA=.false.
#endif
   END FUNCTION xcpot_vx_is_MetaGGA

212 213 214 215 216 217
   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

Matthias Redies's avatar
Matthias Redies committed
218
      xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
Matthias Redies's avatar
Matthias Redies committed
219
      xcpot_exc_is_gga =  ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f03_func_info_get_family(xc_info))
220
#else
221
      xcpot_exc_is_gga=.false.
222
#endif
223
   END FUNCTION xcpot_exc_is_gga
224

Matthias Redies's avatar
Matthias Redies committed
225
   LOGICAL FUNCTION xcpot_exc_is_MetaGGA(xcpot)
Matthias Redies's avatar
Matthias Redies committed
226 227
      IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
Matthias Redies's avatar
Matthias Redies committed
228
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
229
      TYPE(xc_f03_func_info_t)        :: xc_info
Matthias Redies's avatar
Matthias Redies committed
230

Matthias Redies's avatar
Matthias Redies committed
231
      xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
Matthias Redies's avatar
Matthias Redies committed
232
      xcpot_exc_is_MetaGGA=ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA]==xc_f03_func_info_get_family(xc_info))
233
#else
Matthias Redies's avatar
Matthias Redies committed
234
      xcpot_exc_is_MetaGGA =  .False.
Matthias Redies's avatar
Matthias Redies committed
235
#endif
Matthias Redies's avatar
Matthias Redies committed
236
   END FUNCTION xcpot_exc_is_MetaGGA
Matthias Redies's avatar
Matthias Redies committed
237

238 239
   LOGICAL FUNCTION xcpot_is_hybrid(xcpot)
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
240
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
241
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
242
      TYPE(xc_f03_func_info_t)        :: xc_info
Matthias Redies's avatar
Matthias Redies committed
243

Matthias Redies's avatar
Matthias Redies committed
244 245
      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))
246
#else
Matthias Redies's avatar
Matthias Redies committed
247
      xcpot_is_hybrid = .False.
248
#endif
249
   END FUNCTION xcpot_is_hybrid
250

251 252 253
   FUNCTION xcpot_get_exchange_weight(xcpot) RESULT(a_ex)
      USE m_judft
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
254
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
255

Matthias Redies's avatar
Matthias Redies committed
256
      REAL:: a_ex
257
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
258
      a_ex=xc_f03_hyb_exx_coef(xcpot%vxc_func_x)
259
#endif
Matthias Redies's avatar
Matthias Redies committed
260 261 262
   END FUNCTION xcpot_get_exchange_weight

   !***********************************************************************
Matthias Redies's avatar
Matthias Redies committed
263
   SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad, kinED_KS)
264
      USE, INTRINSIC :: IEEE_ARITHMETIC
Matthias Redies's avatar
Matthias Redies committed
265 266 267 268 269 270
      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(:,:)  !
Matthias Redies's avatar
Matthias Redies committed
271
      REAL, INTENT(IN),OPTIONAL:: kinED_KS(:,:)
Matthias Redies's avatar
Matthias Redies committed
272 273
      ! optional arguments for GGA
      TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
274
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
275 276 277
      REAL,ALLOCATABLE  :: vxc_tmp(:,:),vx_tmp(:,:),vsigma(:,:), &
                           tmp_vsig(:,:), tmp_vlapl(:,:), tmp_vtau(:,:), &
                           kinED_libxc(:,:)
278
      integer           :: idx
Matthias Redies's avatar
Matthias Redies committed
279 280 281
      !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
Matthias Redies's avatar
Matthias Redies committed
282

283 284 285 286
      write (*,*) "vgen call"
      if(present(kinED_KS)) then
         write (*,*) "max(abs(kED)", maxval(abs(kinED_KS))
      endif
Matthias Redies's avatar
Matthias Redies committed
287
      IF (xcpot%vc_is_GGA()) THEN
288
         write (*,*) "somehow in GGA branch"
Matthias Redies's avatar
Matthias Redies committed
289 290 291 292 293 294 295 296 297 298 299 300
         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
301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320
         write (*,*) "Is MGGA          = ", xcpot%vx_is_MetaGGA()
         write (*,*) "kinED_KS present = ", present(kinED_KS)
         if(xcpot%vx_is_MetaGGA() .and. present(kinED_KS)) then
            kinED_libXC = transpose(kinED_KS + 0.25 * grad%laplace)

            allocate(tmp_vsig,  mold=vx_tmp)
            allocate(tmp_vlapl, mold=vx_tmp)
            allocate(tmp_vtau,  mold=vx_tmp)
            
            write (*,*) "using = ", trim(xc_f03_func_info_get_name(&
                                      xc_f03_func_get_info(xcpot%vxc_func_x)))

            call xc_f03_mgga_vxc(xcpot%vxc_func_x, size(rh,1), transpose(rh), &
                                 grad%sigma, transpose(grad%laplace), kinED_libxc,&
                                 vx_tmp, tmp_vsig, tmp_vlapl, tmp_vtau)

            idx = find_first_normal(vx_tmp)+1
            vx_tmp(:,:idx) = 0.0
            CALL xc_f03_lda_vxc(initial_lda_func(jspins), idx, TRANSPOSE(rh(:idx,:)), vx_tmp(:,:idx))
            write (*,*) "vx has nan = ", any(ieee_is_nan(vx_tmp))
Matthias Redies's avatar
Matthias Redies committed
321
         else
322 323 324
            write (*,*) "using = ", trim(xc_f03_func_info_get_name(&
                                      xc_f03_func_get_info(initial_lda_func(jspins))))
            CALL xc_f03_lda_vxc(initial_lda_func(jspins), SIZE(rh,1), TRANSPOSE(rh), vx_tmp)
Matthias Redies's avatar
Matthias Redies committed
325
         endif
Matthias Redies's avatar
Matthias Redies committed
326 327
         IF (xcpot%func_vxc_id_c>0) THEN
            CALL xc_f03_lda_vxc(xcpot%vxc_func_c, SIZE(rh,1), TRANSPOSE(rh), vxc_tmp)
328
            write (*,*) "added correlation"
Matthias Redies's avatar
Matthias Redies committed
329 330 331 332 333 334
            vxc_tmp=vxc_tmp+vx_tmp
         ENDIF
      ENDIF
      vx=TRANSPOSE(vx_tmp) 
      vxc=TRANSPOSE(vxc_tmp) 

335
#endif
Matthias Redies's avatar
Matthias Redies committed
336 337
   END SUBROUTINE xcpot_get_vxc

338 339 340 341 342 343 344 345 346 347 348 349 350
   FUNCTION find_first_normal(vec) result(idx)
      USE, INTRINSIC :: IEEE_ARITHMETIC
      implicit none
      real, intent(in)    :: vec(:,:)
      integer             :: idx

      idx = size(vec, dim=2)
      do while(all(.not. ieee_is_nan(vec(:,idx))))
         idx = idx - 1
      enddo
   END FUNCTION find_first_normal


Matthias Redies's avatar
Matthias Redies committed
351
   SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad, kinED_KS, mt_call)
Matthias Redies's avatar
Matthias Redies committed
352
      use m_constants
Matthias Redies's avatar
Matthias Redies committed
353
      IMPLICIT NONE
354 355 356 357
   CLASS(t_xcpot_libxc),INTENT(IN)          :: xcpot
      INTEGER, INTENT (IN)                  :: jspins
      REAL,INTENT (IN)                      :: rh(:,:)  !points,spin
      REAL, INTENT (OUT)                    :: exc(:) !points
358
      ! optional arguments for GGA
359 360
      TYPE(t_gradients),OPTIONAL,INTENT(IN) :: grad
      LOGICAL, OPTIONAL, INTENT(IN)         :: mt_call    
361

362 363 364 365
      ! kinED from Kohn-Sham equations:
      ! tau = sum[phi_i(r)^dag nabla phi_i(r)]
      ! see eq (2) in https://doi.org/10.1063/1.1565316
      ! (-0.5 is applied below)
Matthias Redies's avatar
Matthias Redies committed
366
      REAL, INTENT(IN), OPTIONAL     :: kinED_KS(:,:)
Matthias Redies's avatar
Matthias Redies committed
367

368
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
369
      TYPE(xc_f03_func_info_t)       :: xc_info
Matthias Redies's avatar
Matthias Redies committed
370
      REAL  :: excc(SIZE(exc))
Matthias Redies's avatar
Matthias Redies committed
371
      REAL  :: cut_ratio = 0.1
Matthias Redies's avatar
Matthias Redies committed
372
      INTEGER :: cut_idx
373
      LOGICAL :: is_mt
Matthias Redies's avatar
Matthias Redies committed
374

Matthias Redies's avatar
Matthias Redies committed
375 376
      ! tau = 0.5 * sum[|grad phi_i(r)|²]
      ! see eq (3) in https://doi.org/10.1063/1.1565316
Matthias Redies's avatar
Matthias Redies committed
377
      REAL, ALLOCATABLE              :: kinED_libXC(:,:), pkzb_ratio(:,:), pkzb_zaehler(:,:), pkzb_nenner(:,:)
Matthias Redies's avatar
Matthias Redies committed
378

379
      is_mt = merge(mt_call, .False., present(mt_call))
380
      IF (xcpot%exc_is_gga()) THEN
Matthias Redies's avatar
Matthias Redies committed
381
         IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_exc for a GGA potential without providing derivatives")
Matthias Redies's avatar
Matthias Redies committed
382
         CALL xc_f03_gga_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,exc)
383
         IF (xcpot%func_exc_id_c>0) THEN
Matthias Redies's avatar
Matthias Redies committed
384 385 386
            CALL xc_f03_gga_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,excc)
            exc=exc+excc
         END IF
387
      ELSEIF(xcpot%exc_is_LDA()) THEN  !LDA potentials
Matthias Redies's avatar
Matthias Redies committed
388
         CALL xc_f03_lda_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh), exc)
389
         IF (xcpot%func_exc_id_c>0) THEN
Matthias Redies's avatar
Matthias Redies committed
390 391 392
            CALL xc_f03_lda_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh), excc)
            exc=exc+excc
         END IF
393
      ELSEIF(xcpot%exc_is_MetaGGA()) THEN
Matthias Redies's avatar
Matthias Redies committed
394
         IF(PRESENT(kinED_KS)) THEN 
395
            ! apply correction in  eq (4) in https://doi.org/10.1063/1.1565316
Matthias Redies's avatar
Matthias Redies committed
396
            kinED_libXC = transpose(kinED_KS + 0.25 * grad%laplace)
Matthias Redies's avatar
Matthias Redies committed
397

398 399 400
            !only cut core of muffin tin
            cut_idx = MERGE(NINT(size(rh,1) * cut_ratio), 0, is_mt)

Matthias Redies's avatar
Matthias Redies committed
401 402
            exc  = 0.0
            excc = 0.0
Matthias Redies's avatar
Matthias Redies committed
403 404 405 406
            call xc_f03_mgga_exc(xcpot%exc_func_x, SIZE(rh(cut_idx+1:,:),1),& 
                                                   TRANSPOSE(rh(cut_idx+1:,:)), &
                                                   grad%sigma(:,cut_idx+1:), &
                                                   transpose(grad%laplace(cut_idx+1:,:)), &
Matthias Redies's avatar
Matthias Redies committed
407
                                                   kinED_libXC(:,cut_idx+1:), &
Matthias Redies's avatar
Matthias Redies committed
408 409 410 411 412 413
                                                   exc(cut_idx+1:))

            call xc_f03_gga_exc(xcpot%vxc_func_x, SIZE(rh(:cut_idx,:),1), &
                                                  TRANSPOSE(rh(:cut_idx,:)), &
                                                  grad%sigma(:,:cut_idx), &
                                                  exc(:cut_idx))
Matthias Redies's avatar
Matthias Redies committed
414

Matthias Redies's avatar
Matthias Redies committed
415
            IF (xcpot%func_exc_id_c>0) THEN
Matthias Redies's avatar
Matthias Redies committed
416 417 418 419
               call xc_f03_mgga_exc(xcpot%exc_func_c, SIZE(rh(cut_idx+1:,:),1), &
                                                      TRANSPOSE(rh(cut_idx+1:,:)), &
                                                      grad%sigma(:,cut_idx+1:), &
                                                      transpose(grad%laplace(cut_idx+1:,:)), &
Matthias Redies's avatar
Matthias Redies committed
420
                                                      kinED_libXC(:,cut_idx+1:), &
Matthias Redies's avatar
Matthias Redies committed
421 422 423 424 425 426
                                                      excc(cut_idx+1:))
               
               call xc_f03_gga_exc(xcpot%vxc_func_c, SIZE(rh(:cut_idx,:),1), &
                                                      TRANSPOSE(rh(:cut_idx,:)), &
                                                      grad%sigma(:,:cut_idx), &
                                                      excc(:cut_idx))
Matthias Redies's avatar
Matthias Redies committed
427 428
               exc=exc+excc
            END IF
Matthias Redies's avatar
Matthias Redies committed
429 430 431 432 433 434 435 436
            
         ELSE ! first iteration is GGA
            IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_exc for a MetaGGA potential without providing derivatives")
            CALL xc_f03_gga_exc(xcpot%vxc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,exc)
            IF (xcpot%func_exc_id_c>0) THEN
               CALL xc_f03_gga_exc(xcpot%vxc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,excc)
               exc=exc+excc
            END IF
Matthias Redies's avatar
Matthias Redies committed
437
         ENDIF
Matthias Redies's avatar
Matthias Redies committed
438

439 440
      ELSE
         call juDFT_error("exc is part of a known Family", calledby="xcpot_get_exc@libxc")
Matthias Redies's avatar
Matthias Redies committed
441 442
      ENDIF

443 444
#endif
   END SUBROUTINE xcpot_get_exc
445

446 447 448 449 450 451 452 453 454
   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))
455

456 457 458 459
      grad%gr      = 0.0
      grad%laplace = 0.0
      grad%sigma   = 0.0
      grad%vsigma  = 0.0
460
   END SUBROUTINE xcpot_alloc_gradients
461

Matthias Redies's avatar
Matthias Redies committed
462
#ifdef CPP_LIBXC  
Matthias Redies's avatar
Matthias Redies committed
463 464 465 466 467 468 469
   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

470 471 472
      TYPE(xc_f03_func_t),INTENT(IN)      :: xc_func
      TYPE(xc_f03_func_info_t)            :: xc_info

Matthias Redies's avatar
Matthias Redies committed
473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514
      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
515
   END SUBROUTINE write_xc_info
516 517 518 519 520 521 522 523

   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

524
   function initial_lda_func(jspins) result(lda)
Matthias Redies's avatar
Matthias Redies committed
525 526 527 528
      use, intrinsic :: iso_c_binding
      implicit none
      INTEGER, intent(in) :: jspins
      TYPE(xc_f03_func_t) :: lda
529
      integer             :: ierr
Matthias Redies's avatar
Matthias Redies committed
530 531

      if(jspins == 1) then
532
         CALL xc_f03_func_init(lda, 1, XC_UNPOLARIZED, err=ierr)
Matthias Redies's avatar
Matthias Redies committed
533
      else
534
         CALL xc_f03_func_init(lda, 1, XC_POLARIZED, err=ierr)
Matthias Redies's avatar
Matthias Redies committed
535
      ENDIF
536 537 538 539

      if(ierr /= 0) call juDFT_error("can't find lda_x for init")

   end function initial_lda_func
Matthias Redies's avatar
Matthias Redies committed
540

Matthias Redies's avatar
Matthias Redies committed
541
#endif
542

543
END MODULE m_types_xcpot_libxc