types_xcpot_libxc.F90 21.1 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
Matthias Redies's avatar
Matthias Redies committed
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
      LOGICAL                             :: same_functionals   ! are vxc and exc equal
Matthias Redies's avatar
Matthias Redies committed
61
      INTEGER                             :: errors(4)
Matthias Redies's avatar
Matthias Redies committed
62

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

      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      
Matthias Redies's avatar
Matthias Redies committed
83 84 85
         CALL xc_f03_func_init(xcpot%vxc_func_x, xcpot%func_vxc_id_x, XC_UNPOLARIZED, err=errors(1))
         IF (xcpot%func_vxc_id_c>0) CALL xc_f03_func_init(xcpot%vxc_func_c, xcpot%func_vxc_id_c, &
                                                                 XC_UNPOLARIZED, err=errors(2))
Matthias Redies's avatar
Matthias Redies committed
86 87

         ! energy functionals
Matthias Redies's avatar
Matthias Redies committed
88 89 90
         CALL xc_f03_func_init(xcpot%exc_func_x, xcpot%func_exc_id_x, XC_UNPOLARIZED, err=errors(3))
         IF (xcpot%func_exc_id_c>0) CALL xc_f03_func_init(xcpot%exc_func_c, xcpot%func_exc_id_c, &
                                                                  XC_UNPOLARIZED, err=errors(4))
Matthias Redies's avatar
Matthias Redies committed
91
         
Matthias Redies's avatar
Matthias Redies committed
92
      ELSE
Matthias Redies's avatar
Matthias Redies committed
93
         ! potential functionals  
Matthias Redies's avatar
Matthias Redies committed
94 95 96
         CALL xc_f03_func_init(xcpot%vxc_func_x, xcpot%func_vxc_id_x, XC_POLARIZED, err=errors(1))
         IF (xcpot%func_vxc_id_c>0) CALL xc_f03_func_init(xcpot%vxc_func_c, xcpot%func_vxc_id_c, &
                                                                  XC_POLARIZED, err=errors(2))
Matthias Redies's avatar
Matthias Redies committed
97 98

         !energy functionals
Matthias Redies's avatar
Matthias Redies committed
99 100 101
         CALL xc_f03_func_init(xcpot%exc_func_x, xcpot%func_exc_id_x, XC_POLARIZED, err=errors(3))
         IF (xcpot%func_exc_id_c>0) CALL xc_f03_func_init(xcpot%exc_func_c, xcpot%func_exc_id_c, &
                                                                  XC_POLARIZED, err=errors(4))
Matthias Redies's avatar
Matthias Redies committed
102 103
      END IF

Matthias Redies's avatar
Matthias Redies committed
104 105 106 107 108
      IF(errors(1) /= 0) call juDFT_error("Exchange potential functional not in LibXC")
      IF(errors(2) /= 0) call juDFT_error("Correlation potential functional not in LibXC")
      IF(errors(3) /= 0) call juDFT_error("Exchange energy functional not in LibXC")
      IF(errors(4) /= 0) call juDFT_error("Correlation energy functional not in LibXC")

109
      !check if any potental is a MetaGGA
Matthias Redies's avatar
Matthias Redies committed
110 111
      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")
112 113
      ENDIF

Matthias Redies's avatar
Matthias Redies committed
114 115 116 117 118 119 120 121
      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

122 123
      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
124 125 126 127 128 129 130
      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
131 132
      ELSE
         write (*,*) "Using same functional for VXC and EXC"
Matthias Redies's avatar
Matthias Redies committed
133
      END IF
Matthias Redies's avatar
Matthias Redies committed
134 135

      xcpot%kinED%set = .False.
136
#else
Matthias Redies's avatar
Matthias Redies committed
137
      CALL judft_error("You specified a libxc-exchange correlation potential but FLEUR is not linked against libxc", &
Matthias Redies's avatar
Matthias Redies committed
138
         hint="Please recompile FLEUR with libxc support")
139
#endif
Matthias Redies's avatar
Matthias Redies committed
140
   END SUBROUTINE xcpot_init
141

Matthias Redies's avatar
Matthias Redies committed
142
   ! LDA
143
   LOGICAL FUNCTION xcpot_vx_is_LDA(xcpot)
Matthias Redies's avatar
Matthias Redies committed
144 145 146 147 148 149
      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)
150
      xcpot_vx_is_LDA =  XC_FAMILY_LDA == xc_f03_func_info_get_family(xc_info)
Matthias Redies's avatar
Matthias Redies committed
151
#else
152
      xcpot_vx_is_LDA = .false.
Matthias Redies's avatar
Matthias Redies committed
153
#endif
154 155 156 157 158 159 160
   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
161

162 163 164 165 166 167 168
      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
169 170 171 172 173 174
   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
175 176
      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
177 178 179 180 181 182
#else
      xcpot_exc_is_LDA = .false.
#endif
   END FUNCTION xcpot_exc_is_LDA

   ! GGA
