types_xcpot_libxc.F90 8.22 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 25 26 27 28 29 30 31 32 33 34 35
   CONTAINS
     PROCEDURE        :: is_gga=>xcpot_is_gga
     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

36
  SUBROUTINE xcpot_init(xcpot,jspins,id_x,id_c)
37 38 39
    USE m_judft
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(OUT)    :: xcpot
40 41
    INTEGER,INTENT(IN)                  :: jspins,id_x,id_c
    
42 43
#ifdef CPP_LIBXC
    INTEGER :: err
44 45 46 47
    xcpot%jspins=jspins
    xcpot%func_id_x=id_x
    xcpot%func_id_c=id_c
 
48
    IF (jspins==1) THEN
49 50
       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)
51
    ELSE
52 53
       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)
54
    END IF
55
    xcpot%xc_info_x=xc_f03_func_get_info(xcpot%xc_func_x)
56 57 58 59 60 61 62
    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
63 64 65 66 67 68 69 70 71 72 73

#else
    CALL judft_error("You specified a libxc-exchange correlation potential but FLEUR is not linked against libxc",hint="Please recompile FLEUR with libxc support")
#endif
  END SUBROUTINE xcpot_init


  LOGICAL FUNCTION xcpot_is_gga(xcpot)
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC    
74
    xcpot_is_gga=ANY((/XC_FAMILY_GGA, XC_FAMILY_HYB_GGA/)==xc_f03_func_info_get_family(xcpot%xc_info_x))
75 76 77 78 79 80 81
#endif
  END FUNCTION xcpot_is_gga

  LOGICAL FUNCTION xcpot_is_hybrid(xcpot)
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
#ifdef CPP_LIBXC
82
    xcpot_is_hybrid=ANY((/XC_FAMILY_HYB_MGGA, XC_FAMILY_HYB_GGA/)==xc_f03_func_info_get_family(xcpot%xc_info_x))
83 84 85 86 87 88 89 90 91 92
#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    
93
    a_ex=xc_f03_hyb_exx_coef(xcpot%xc_func_x)
94 95 96 97
#endif
  END FUNCTION xcpot_get_exchange_weight

  !***********************************************************************
98
  SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad)
99 100 101 102
    !***********************************************************************
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
    INTEGER, INTENT (IN)     :: jspins
103 104 105
    REAL,INTENT (IN)         :: rh(:,:)   !Dimensions here
    REAL, INTENT (OUT)       :: vx (:,:)  !points,spin
    REAL, INTENT (OUT  )     :: vxc(:,:)  !
106
    ! optional arguments for GGA
107
    TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
108
#ifdef CPP_LIBXC    
109
    REAL,ALLOCATABLE::vxc_tmp(:,:),vx_tmp(:,:),vsigma(:,:)
110
    !libxc uses the spin as a first index, hence we have to transpose....
111 112
    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
113 114
    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
115 116 117
       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)
118
       IF (xcpot%func_id_c>0) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
119
          CALL xc_f03_gga_vxc(xcpot%xc_func_c, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,vxc_tmp,grad%vsigma)
120
          grad%vsigma=grad%vsigma+vsigma
121
          vxc_tmp=vxc_tmp+vx_tmp
Daniel Wortmann's avatar
Daniel Wortmann committed
122 123
       ELSE     
         vxc_tmp=vx_tmp
124
       ENDIF
125
    ELSE  !LDA potentials
126 127 128 129 130
       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
131
    ENDIF
132
    vx=TRANSPOSE(vx_tmp) 
133
    vxc=TRANSPOSE(vxc_tmp) 
134
    
135 136 137
#endif
  END SUBROUTINE xcpot_get_vxc

138
  
139 140 141 142 143 144
  !***********************************************************************
  SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad)
    !***********************************************************************
    IMPLICIT NONE
    CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
    INTEGER, INTENT (IN)     :: jspins
145 146
    REAL,INTENT (IN) :: rh(:,:)  !points,spin
    REAL, INTENT (OUT) :: exc(:) !points
147 148
    ! optional arguments for GGA      
    TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad
149 150

    REAL  :: excc(SIZE(exc))
151 152 153
#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")
154
       CALL xc_f03_gga_exc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,exc)
155 156 157 158
       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
159
    ELSE  !LDA potentials
160
       CALL xc_f03_lda_exc(xcpot%xc_func_x, SIZE(rh,1), TRANSPOSE(rh), exc)
161 162 163 164
       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
165
    ENDIF
166
    
167 168 169 170 171
#endif
  END SUBROUTINE xcpot_get_exc

  SUBROUTINE xcpot_alloc_gradients(ngrid,jspins,grad)
    INTEGER, INTENT (IN)         :: jspins,ngrid
172
    TYPE(t_gradients),INTENT(INOUT):: grad
173
    !For libxc we only need the sigma array...
174
    IF (ALLOCATED(grad%sigma)) DEALLOCATE(grad%sigma,grad%gr,grad%laplace,grad%vsigma)
175
    ALLOCATE(grad%sigma(MERGE(1,3,jspins==1),ngrid))
176
    ALLOCATE(grad%gr(3,ngrid,jspins))
177 178 179
    ALLOCATE(grad%laplace(ngrid,jspins))
    ALLOCATE(grad%vsigma,mold=grad%sigma)

180 181 182
  END SUBROUTINE xcpot_alloc_gradients


Daniel Wortmann's avatar
Daniel Wortmann committed
183
#ifdef CPP_LIBXC
184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224
  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
225
#endif
226
  
227
END MODULE m_types_xcpot_libxc