types_xcpot_libxc.F90 9.07 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
!--------------------------------------------------------------------------------
! 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
#ifdef CPP_LIBXC  
  USE xc_f03_lib_m
#endif
  USE m_types_xcpot
  USE m_judft
  IMPLICIT NONE
  PRIVATE

  TYPE,EXTENDS(t_xcpot):: t_xcpot_libxc 
#ifdef CPP_LIBXC    
19 20
     TYPE(xc_f03_func_t) :: xc_func_x,xc_func_c
     TYPE(xc_f03_func_info_t) :: xc_info_x,xc_info_c
21
#endif
22
     INTEGER          :: func_id_c,func_id_x,jspins
23 24
   CONTAINS
     PROCEDURE        :: is_gga=>xcpot_is_gga
Matthias Redies's avatar
Matthias Redies committed
25
     PROCEDURE        :: is_MetaGGA=>xcpot_is_MetaGGA
26 27 28 29 30 31 32 33 34 35 36
     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
CONTAINS

37
  SUBROUTINE xcpot_init(xcpot,jspins,id_x,id_c)
38 39 40
    USE m_judft
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(OUT)    :: xcpot
41 42
    INTEGER,INTENT(IN)                  :: jspins,id_x,id_c
    
43 44
#ifdef CPP_LIBXC
    INTEGER :: err
45 46 47
    xcpot%jspins=jspins
    xcpot%func_id_x=id_x
    xcpot%func_id_c=id_c
48 49 50 51 52

    if(xcpot%func_id_x == 0 .or. xcpot%func_id_c == 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) //& 
53
              '  <libXC  exchange="1" correlation="1" /> ' // ACHAR(10) //& 
54 55 56
              '</xcFunctional> ')
    endif 

57
    IF (jspins==1) THEN
58 59
       CALL xc_f03_func_init(xcpot%xc_func_x, xcpot%func_id_x, XC_UNPOLARIZED)
       IF (xcpot%func_id_c>0) CALL xc_f03_func_init(xcpot%xc_func_c, xcpot%func_id_c, XC_UNPOLARIZED)
60
    ELSE
61 62
       CALL xc_f03_func_init(xcpot%xc_func_x, xcpot%func_id_x, XC_POLARIZED)
       IF (xcpot%func_id_c>0) CALL xc_f03_func_init(xcpot%xc_func_c, xcpot%func_id_c, XC_POLARIZED)
63
    END IF
64
    xcpot%xc_info_x=xc_f03_func_get_info(xcpot%xc_func_x)
65 66 67 68 69 70 71
    CALL priv_write_info(xcpot%xc_info_x)
    IF (xcpot%func_id_c>0) THEN
       xcpot%xc_info_c=xc_f03_func_get_info(xcpot%xc_func_c)
       CALL priv_write_info(xcpot%xc_info_c)
    ELSE
       WRITE(*,*) "No Correlation functional"
    END IF
72 73

#else
Matthias Redies's avatar
Matthias Redies committed
74 75
    CALL judft_error("You specified a libxc-exchange correlation potential but FLEUR is not linked against libxc", &
                      hint="Please recompile FLEUR with libxc support")
76 77 78 79 80 81 82 83
#endif
  END SUBROUTINE xcpot_init


  LOGICAL FUNCTION xcpot_is_gga(xcpot)
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC    
84
    xcpot_is_gga=ANY((/XC_FAMILY_GGA, XC_FAMILY_HYB_GGA/)==xc_f03_func_info_get_family(xcpot%xc_info_x))
85 86
#else
    xcpot_is_gga=.false.
87 88 89
#endif
  END FUNCTION xcpot_is_gga

Matthias Redies's avatar
Matthias Redies committed
90 91 92 93 94
  LOGICAL FUNCTION xcpot_is_MetaGGA(xcpot)
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC    
    xcpot_is_MetaGGA=ANY((/XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA/)==xc_f03_func_info_get_family(xcpot%xc_info_x))
95 96
#else
    xcpot_is_MetaGGA=.false.
Matthias Redies's avatar
Matthias Redies committed
97 98 99
#endif
  END FUNCTION xcpot_is_MetaGGA

100 101 102 103
  LOGICAL FUNCTION xcpot_is_hybrid(xcpot)
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC
104
    xcpot_is_hybrid=ANY((/XC_FAMILY_HYB_MGGA, XC_FAMILY_HYB_GGA/)==xc_f03_func_info_get_family(xcpot%xc_info_x))
105 106
#else
    xcpot_is_hybrid=.false.