183
   LOGICAL FUNCTION xcpot_vc_is_gga(xcpot)
Matthias Redies's avatar
Matthias Redies committed
184
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
185
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
Matthias Redies's avatar
Matthias Redies committed
186
#ifdef CPP_LIBXC
187 188
      TYPE(xc_f03_func_info_t)        :: xc_info

189 190 191 192 193 194 195 196 197 198 199 200 201 202 203
      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))
204
#else
205
      xcpot_vx_is_gga=.false.
206
#endif
207
   END FUNCTION xcpot_vx_is_gga
208

Matthias Redies's avatar
Matthias Redies committed
209 210 211 212 213 214
   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
215
      xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
Matthias Redies's avatar
Matthias Redies committed
216 217 218 219 220 221
      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

222 223 224 225 226 227
   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
228
      xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
Matthias Redies's avatar
Matthias Redies committed
229
      xcpot_exc_is_gga =  ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f03_func_info_get_family(xc_info))
Matthias Redies's avatar
Matthias Redies committed
230
#else
231
      xcpot_exc_is_gga=.false.
Matthias Redies's avatar
Matthias Redies committed
232
#endif
233
   END FUNCTION xcpot_exc_is_gga
234

Matthias Redies's avatar
Matthias Redies committed
235
   LOGICAL FUNCTION xcpot_exc_is_MetaGGA(xcpot)
Matthias Redies's avatar
Matthias Redies committed
236 237
      IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
Matthias Redies's avatar
Matthias Redies committed
238
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
239
      TYPE(xc_f03_func_info_t)        :: xc_info
Matthias Redies's avatar
Matthias Redies committed
240

Matthias Redies's avatar
Matthias Redies committed
241
      xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
Matthias Redies's avatar
Matthias Redies committed
242
      xcpot_exc_is_MetaGGA=ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA]==xc_f03_func_info_get_family(xc_info))
243
#else
Matthias Redies's avatar
Matthias Redies committed
244
      xcpot_exc_is_MetaGGA =  .False.
Matthias Redies's avatar
Matthias Redies committed
245
#endif
Matthias Redies's avatar
Matthias Redies committed
246
   END FUNCTION xcpot_exc_is_MetaGGA
Matthias Redies's avatar
Matthias Redies committed
247

Matthias Redies's avatar
Matthias Redies committed
248 249
   LOGICAL FUNCTION xcpot_is_hybrid(xcpot)
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
250
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
251
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
252
      TYPE(xc_f03_func_info_t)        :: xc_info
Matthias Redies's avatar
Matthias Redies committed
253

Matthias Redies's avatar
Matthias Redies committed
254 255
      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))
256
#else
Matthias Redies's avatar
Matthias Redies committed
257
      xcpot_is_hybrid = .False.
258
#endif
Matthias Redies's avatar
Matthias Redies committed
259
   END FUNCTION xcpot_is_hybrid
260

Matthias Redies's avatar
Matthias Redies committed
261 262 263
   FUNCTION xcpot_get_exchange_weight(xcpot) RESULT(a_ex)
      USE m_judft
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
264
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
265

Matthias Redies's avatar
Matthias Redies committed
266
      REAL:: a_ex
267
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
268
      a_ex=xc_f03_hyb_exx_coef(xcpot%vxc_func_x)
269
#endif
Matthias Redies's avatar
Matthias Redies committed
270 271 272
   END FUNCTION xcpot_get_exchange_weight

   !***********************************************************************
Matthias Redies's avatar
Matthias Redies committed
273
   SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad, kinED_KS)
274
      USE, INTRINSIC :: IEEE_ARITHMETIC
Matthias Redies's avatar
Matthias Redies committed
275 276 277 278 279 280
      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
281
      REAL, INTENT(IN),OPTIONAL:: kinED_KS(:,:)
Matthias Redies's avatar
Matthias Redies committed
282 283
      ! optional arguments for GGA
      TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
284
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
285 286 287
      REAL,ALLOCATABLE  :: vxc_tmp(:,:),vx_tmp(:,:),vsigma(:,:), &
                           tmp_vsig(:,:), tmp_vlapl(:,:), tmp_vtau(:,:), &
                           kinED_libxc(:,:)
288
      integer           :: idx
Matthias Redies's avatar
Matthias Redies committed
289 290 291
      !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
292

Matthias Redies's avatar
Matthias Redies committed
293
      IF (xcpot%vc_is_GGA()) THEN
Matthias Redies's avatar
Matthias Redies committed
294 295 296 297 298 299 300 301 302 303 304 305
         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
306 307 308 309 310 311 312 313 314 315 316 317 318 319
         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)
            
            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))
Matthias Redies's avatar
Matthias Redies committed
320
         else
321
            CALL xc_f03_lda_vxc(initial_lda_func(jspins), SIZE(rh,1), TRANSPOSE(rh), vx_tmp)
