types_xcpot_libxc.F90 13 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 12 13 14
#endif
  USE m_types_xcpot
  USE m_judft
  IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
15
  PRIVATE :: write_xc_info
16 17 18

  TYPE,EXTENDS(t_xcpot):: t_xcpot_libxc 
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
19 20
     TYPE(xc_f03_func_t)      :: vxc_func_x, vxc_func_c
     TYPE(xc_f03_func_t)      :: exc_func_x, exc_func_c
21
#endif
Matthias Redies's avatar
Matthias Redies committed
22 23 24
     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
25
   CONTAINS
Matthias Redies's avatar
Matthias Redies committed
26 27 28 29 30 31 32 33 34
     PROCEDURE        :: is_gga              => xcpot_is_gga
     PROCEDURE        :: is_MetaGGA          => xcpot_is_MetaGGA
     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
35 36 37 38
  END TYPE t_xcpot_libxc
  PUBLIC t_xcpot_libxc
CONTAINS

Matthias Redies's avatar
Matthias Redies committed
39
  SUBROUTINE xcpot_init(xcpot,jspins,vxc_id_x,vxc_id_c,exc_id_x,exc_id_c)
40 41 42
    USE m_judft
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(OUT)    :: xcpot
Matthias Redies's avatar
Matthias Redies committed
43 44 45 46
    INTEGER, INTENT(IN)                 :: jspins
    INTEGER, INTENT(IN)                 :: vxc_id_x, vxc_id_c ! potential functional
    INTEGER, INTENT(IN), OPTIONAL       :: exc_id_x, exc_id_c ! energy functionals
    LOGICAL                             :: same_functionals   ! are vxc and exc equal
47
    
48 49
#ifdef CPP_LIBXC
    INTEGER :: err
Matthias Redies's avatar
Matthias Redies committed
50 51 52 53 54 55 56 57 58 59 60 61
    xcpot%jspins = jspins

    xcpot%func_vxc_id_x = vxc_id_x
    xcpot%func_vxc_id_c = vxc_id_c

    xcpot%func_exc_id_x = merge(exc_id_x, vxc_id_x, PRESENT(exc_id_x))
    xcpot%func_exc_id_x = merge(exc_id_c, vxc_id_c, PRESENT(exc_id_c))

    same_functionals =     (xcpot%func_vxc_id_x == xcpot%func_exc_id_x) &
                     .and. (xcpot%func_vxc_id_c == xcpot%func_exc_id_c)

    if(xcpot%func_vxc_id_x == 0 .or. xcpot%func_exc_id_x == 0 ) then
62 63 64 65

      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) //& 
66
              '  <libXC  exchange="1" correlation="1" /> ' // ACHAR(10) //& 
67 68 69
              '</xcFunctional> ')
    endif 

70
    IF (jspins==1) THEN
Matthias Redies's avatar
Matthias Redies committed
71 72 73 74 75 76 77
      ! 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)
78
    ELSE
Matthias Redies's avatar
Matthias Redies committed
79 80 81 82 83 84 85
      ! 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)
86
    END IF
Matthias Redies's avatar
Matthias Redies committed
87 88 89 90 91

    CALL write_xc_info(xcpot%vxc_func_x)

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

Matthias Redies's avatar
Matthias Redies committed
96 97 98 99 100 101 102 103
    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
    END IF
104
#else
Matthias Redies's avatar
Matthias Redies committed
105 106
      CALL judft_error("You specified a libxc-exchange correlation potential but FLEUR is not linked against libxc", &
                       hint="Please recompile FLEUR with libxc support")
107
#endif
Matthias Redies's avatar
Matthias Redies committed
108
   END SUBROUTINE xcpot_init
109

Matthias Redies's avatar
Matthias Redies committed
110 111 112 113 114 115 116 117 118
   LOGICAL FUNCTION xcpot_is_gga(xcpot)
      IMPLICIT NONE
      CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC
      xcpot_is_gga=ANY((/XC_FAMILY_GGA, XC_FAMILY_HYB_GGA/)==xc_f03_func_info_get_family(xcpot%xc_info_x))
