types_xcpot_libxc.F90 11.7 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
   USE m_types_xcpot
   USE m_judft
   IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
15
16
17


#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
18
   PRIVATE :: write_xc_info
Matthias Redies's avatar
Matthias Redies committed
19
#endif
20

Matthias Redies's avatar
Matthias Redies committed
21
   TYPE,EXTENDS(t_xcpot):: t_xcpot_libxc 
22
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
23
24
      TYPE(xc_f03_func_t)      :: vxc_func_x, vxc_func_c
      TYPE(xc_f03_func_t)      :: exc_func_x, exc_func_c
25
#endif
Matthias Redies's avatar
Matthias Redies committed
26
27
28
      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
29
   CONTAINS
Matthias Redies's avatar
Matthias Redies committed
30
31
32
33
34
35
36
37
38
39
40
      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
41
42
CONTAINS

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

52
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
53
54
55
56
57
      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
58
      
59
60
      xcpot%func_exc_id_x = exc_id_x
      xcpot%func_exc_id_c = exc_id_c
Matthias Redies's avatar
Matthias Redies committed
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

      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
         ! 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

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

Matthias Redies's avatar
Matthias Redies committed
114
115
   LOGICAL FUNCTION xcpot_is_gga(xcpot)
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
116
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
Matthias Redies's avatar
Matthias Redies committed
117
#ifdef CPP_LIBXC
118
119
120
121
      TYPE(xc_f03_func_info_t)        :: xc_info

      xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
      xcpot_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
122
123
124
125
#else
      xcpot_is_gga=.false.
#endif
   END FUNCTION xcpot_is_gga
126

Matthias Redies's avatar
Matthias Redies committed
127
128
129
   LOGICAL FUNCTION xcpot_is_MetaGGA(xcpot)
      IMPLICIT NONE
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
Matthias Redies's avatar
Matthias Redies committed
130
#ifdef CPP_LIBXC    
Matthias Redies's avatar
Matthias Redies committed
131
      TYPE(xc_f03_func_info_t)        :: xc_info
Matthias Redies's avatar
Matthias Redies committed
132

Matthias Redies's avatar
Matthias Redies committed
133
134
      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))
135
#else
Matthias Redies's avatar
Matthias Redies committed
136
      xcpot_is_MetaGGA =  .False.
Matthias Redies's avatar
Matthias Redies committed
137
#endif
Matthias Redies's avatar
Matthias Redies committed
138
   END FUNCTION xcpot_is_MetaGGA
Matthias Redies's avatar
Matthias Redies committed
139

Matthias Redies's avatar
Matthias Redies committed
140
141
   LOGICAL FUNCTION xcpot_is_hybrid(xcpot)
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
142
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
143
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
144
      TYPE(xc_f03_func_info_t)        :: xc_info
Matthias Redies's avatar
Matthias Redies committed
145

Matthias Redies's avatar
Matthias Redies committed
146
147
      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))
148
#else
Matthias Redies's avatar
Matthias Redies committed
149
      xcpot_is_hybrid = .False.
150
#endif
Matthias Redies's avatar
Matthias Redies committed
151
   END FUNCTION xcpot_is_hybrid
152

Matthias Redies's avatar
Matthias Redies committed
153
154
155
   FUNCTION xcpot_get_exchange_weight(xcpot) RESULT(a_ex)
      USE m_judft
      IMPLICIT NONE
Matthias Redies's avatar
Matthias Redies committed
156
   CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
157

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

202
#endif
Matthias Redies's avatar
Matthias Redies committed
203
204
205
206
207
208
209
210
211
   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
212
      ! optional arguments for GGA
Matthias Redies's avatar
Matthias Redies committed
213
214
215
      TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad

      REAL  :: excc(SIZE(exc))
216
#ifdef CPP_LIBXC
Matthias Redies's avatar
Matthias Redies committed
217
218
219
      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)
220
         IF (xcpot%func_exc_id_c>0) THEN
Matthias Redies's avatar
Matthias Redies committed
221
222
223
224
225
            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)
226
         IF (xcpot%func_exc_id_c>0) THEN
Matthias Redies's avatar
Matthias Redies committed
227
228
229
230
231
            CALL xc_f03_lda_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh), excc)
            exc=exc+excc
         END IF
      ENDIF

Matthias Redies's avatar
Matthias Redies committed
232
233
#endif
   END SUBROUTINE xcpot_get_exc
234

Matthias Redies's avatar
Matthias Redies committed
235
236
237
238
239
240
241
242
243
   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))
244

Matthias Redies's avatar
Matthias Redies committed
245
   END SUBROUTINE xcpot_alloc_gradients
Matthias Redies's avatar
Matthias Redies committed
246

Matthias Redies's avatar
Matthias Redies committed
247
#ifdef CPP_LIBXC  
Matthias Redies's avatar
Matthias Redies committed
248
249
250
251
252
253
254
   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
255
256
257
      TYPE(xc_f03_func_t),INTENT(IN)      :: xc_func
      TYPE(xc_f03_func_info_t)            :: xc_info

Matthias Redies's avatar
Matthias Redies committed
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
      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
300
   END SUBROUTINE write_xc_info
Matthias Redies's avatar
Matthias Redies committed
301
#endif
Matthias Redies's avatar
Matthias Redies committed
302

303
END MODULE m_types_xcpot_libxc