Matthias Redies's avatar
Matthias Redies committed
322
         endif
Matthias Redies's avatar
Matthias Redies committed
323 324 325 326 327 328 329 330
         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) 

331
#endif
Matthias Redies's avatar
Matthias Redies committed
332 333
   END SUBROUTINE xcpot_get_vxc

334 335 336 337 338 339 340 341 342
   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
343
         if(idx == 0) exit
344 345 346 347
      enddo
   END FUNCTION find_first_normal


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

Matthias Redies's avatar
Matthias Redies committed
359 360 361 362
      ! 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
363
      REAL, INTENT(IN), OPTIONAL     :: kinED_KS(:,:)
Matthias Redies's avatar
Matthias Redies committed
364

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

Matthias Redies's avatar
Matthias Redies committed
372 373
      ! 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
374
      REAL, ALLOCATABLE              :: kinED_libXC(:,:), pkzb_ratio(:,:), pkzb_zaehler(:,:), pkzb_nenner(:,:)
Matthias Redies's avatar
Matthias Redies committed
375

376
      is_mt = merge(mt_call, .False., present(mt_call))
377
      IF (xcpot%exc_is_gga()) THEN
Matthias Redies's avatar
Matthias Redies committed
378
         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
379
         CALL xc_f03_gga_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,exc)
380
         IF (xcpot%func_exc_id_c>0) THEN
Matthias Redies's avatar
Matthias Redies committed
381 382 383
            CALL xc_f03_gga_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,excc)
            exc=exc+excc
         END IF
384
      ELSEIF(xcpot%exc_is_LDA()) THEN  !LDA potentials
Matthias Redies's avatar
Matthias Redies committed
385
         CALL xc_f03_lda_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh), exc)
386
         IF (xcpot%func_exc_id_c>0) THEN
Matthias Redies's avatar
Matthias Redies committed
387 388 389
            CALL xc_f03_lda_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh), excc)
            exc=exc+excc
         END IF
390
      ELSEIF(xcpot%exc_is_MetaGGA()) THEN
Matthias Redies's avatar
Matthias Redies committed
391
         IF(PRESENT(kinED_KS)) THEN 
Matthias Redies's avatar
Matthias Redies committed
392
            ! apply correction in  eq (4) in https://doi.org/10.1063/1.1565316
Matthias Redies's avatar
Matthias Redies committed
393
            kinED_libXC = transpose(kinED_KS + 0.25 * grad%laplace)
Matthias Redies's avatar
Matthias Redies committed
394

395 396 397
            !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
398 399
            exc  = 0.0
            excc = 0.0
Matthias Redies's avatar
Matthias Redies committed
400 401 402 403
            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
404
                                                   kinED_libXC(:,cut_idx+1:), &
Matthias Redies's avatar
Matthias Redies committed
405 406 407 408 409 410
                                                   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
cut off  
Matthias Redies committed
411

Matthias Redies's avatar
Matthias Redies committed
412
            IF (xcpot%func_exc_id_c>0) THEN
Matthias Redies's avatar
Matthias Redies committed
413 414 415 416
               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
417
                                                      kinED_libXC(:,cut_idx+1:), &
Matthias Redies's avatar
Matthias Redies committed
418 419 420 421 422 423
                                                      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
424 425
               exc=exc+excc
            END IF
Matthias Redies's avatar
Matthias Redies committed
426 427 428 429 430 431 432 433
            
         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
434
         ENDIF
Matthias Redies's avatar
Matthias Redies committed
435

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

Matthias Redies's avatar
Matthias Redies committed
440 441
#endif
   END SUBROUTINE xcpot_get_exc
442

Matthias Redies's avatar
Matthias Redies committed
443 444 445 446 447 448 449 450 451
   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))
452

453 454 455 456
      grad%gr      = 0.0
      grad%laplace = 0.0
      grad%sigma   = 0.0
      grad%vsigma  = 0.0
Matthias Redies's avatar
Matthias Redies committed
457
   END SUBROUTINE xcpot_alloc_gradients
Matthias Redies's avatar
Matthias Redies committed
458

Matthias Redies's avatar
Matthias Redies committed
459
#ifdef CPP_LIBXC  
Matthias Redies's avatar
Matthias Redies committed
460 461 462 463 464 465 466
   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

Matthias Redies's avatar
Matthias Redies committed
467 468 469
      TYPE(xc_f03_func_t),INTENT(IN)      :: xc_func
      TYPE(xc_f03_func_info_t)            :: xc_info

Matthias Redies's avatar
Matthias Redies committed
470 471 472 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
      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
Matthias Redies's avatar
Matthias Redies committed
512
   END SUBROUTINE write_xc_info
513 514 515 516 517 518 519 520

   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

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

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

      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
537

Matthias Redies's avatar
Matthias Redies committed
538
#endif
Matthias Redies's avatar
Matthias Redies committed
539

540
END MODULE m_types_xcpot_libxc