types_xcpot_libxc.F90 17.4 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 30 31
      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


32
   CONTAINS
Matthias Redies's avatar
Matthias Redies committed
33 34
      PROCEDURE        :: vxc_is_LDA          => xcpot_vxc_is_LDA
      PROCEDURE        :: exc_is_LDA          => xcpot_exc_is_LDA
35 36
      PROCEDURE        :: vxc_is_gga          => xcpot_vxc_is_gga
      PROCEDURE        :: exc_is_gga          => xcpot_exc_is_gga
Matthias Redies's avatar
Matthias Redies committed
37
      PROCEDURE        :: exc_is_MetaGGA      => xcpot_exc_is_MetaGGA
Matthias Redies's avatar
Matthias Redies committed
38 39 40 41 42 43 44 45 46
      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
47 48
CONTAINS

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

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

      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
84
         
Matthias Redies's avatar
Matthias Redies committed
85
      ELSE
Matthias Redies's avatar
Matthias Redies committed
86 87 88 89 90 91 92 93 94
         ! 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

95
      !check if any potental is a MetaGGA
Matthias Redies's avatar
Matthias Redies committed
96
      IF(ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA] == xc_get_family(xcpot%vxc_func_x))) THEN
97
         call juDFT_error("vxc_x: MetaGGA is not implemented for potentials")
Matthias Redies's avatar
Matthias Redies committed
98 99
      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
100 101 102 103
            call juDFT_error("vxc_x: MetaGGA is not implemented for potentials")
         ENDIF
      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
124
#else
Matthias Redies's avatar
Matthias Redies committed
125
      CALL judft_error("You specified a libxc-exchange correlation potential but FLEUR is not linked against libxc", &
Matthias Redies's avatar
Matthias Redies committed
126
         hint="Please recompile FLEUR with libxc support")
127
#endif
Matthias Redies's avatar
Matthias Redies committed
128
   END SUBROUTINE xcpot_init
129

Matthias Redies's avatar
Matthias Redies committed
130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
   ! 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
158
   LOGICAL FUNCTION xcpot_vxc_is_gga(xcpot)
Matthias Redies's avatar
Matthias Redies committed
159
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
160
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
Matthias Redies's avatar
Matthias Redies committed
161
#ifdef CPP_LIBXC
162 163
      TYPE(xc_f03_func_info_t)        :: xc_info

164
      xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
Matthias Redies's avatar
Matthias Redies committed
165
      xcpot_vxc_is_gga =  ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f03_func_info_get_family(xc_info))
166 167 168 169 170 171 172 173 174 175 176
#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

Matthias Redies's avatar
Matthias Redies committed
177
      xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
Matthias Redies's avatar
Matthias Redies committed
178
      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
179
#else
180
      xcpot_exc_is_gga=.false.
Matthias Redies's avatar
Matthias Redies committed
181
#endif
182
   END FUNCTION xcpot_exc_is_gga
183

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

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

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

Matthias Redies's avatar
Matthias Redies committed
203 204
      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))
205
#else
Matthias Redies's avatar
Matthias Redies committed
206
      xcpot_is_hybrid = .False.
207
#endif
Matthias Redies's avatar
Matthias Redies committed
208
   END FUNCTION xcpot_is_hybrid
209

Matthias Redies's avatar
Matthias Redies committed
210 211 212
   FUNCTION xcpot_get_exchange_weight(xcpot) RESULT(a_ex)
      USE m_judft
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
213
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
214

Matthias Redies's avatar
Matthias Redies committed
215
      REAL:: a_ex
216
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
217
      a_ex=xc_f03_hyb_exx_coef(xcpot%vxc_func_x)
218
#endif
Matthias Redies's avatar
Matthias Redies committed
219 220 221 222 223 224 225 226 227 228 229 230 231
   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
232
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
233 234 235 236
      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
Matthias Redies's avatar
Matthias Redies committed
237

238
      IF (xcpot%needs_grad()) THEN
Matthias Redies's avatar
Matthias Redies committed
239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259
         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) 

260
#endif
Matthias Redies's avatar
Matthias Redies committed
261 262 263
   END SUBROUTINE xcpot_get_vxc


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

Matthias Redies's avatar
Matthias Redies committed
274 275 276 277 278
      ! 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)
      REAL, INTENT(IN), OPTIONAL     :: kinEnergyDen_KS(:,:)
Matthias Redies's avatar
Matthias Redies committed
279
      character(len=200)             :: filename
Matthias Redies's avatar
Matthias Redies committed
280

281
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
282
      TYPE(xc_f03_func_info_t)       :: xc_info