107 108 109 110 111 112 113 114 115 116
#endif
  END FUNCTION xcpot_is_hybrid

  FUNCTION xcpot_get_exchange_weight(xcpot) RESULT(a_ex)
    USE m_judft
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot

    REAL:: a_ex
#ifdef CPP_LIBXC    
117
    a_ex=xc_f03_hyb_exx_coef(xcpot%xc_func_x)
118 119 120 121
#endif
  END FUNCTION xcpot_get_exchange_weight

  !***********************************************************************
122
  SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad)
123 124 125 126
    !***********************************************************************
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
    INTEGER, INTENT (IN)     :: jspins
127 128 129
    REAL,INTENT (IN)         :: rh(:,:)   !Dimensions here
    REAL, INTENT (OUT)       :: vx (:,:)  !points,spin
    REAL, INTENT (OUT  )     :: vxc(:,:)  !
130
    ! optional arguments for GGA
131
    TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
132
#ifdef CPP_LIBXC    
133
    REAL,ALLOCATABLE::vxc_tmp(:,:),vx_tmp(:,:),vsigma(:,:)
134
    !libxc uses the spin as a first index, hence we have to transpose....
135 136
    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
137 138
    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
139 140 141
       ALLOCATE(vsigma,mold=grad%vsigma)
       !where(abs(grad%sigma)<1E-9) grad%sigma=1E-9
       CALL xc_f03_gga_vxc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vx_tmp,vsigma)
142
       IF (xcpot%func_id_c>0) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
143
          CALL xc_f03_gga_vxc(xcpot%xc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vxc_tmp,grad%vsigma)
144
          grad%vsigma=grad%vsigma+vsigma
145
          vxc_tmp=vxc_tmp+vx_tmp
Daniel Wortmann's avatar
Daniel Wortmann committed
146 147
       ELSE     
         vxc_tmp=vx_tmp
148
       ENDIF
149
    ELSE  !LDA potentials
150 151 152 153 154
       CALL xc_f03_lda_vxc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh), vx_tmp)
       IF (xcpot%func_id_c>0) THEN
          CALL xc_f03_lda_vxc(xcpot%xc_func_c, SIZE(rh,1), TRANSPOSE(rh), vxc_tmp)
          vxc_tmp=vxc_tmp+vx_tmp
       ENDIF
155
    ENDIF
156
    vx=TRANSPOSE(vx_tmp) 
157
    vxc=TRANSPOSE(vxc_tmp) 
158
    
159 160 161
#endif
  END SUBROUTINE xcpot_get_vxc

162
  
163 164 165 166 167 168
  !***********************************************************************
  SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad)
    !***********************************************************************
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
    INTEGER, INTENT (IN)     :: jspins
169 170
    REAL,INTENT (IN) :: rh(:,:)  !points,spin
    REAL, INTENT (OUT) :: exc(:) !points
171 172
    ! optional arguments for GGA      
    TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad
173 174

    REAL  :: excc(SIZE(exc))
175 176 177
#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")
178
       CALL xc_f03_gga_exc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,exc)
179 180 181 182
       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
183
    ELSE  !LDA potentials
184
       CALL xc_f03_lda_exc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh), exc)
185 186 187 188
       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
189
    ENDIF
190
    
191 192 193 194 195
#endif
  END SUBROUTINE xcpot_get_exc

  SUBROUTINE xcpot_alloc_gradients(ngrid,jspins,grad)
    INTEGER, INTENT (IN)         :: jspins,ngrid
196
    TYPE(t_gradients),INTENT(INOUT):: grad
197
    !For libxc we only need the sigma array...
198
    IF (ALLOCATED(grad%sigma)) DEALLOCATE(grad%sigma,grad%gr,grad%laplace,grad%vsigma)
199
    ALLOCATE(grad%sigma(MERGE(1,3,jspins==1),ngrid))
200
    ALLOCATE(grad%gr(3,ngrid,jspins))
201
    ALLOCATE(grad%laplace(ngrid,jspins))
202
    ALLOCATE(grad%vsigma(MERGE(1,3,jspins==1),ngrid))
203

204 205 206
  END SUBROUTINE xcpot_alloc_gradients


Daniel Wortmann's avatar
Daniel Wortmann committed
207
#ifdef CPP_LIBXC
208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248
  SUBROUTINE priv_write_info(xc_info)
    IMPLICIT NONE
    TYPE(xc_f03_func_info_t),INTENT(IN) :: xc_info
    INTEGER :: i
    CHARACTER(len=120) :: kind, family
    
    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
    
    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)
    
    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 priv_write_info
Daniel Wortmann's avatar
Daniel Wortmann committed
249
#endif
250
  
251
END MODULE m_types_xcpot_libxc