#else
      xcpot_is_gga=.false.
#endif
   END FUNCTION xcpot_is_gga
119

Matthias Redies's avatar
Matthias Redies committed
120 121 122 123 124
   LOGICAL FUNCTION xcpot_is_LDA(xcpot)
   IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC
   xcpot_is_LDA= (XC_FAMILY_LDA==xc_f03_func_info_get_family(xcpot%xc_info_x))
125
#else
Matthias Redies's avatar
Matthias Redies committed
126
   xcpot_is_LDA=.false.
Matthias Redies's avatar
Matthias Redies committed
127

Matthias Redies's avatar
Matthias Redies committed
128
   #endif
Matthias Redies's avatar
Matthias Redies committed
129
END FUNCTION xcpot_is_LDA
130

Matthias Redies's avatar
Matthias Redies committed
131 132 133 134
  LOGICAL FUNCTION xcpot_is_MetaGGA(xcpot)
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
135 136 137 138
    TYPE(xc_f03_func_info_t)        :: xc_info

    xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
    xcpot_is_MetaGGA=ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA]==xc_f03_func_info_get_family(xc_info))
139
#else
Matthias Redies's avatar
Matthias Redies committed
140
   xcpot_is_MetaGGA =  .False.
Matthias Redies's avatar
Matthias Redies committed
141
#endif
Matthias Redies's avatar
Matthias Redies committed
142
   END FUNCTION xcpot_is_MetaGGA
Matthias Redies's avatar
Matthias Redies committed
143

Matthias Redies's avatar
Matthias Redies committed
144 145 146
   LOGICAL FUNCTION xcpot_is_hybrid(xcpot)
      IMPLICIT NONE
      CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
147
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
148 149 150 151
    TYPE(xc_f03_func_info_t)        :: xc_info

    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))
152
#else
Matthias Redies's avatar
Matthias Redies committed
153
   xcpot_is_hybrid = .False.
154
#endif
Matthias Redies's avatar
Matthias Redies committed
155
   END FUNCTION xcpot_is_hybrid
156

Matthias Redies's avatar
Matthias Redies committed
157 158 159 160
   FUNCTION xcpot_get_exchange_weight(xcpot) RESULT(a_ex)
      USE m_judft
      IMPLICIT NONE
      CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
161 162 163

    REAL:: a_ex
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
164
    a_ex=xc_f03_hyb_exx_coef(xcpot%vxc_func_x)
165 166 167 168
#endif
  END FUNCTION xcpot_get_exchange_weight

  !***********************************************************************
169
  SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad)
170 171 172 173
    !***********************************************************************
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
    INTEGER, INTENT (IN)     :: jspins
174 175 176
    REAL,INTENT (IN)         :: rh(:,:)   !Dimensions here
    REAL, INTENT (OUT)       :: vx (:,:)  !points,spin
    REAL, INTENT (OUT  )     :: vxc(:,:)  !
177
    ! optional arguments for GGA
178
    TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
179
#ifdef CPP_LIBXC    
180
    REAL,ALLOCATABLE::vxc_tmp(:,:),vx_tmp(:,:),vsigma(:,:)
181
    !libxc uses the spin as a first index, hence we have to transpose....
182 183
    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
184 185
    IF (xcpot%is_gga()) THEN
       IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives")
Daniel Wortmann's avatar
Daniel Wortmann committed
186 187
       ALLOCATE(vsigma,mold=grad%vsigma)
       !where(abs(grad%sigma)<1E-9) grad%sigma=1E-9
Matthias Redies's avatar
Matthias Redies committed
188 189 190
       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)
191
          grad%vsigma=grad%vsigma+vsigma
192
          vxc_tmp=vxc_tmp+vx_tmp
Daniel Wortmann's avatar
Daniel Wortmann committed
193 194
       ELSE     
         vxc_tmp=vx_tmp
195
       ENDIF
196
    ELSE  !LDA potentials
Matthias Redies's avatar
Matthias Redies committed
197 198 199
       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)
200 201
          vxc_tmp=vxc_tmp+vx_tmp
       ENDIF
