types_xcpot_libxc.F90 13.3 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 14 15
   USE m_types_xcpot
   USE m_judft
   IMPLICIT NONE
   PRIVATE :: write_xc_info
16

Matthias Redies's avatar
Matthias Redies committed
17
   TYPE,EXTENDS(t_xcpot):: t_xcpot_libxc 
18
#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 35 36
      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
   END TYPE t_xcpot_libxc
   PUBLIC t_xcpot_libxc
37 38
CONTAINS

Matthias Redies's avatar
Matthias Redies committed
39 40 41 42 43 44 45 46 47
   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
      INTEGER, INTENT(IN), OPTIONAL       :: exc_id_x, exc_id_c ! energy functionals
      LOGICAL                             :: same_functionals   ! are vxc and exc equal

48
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
      INTEGER :: err
      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

         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
78
      ELSE
Matthias Redies's avatar
Matthias Redies committed
79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
         ! 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

      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

      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
      CALL judft_error("You specified a libxc-exchange correlation potential but FLEUR is not linked against libxc", &
Matthias Redies's avatar
Matthias Redies committed
106
         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
   LOGICAL FUNCTION xcpot_is_gga(xcpot)
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
112
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
Matthias Redies's avatar
Matthias Redies committed
113 114 115 116 117 118
#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
   LOGICAL FUNCTION xcpot_is_LDA(xcpot)
Matthias Redies's avatar
Matthias Redies committed
121
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
122 123
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
124
      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 129
#endif
   END FUNCTION xcpot_is_LDA
130

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

Matthias Redies's avatar
Matthias Redies committed
137 138
      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
   LOGICAL FUNCTION xcpot_is_hybrid(xcpot)
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
146
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
147
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
148
      TYPE(xc_f03_func_info_t)        :: xc_info
Matthias Redies's avatar
Matthias Redies committed
149

Matthias Redies's avatar
Matthias Redies committed
150 151
      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
   FUNCTION xcpot_get_exchange_weight(xcpot) RESULT(a_ex)
      USE m_judft
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
160
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
161

Matthias Redies's avatar
Matthias Redies committed
162
      REAL:: a_ex
163
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
164
      a_ex=xc_f03_hyb_exx_coef(xcpot%vxc_func_x)
165
#endif
Matthias Redies's avatar
Matthias Redies committed
166 167 168 169 170 171 172 173 174 175 176 177 178
   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
179
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
      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
      IF (xcpot%is_gga()) THEN
         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) 

206
#endif
Matthias Redies's avatar
Matthias Redies committed
207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
   END SUBROUTINE xcpot_get_vxc


   !***********************************************************************
   SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad)
      !***********************************************************************
      IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
      INTEGER, INTENT (IN)     :: jspins
      REAL,INTENT (IN) :: rh(:,:)  !points,spin
      REAL, INTENT (OUT) :: exc(:) !points
      ! optional arguments for GGA      
      TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad

      REAL  :: excc(SIZE(exc))
222
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
223 224 225 226 227 228 229 230 231 232 233 234 235 236 237
      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%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)
            exc=exc+excc
         END IF
      ELSE  !LDA potentials
         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)
            exc=exc+excc
         END IF
      ENDIF

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
   !***********************************************************************
   SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad)
      !***********************************************************************
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
245
   CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
Matthias Redies's avatar
Matthias Redies committed
246 247 248 249 250
      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
Matthias Redies's avatar
Matthias Redies committed
283
#ifdef CPP_LIBXC  
Matthias Redies's avatar
Matthias Redies committed
284 285 286 287 288 289 290 291 292 293 294 295 296 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 324 325 326 327 328 329 330 331 332 333 334 335
   SUBROUTINE write_xc_info(xc_func, is_E_func)
      IMPLICIT NONE
      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

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

338
END MODULE m_types_xcpot_libxc