Matthias Redies's avatar
Matthias Redies committed
283
      REAL  :: excc(SIZE(exc))
Matthias Redies's avatar
Matthias Redies committed
284
      REAL  :: cut_ratio = 0.1
Matthias Redies's avatar
Matthias Redies committed
285
      INTEGER :: cut_idx
Matthias Redies's avatar
Matthias Redies committed
286

Matthias Redies's avatar
Matthias Redies committed
287 288
      ! 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
289
      REAL, ALLOCATABLE              :: kinEnergyDen_libXC(:,:), pkzb_ratio(:,:), pkzb_zaehler(:,:), pkzb_nenner(:,:)
Matthias Redies's avatar
Matthias Redies committed
290

Matthias Redies's avatar
Matthias Redies committed
291

292
      IF (xcpot%exc_is_gga()) THEN
Matthias Redies's avatar
Matthias Redies committed
293
         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
294
         CALL xc_f03_gga_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,exc)
295
         IF (xcpot%func_exc_id_c>0) THEN
Matthias Redies's avatar
Matthias Redies committed
296 297 298
            CALL xc_f03_gga_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,excc)
            exc=exc+excc
         END IF
299
      ELSEIF(xcpot%exc_is_LDA()) THEN  !LDA potentials
Matthias Redies's avatar
Matthias Redies committed
300
         CALL xc_f03_lda_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh), exc)
301
         IF (xcpot%func_exc_id_c>0) THEN
Matthias Redies's avatar
Matthias Redies committed
302 303 304
            CALL xc_f03_lda_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh), excc)
            exc=exc+excc
         END IF
305
      ELSEIF(xcpot%exc_is_MetaGGA()) THEN
Matthias Redies's avatar
Matthias Redies committed
306
         IF(PRESENT(kinEnergyDen_KS)) THEN 
Matthias Redies's avatar
Matthias Redies committed
307
            cut_idx = NINT(size(rh,1) * cut_ratio)
Matthias Redies's avatar
Matthias Redies committed
308
            ! apply correction in  eq (4) in https://doi.org/10.1063/1.1565316
Matthias Redies's avatar
merge  
Matthias Redies committed
309
            kinEnergyDen_libXC = transpose(kinEnergyDen_KS + 0.25 * grad%laplace)
Matthias Redies's avatar
Matthias Redies committed
310 311 312

            exc  = 0.0
            excc = 0.0
Matthias Redies's avatar
Matthias Redies committed
313 314 315 316 317 318 319 320 321 322 323
            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:,:)), &
                                                   kinEnergyDen_libXC(:,cut_idx+1:), &
                                                   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
324

Matthias Redies's avatar
Matthias Redies committed
325
            IF (xcpot%func_exc_id_c>0) THEN
Matthias Redies's avatar
Matthias Redies committed
326 327 328 329 330 331 332 333 334 335 336 337 338
               !CALL xc_f03_mgga_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh), grad%sigma, &
                                    !transpose(grad%laplace), kinEnergyDen_libXC, excc)
               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:,:)), &
                                                      kinEnergyDen_libXC(:,cut_idx+1:), &
                                                      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
339 340
               exc=exc+excc
            END IF
Matthias Redies's avatar
Matthias Redies committed
341 342 343 344 345 346 347 348
            
         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
349
         ENDIF
Matthias Redies's avatar
Matthias Redies committed
350

351 352
      ELSE
         call juDFT_error("exc is part of a known Family", calledby="xcpot_get_exc@libxc")
Matthias Redies's avatar
Matthias Redies committed
353 354
      ENDIF

Matthias Redies's avatar
Matthias Redies committed
355 356
#endif
   END SUBROUTINE xcpot_get_exc
357

Matthias Redies's avatar
Matthias Redies committed
358 359 360 361 362 363 364 365 366
   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))
367

Matthias Redies's avatar
Matthias Redies committed
368
   END SUBROUTINE xcpot_alloc_gradients
Matthias Redies's avatar
Matthias Redies committed
369

Matthias Redies's avatar
Matthias Redies committed
370
#ifdef CPP_LIBXC  
Matthias Redies's avatar
Matthias Redies committed
371 372 373 374 375 376 377
   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
378 379 380
      TYPE(xc_f03_func_t),INTENT(IN)      :: xc_func
      TYPE(xc_f03_func_info_t)            :: xc_info

Matthias Redies's avatar
Matthias Redies committed
381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422
      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
423
   END SUBROUTINE write_xc_info
424 425 426 427 428 429 430 431

   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
432
#endif
Matthias Redies's avatar
Matthias Redies committed
433

434
END MODULE m_types_xcpot_libxc