202
    ENDIF
203
    vx=TRANSPOSE(vx_tmp) 
204
    vxc=TRANSPOSE(vxc_tmp) 
205
    
206 207 208
#endif
  END SUBROUTINE xcpot_get_vxc

209
  
210 211 212 213 214 215
  !***********************************************************************
  SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad)
    !***********************************************************************
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
    INTEGER, INTENT (IN)     :: jspins
216 217
    REAL,INTENT (IN) :: rh(:,:)  !points,spin
    REAL, INTENT (OUT) :: exc(:) !points
218 219
    ! optional arguments for GGA      
    TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad
220 221

    REAL  :: excc(SIZE(exc))
222 223 224
#ifdef CPP_LIBXC    
    IF (xcpot%is_gga()) THEN
       IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives")
Matthias Redies's avatar
Matthias Redies committed
225 226 227
       CALL xc_f03_gga_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,exc)
       IF (xcpot%func_vxc_id_c>0) THEN
          CALL xc_f03_gga_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,excc)
228 229
          exc=exc+excc
       END IF
230
    ELSE  !LDA potentials
Matthias Redies's avatar
Matthias Redies committed
231 232 233
       CALL xc_f03_lda_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh), exc)
       IF (xcpot%func_vxc_id_c>0) THEN
          CALL xc_f03_lda_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh), excc)
234 235
          exc=exc+excc
       END IF
236
    ENDIF
237
    
238
#endif
Matthias Redies's avatar
Matthias Redies committed
239
   END SUBROUTINE xcpot_get_vxc
240

Matthias Redies's avatar
Matthias Redies committed
241 242 243 244 245 246 247 248 249 250
   !***********************************************************************
   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
      ! optional arguments for GGA
      TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad
251

Matthias Redies's avatar
Matthias Redies committed
252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267
      REAL  :: excc(SIZE(exc))
#ifdef CPP_LIBXC
      IF (xcpot%is_gga()) THEN
         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%xc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,exc)
         IF (xcpot%func_id_c>0) THEN
            CALL xc_f03_gga_exc(xcpot%xc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,excc)
            exc=exc+excc
         END IF
      ELSE  !LDA potentials
         CALL xc_f03_lda_exc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh), exc)
         IF (xcpot%func_id_c>0) THEN
            CALL xc_f03_lda_exc(xcpot%xc_func_c, SIZE(rh,1), TRANSPOSE(rh), excc)
            exc=exc+excc
         END IF
      ENDIF
268

Matthias Redies's avatar
Matthias Redies committed
269 270
#endif
   END SUBROUTINE xcpot_get_exc
271

Matthias Redies's avatar
Matthias Redies committed
272 273 274 275 276 277 278 279 280
   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))
281

Matthias Redies's avatar
Matthias Redies committed
282
   END SUBROUTINE xcpot_alloc_gradients
283

Matthias Redies's avatar
Matthias Redies committed
284 285
#ifdef CPP_LIBXC  
  SUBROUTINE write_xc_info(xc_func, is_E_func)
286
    IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
287 288 289 290 291 292
    TYPE(xc_f03_func_t),INTENT(IN)      :: xc_func
    LOGICAL,INTENT(IN),OPTIONAL         :: is_E_func
    TYPE(xc_f03_func_info_t)            :: xc_info
    INTEGER                             :: i
    CHARACTER(len=120)                  :: kind, family
    LOGICAL                             :: is_energy_func
293
    
Matthias Redies's avatar
Matthias Redies committed
294 295 296
    xc_info = xc_f03_func_get_info(xc_func)
    is_energy_func = merge(is_E_func, .False., PRESENT(is_E_func))

297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323
    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
    
Matthias Redies's avatar
Matthias Redies committed
324 325 326 327 328 329 330 331
    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

332 333 334 335
    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
336
  END SUBROUTINE write_xc_info
Daniel Wortmann's avatar
Daniel Wortmann committed
337
#endif
Matthias Redies's avatar
Matthias Redies committed
338

339
END MODULE m_types_xcpot_libxc