Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : !>This module contains the xcpot-type providing an interface to libxc
8 : MODULE m_types_xcpot_libxc
9 : #ifdef CPP_LIBXC
10 : USE xc_f90_lib_m
11 : #endif
12 : USE m_types_xcpot
13 : USE m_judft
14 : use m_types_misc
15 : IMPLICIT NONE
16 :
17 : #ifdef CPP_LIBXC
18 : PRIVATE :: write_xc_info
19 : #endif
20 :
21 : TYPE,EXTENDS(t_xcpot):: t_xcpot_libxc
22 : #ifdef CPP_LIBXC
23 : TYPE(xc_f90_func_t) :: vxc_func_x, vxc_func_c
24 : TYPE(xc_f90_func_t) :: exc_func_x, exc_func_c
25 : #endif
26 : INTEGER :: jspins
27 :
28 : CONTAINS
29 : !PROCEDURE :: vxc_is_LDA => xcpot_vxc_is_LDA
30 : !PROCEDURE :: vxc_is_gga => xcpot_vxc_is_gga
31 :
32 : PROCEDURE :: vx_is_LDA => xcpot_vx_is_LDA
33 : PROCEDURE :: vx_is_GGA => xcpot_vx_is_GGA
34 : PROCEDURE :: vx_is_MetaGGA => xcpot_vx_is_MetaGGA
35 :
36 : PROCEDURE :: vc_is_LDA => xcpot_vc_is_LDA
37 : PROCEDURE :: vc_is_GGA => xcpot_vc_is_GGA
38 :
39 : PROCEDURE :: exc_is_LDA => xcpot_exc_is_LDA
40 : PROCEDURE :: exc_is_gga => xcpot_exc_is_gga
41 : PROCEDURE :: exc_is_MetaGGA => xcpot_exc_is_MetaGGA
42 :
43 : PROCEDURE :: is_hybrid => xcpot_is_hybrid
44 : PROCEDURE :: get_exchange_weight => xcpot_get_exchange_weight
45 : PROCEDURE :: get_vxc => xcpot_get_vxc
46 : PROCEDURE :: get_exc => xcpot_get_exc
47 : PROCEDURE :: get_fxc => xcpot_get_fxc
48 : PROCEDURE, NOPASS :: alloc_gradients => xcpot_alloc_gradients
49 : !Not overloeaded...
50 : PROCEDURE :: init => xcpot_init
51 : PROCEDURE,NOPASS :: apply_cutoffs
52 : END TYPE t_xcpot_libxc
53 : PUBLIC t_xcpot_libxc
54 : CONTAINS
55 5 : subroutine apply_cutoffs(density_cutoff,rh,grad)
56 : real,intent(INOUT) :: rh(:,:)
57 : real,INTENT(IN) :: density_cutoff
58 : type(t_gradients),INTENT(INOUT),OPTIONAL :: grad
59 :
60 :
61 :
62 : integer:: i,j
63 13 : DO j=1,size(rh,2)
64 99565 : DO i=1,size(rh,1)
65 99560 : if (abs(rh(i,j))<density_cutoff) THEN
66 0 : rh(i,j)=density_cutoff
67 0 : if (present(grad)) Then
68 0 : if (allocated(grad%sigma)) grad%sigma(:,i)=0.0 !if one spin is small, apply cutoff to all gradients!
69 0 : if (allocated(grad%gr)) grad%gr(:,i,j)=0.0
70 0 : if (allocated(grad%laplace)) grad%laplace(i,j)=0.0
71 : endif
72 : endif
73 : ENDDO
74 : ENDDO
75 :
76 5 : end subroutine
77 :
78 6 : SUBROUTINE xcpot_init(xcpot, func_vxc_id_x, func_vxc_id_c, func_exc_id_x, func_exc_id_c, jspins)
79 : USE m_judft
80 : IMPLICIT NONE
81 : CLASS(t_xcpot_libxc), INTENT(INOUT) :: xcpot
82 : INTEGER, INTENT(IN) :: jspins, func_vxc_id_x, func_vxc_id_c, func_exc_id_x, func_exc_id_c
83 : LOGICAL :: same_functionals ! are vxc and exc equal
84 : INTEGER :: errors(4)
85 :
86 : #ifdef CPP_LIBXC
87 30 : errors = -1
88 6 : xcpot%jspins = jspins
89 6 : xcpot%func_vxc_id_x = func_vxc_id_x
90 6 : xcpot%func_exc_id_x = func_exc_id_x
91 6 : xcpot%func_vxc_id_c = func_vxc_id_c
92 6 : xcpot%func_exc_id_c = func_exc_id_c
93 :
94 6 : IF (xcpot%func_vxc_id_x == 0 .OR. xcpot%func_exc_id_x == 0) THEN
95 : CALL judft_error("LibXC exchange- and correlation-function indicies need to be set" &
96 : , hint='Try this: '//ACHAR(10)// &
97 : '<xcFunctional name="libxc" relativisticCorrections="F">'//ACHAR(10)// &
98 : ' <libXC exchange="1" correlation="1" /> '//ACHAR(10)// &
99 0 : '</xcFunctional> ')
100 : ENDIF
101 :
102 6 : IF (jspins==1) THEN
103 : ! potential functionals
104 4 : CALL xc_f90_func_init(xcpot%vxc_func_x, xcpot%func_vxc_id_x, XC_UNPOLARIZED, err=errors(1))
105 4 : IF (xcpot%func_vxc_id_c>0) CALL xc_f90_func_init(xcpot%vxc_func_c, xcpot%func_vxc_id_c, &
106 4 : XC_UNPOLARIZED, err=errors(2))
107 :
108 : ! energy functionals
109 4 : CALL xc_f90_func_init(xcpot%exc_func_x, xcpot%func_exc_id_x, XC_UNPOLARIZED, err=errors(3))
110 4 : IF (xcpot%func_exc_id_c>0) CALL xc_f90_func_init(xcpot%exc_func_c, xcpot%func_exc_id_c, &
111 4 : XC_UNPOLARIZED, err=errors(4))
112 :
113 : ELSE
114 : ! potential functionals
115 2 : CALL xc_f90_func_init(xcpot%vxc_func_x, xcpot%func_vxc_id_x, XC_POLARIZED, err=errors(1))
116 2 : IF (xcpot%func_vxc_id_c>0) CALL xc_f90_func_init(xcpot%vxc_func_c, xcpot%func_vxc_id_c, &
117 2 : XC_POLARIZED, err=errors(2))
118 :
119 : !energy functionals
120 2 : CALL xc_f90_func_init(xcpot%exc_func_x, xcpot%func_exc_id_x, XC_POLARIZED, err=errors(3))
121 2 : IF (xcpot%func_exc_id_c>0) CALL xc_f90_func_init(xcpot%exc_func_c, xcpot%func_exc_id_c, &
122 2 : XC_POLARIZED, err=errors(4))
123 : END IF
124 :
125 : !IF(errors(1) /= 0) call juDFT_error("Exchange potential functional not in LibXC")
126 : !IF(errors(2) /= 0) call juDFT_error("Correlation potential functional not in LibXC")
127 : !IF(errors(3) /= 0) call juDFT_error("Exchange energy functional not in LibXC")
128 : !IF(errors(4) /= 0) call juDFT_error("Correlation energy functional not in LibXC")
129 :
130 : !check if any potental is a MetaGGA
131 12 : IF (ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA] == xc_get_family(xcpot%vxc_func_x))) THEN
132 0 : CALL juDFT_error("vxc_x: MetaGGA is not implemented for potentials")
133 6 : ELSEIF (xcpot%func_vxc_id_c > 0) THEN
134 12 : IF (ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA] == xc_get_family(xcpot%vxc_func_c))) THEN
135 0 : CALL juDFT_error("vxc_x: MetaGGA is not implemented for potentials")
136 : ENDIF
137 : ENDIF
138 :
139 6 : CALL write_xc_info(xcpot%vxc_func_x)
140 :
141 6 : IF (xcpot%func_vxc_id_c > 0) THEN
142 6 : CALL write_xc_info(xcpot%vxc_func_c)
143 : ELSE
144 0 : WRITE (*, *) "No Correlation functional"
145 : END IF
146 :
147 : same_functionals = (xcpot%func_vxc_id_x == xcpot%func_exc_id_x) &
148 6 : .AND. (xcpot%func_vxc_id_c == xcpot%func_exc_id_c)
149 : IF (.NOT. same_functionals) THEN
150 0 : CALL write_xc_info(xcpot%exc_func_x)
151 0 : IF (xcpot%func_exc_id_c > 0) THEN
152 0 : CALL write_xc_info(xcpot%exc_func_c)
153 : ELSE
154 0 : WRITE (*, *) "No Correlation functional for TotalE"
155 : ENDIF
156 : ELSE
157 6 : WRITE (*, *) "Using same functional for VXC and EXC"
158 : END IF
159 : #else
160 : CALL judft_error("You specified a libxc-exchange correlation potential but FLEUR is not linked against libxc", &
161 : hint="Please recompile FLEUR with libxc support")
162 : #endif
163 6 : END SUBROUTINE xcpot_init
164 :
165 : ! LDA
166 0 : LOGICAL FUNCTION xcpot_vx_is_LDA(xcpot)
167 : IMPLICIT NONE
168 : CLASS(t_xcpot_libxc), INTENT(IN):: xcpot
169 : #ifdef CPP_LIBXC
170 : TYPE(xc_f90_func_info_t) :: xc_info
171 :
172 0 : xc_info = xc_f90_func_get_info(xcpot%vxc_func_x)
173 0 : xcpot_vx_is_LDA = XC_FAMILY_LDA == xc_f90_func_info_get_family(xc_info)
174 : #else
175 : xcpot_vx_is_LDA = .false.
176 : #endif
177 0 : END FUNCTION xcpot_vx_is_LDA
178 :
179 0 : LOGICAL FUNCTION xcpot_vc_is_LDA(xcpot)
180 : IMPLICIT NONE
181 : CLASS(t_xcpot_libxc), INTENT(IN):: xcpot
182 : #ifdef CPP_LIBXC
183 : TYPE(xc_f90_func_info_t) :: xc_info
184 :
185 0 : xc_info = xc_f90_func_get_info(xcpot%vxc_func_c)
186 0 : xcpot_vc_is_LDA = XC_FAMILY_LDA == xc_f90_func_info_get_family(xc_info)
187 : #else
188 : xcpot_vc_is_LDA = .false.
189 : #endif
190 0 : END FUNCTION xcpot_vc_is_LDA
191 :
192 3 : LOGICAL FUNCTION xcpot_exc_is_LDA(xcpot)
193 : IMPLICIT NONE
194 : CLASS(t_xcpot_libxc), INTENT(IN):: xcpot
195 : #ifdef CPP_LIBXC
196 : TYPE(xc_f90_func_info_t) :: xc_info
197 :
198 3 : xc_info = xc_f90_func_get_info(xcpot%exc_func_x)
199 3 : xcpot_exc_is_LDA = (XC_FAMILY_LDA == xc_f90_func_info_get_family(xc_info))
200 : #else
201 : xcpot_exc_is_LDA = .false.
202 : #endif
203 3 : END FUNCTION xcpot_exc_is_LDA
204 :
205 : ! GGA
206 191 : LOGICAL FUNCTION xcpot_vc_is_gga(xcpot)
207 : IMPLICIT NONE
208 : CLASS(t_xcpot_libxc), INTENT(IN):: xcpot
209 : #ifdef CPP_LIBXC
210 : TYPE(xc_f90_func_info_t) :: xc_info
211 :
212 191 : xc_info = xc_f90_func_get_info(xcpot%vxc_func_c)
213 191 : xcpot_vc_is_gga = ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f90_func_info_get_family(xc_info))
214 : #else
215 : xcpot_vc_is_gga = .false.
216 : #endif
217 191 : END FUNCTION xcpot_vc_is_gga
218 :
219 0 : LOGICAL FUNCTION xcpot_vx_is_gga(xcpot)
220 : IMPLICIT NONE
221 : CLASS(t_xcpot_libxc), INTENT(IN):: xcpot
222 : #ifdef CPP_LIBXC
223 : TYPE(xc_f90_func_info_t) :: xc_info
224 :
225 0 : xc_info = xc_f90_func_get_info(xcpot%vxc_func_x)
226 0 : xcpot_vx_is_gga = ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f90_func_info_get_family(xc_info))
227 : #else
228 : xcpot_vx_is_gga = .false.
229 : #endif
230 0 : END FUNCTION xcpot_vx_is_gga
231 :
232 26 : LOGICAL FUNCTION xcpot_vx_is_MetaGGA(xcpot)
233 : IMPLICIT NONE
234 : CLASS(t_xcpot_libxc), INTENT(IN):: xcpot
235 : #ifdef CPP_LIBXC
236 : TYPE(xc_f90_func_info_t) :: xc_info
237 :
238 26 : xc_info = xc_f90_func_get_info(xcpot%vxc_func_x)
239 26 : xcpot_vx_is_MetaGGA = ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA]==xc_f90_func_info_get_family(xc_info))
240 : #else
241 : xcpot_vx_is_MetaGGA = .false.
242 : #endif
243 26 : END FUNCTION xcpot_vx_is_MetaGGA
244 :
245 14 : LOGICAL FUNCTION xcpot_exc_is_gga(xcpot)
246 : IMPLICIT NONE
247 : CLASS(t_xcpot_libxc), INTENT(IN):: xcpot
248 : #ifdef CPP_LIBXC
249 : TYPE(xc_f90_func_info_t) :: xc_info
250 :
251 14 : xc_info = xc_f90_func_get_info(xcpot%exc_func_x)
252 14 : xcpot_exc_is_gga = ANY([XC_FAMILY_GGA, XC_FAMILY_HYB_GGA]==xc_f90_func_info_get_family(xc_info))
253 : #else
254 : xcpot_exc_is_gga = .false.
255 : #endif
256 14 : END FUNCTION xcpot_exc_is_gga
257 :
258 40 : LOGICAL FUNCTION xcpot_exc_is_MetaGGA(xcpot)
259 : IMPLICIT NONE
260 : CLASS(t_xcpot_libxc),INTENT(IN):: xcpot
261 : #ifdef CPP_LIBXC
262 : TYPE(xc_f90_func_info_t) :: xc_info
263 :
264 40 : xc_info = xc_f90_func_get_info(xcpot%exc_func_x)
265 40 : xcpot_exc_is_MetaGGA=ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA]==xc_f90_func_info_get_family(xc_info))
266 : #else
267 : xcpot_exc_is_MetaGGA = .False.
268 : #endif
269 40 : END FUNCTION xcpot_exc_is_MetaGGA
270 :
271 21 : LOGICAL FUNCTION xcpot_is_hybrid(xcpot)
272 : IMPLICIT NONE
273 : CLASS(t_xcpot_libxc), INTENT(IN):: xcpot
274 : #ifdef CPP_LIBXC
275 : TYPE(xc_f90_func_info_t) :: xc_info
276 :
277 21 : xc_info = xc_f90_func_get_info(xcpot%vxc_func_x)
278 21 : xcpot_is_hybrid=ANY([XC_FAMILY_HYB_MGGA, XC_FAMILY_HYB_GGA]==xc_f90_func_info_get_family(xc_info))
279 : #else
280 : xcpot_is_hybrid = .False.
281 : #endif
282 21 : END FUNCTION xcpot_is_hybrid
283 :
284 26 : FUNCTION xcpot_get_exchange_weight(xcpot) RESULT(a_ex)
285 : USE m_judft
286 : IMPLICIT NONE
287 : CLASS(t_xcpot_libxc), INTENT(IN):: xcpot
288 :
289 : REAL:: a_ex
290 : #ifdef CPP_LIBXC
291 26 : a_ex=xc_f90_hyb_exx_coef(xcpot%vxc_func_x)
292 : #endif
293 26 : END FUNCTION xcpot_get_exchange_weight
294 :
295 : !***********************************************************************
296 61 : SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad, kinenergyden_ks)
297 : USE, INTRINSIC :: IEEE_ARITHMETIC
298 : use iso_c_binding
299 : IMPLICIT NONE
300 : CLASS(t_xcpot_libxc), INTENT(IN) :: xcpot
301 : INTEGER, INTENT(IN) :: jspins
302 : REAL, INTENT(IN) :: rh(:, :) !Dimensions here
303 : REAL, INTENT(OUT) :: vx(:, :) !points,spin
304 : REAL, INTENT(OUT) :: vxc(:, :) !
305 : ! optional arguments for GGA
306 : TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
307 : REAL, INTENT(IN), OPTIONAL :: kinenergyden_ks(:, :)
308 : #ifdef CPP_LIBXC
309 61 : REAL,ALLOCATABLE :: vxc_tmp(:,:),vx_tmp(:,:),vsigma(:,:), &
310 : tmp_vsig(:,:), tmp_vlapl(:,:), tmp_vtau(:,:), &
311 : kinED_libxc(:,:)
312 : integer(kind=c_size_t) :: idx
313 : !libxc uses the spin as a first index, hence we have to transpose....
314 3903896 : ALLOCATE (vxc_tmp(SIZE(vxc, 2), SIZE(vxc, 1))); vxc_tmp = 0.0
315 3903896 : ALLOCATE (vx_tmp(SIZE(vx, 2), SIZE(vx, 1))); vx_tmp = 0.0
316 :
317 61 : IF (xcpot%needs_grad()) THEN
318 36 : IF (.NOT. PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives")
319 144 : ALLOCATE (vsigma, mold=grad%vsigma)
320 : !where(abs(grad%sigma)<1E-9) grad%sigma=1E-9
321 36 : CALL xc_f90_gga_vxc(xcpot%vxc_func_x, SIZE(rh, 1, kind=c_size_t), TRANSPOSE(rh), grad%sigma, vx_tmp, vsigma)
322 36 : IF (xcpot%func_vxc_id_c > 0) THEN
323 36 : CALL xc_f90_gga_vxc(xcpot%vxc_func_c, SIZE(rh, 1, kind=c_size_t), TRANSPOSE(rh), grad%sigma, vxc_tmp, grad%vsigma)
324 3026060 : grad%vsigma = grad%vsigma + vsigma
325 2421275 : vxc_tmp = vxc_tmp + vx_tmp
326 : ELSE
327 0 : vxc_tmp = vx_tmp
328 : ENDIF
329 : ELSE !LDA potentials
330 25 : CALL xc_f90_lda_vxc(xcpot%vxc_func_x, SIZE(rh, 1, kind=c_size_t), TRANSPOSE(rh), vx_tmp)
331 25 : IF (xcpot%func_vxc_id_c > 0) THEN
332 25 : CALL xc_f90_lda_vxc(xcpot%vxc_func_c, SIZE(rh, 1, kind=c_size_t), TRANSPOSE(rh), vxc_tmp)
333 1482438 : vxc_tmp = vxc_tmp + vx_tmp
334 : ENDIF
335 : ENDIF
336 2254394 : vx = TRANSPOSE(vx_tmp)
337 2254394 : vxc = TRANSPOSE(vxc_tmp)
338 :
339 : #endif
340 61 : END SUBROUTINE xcpot_get_vxc
341 :
342 14 : SUBROUTINE xcpot_get_exc(xcpot, jspins, rh, exc, grad, kinEnergyDen_KS, mt_call)
343 61 : use m_constants
344 : use ISO_C_BINDING
345 : IMPLICIT NONE
346 : CLASS(t_xcpot_libxc), INTENT(IN) :: xcpot
347 : INTEGER, INTENT(IN) :: jspins
348 : REAL, INTENT(IN) :: rh(:, :) !points,spin
349 : REAL, INTENT(OUT) :: exc(:) !points
350 : ! optional arguments for GGA
351 : TYPE(t_gradients), OPTIONAL, INTENT(IN) :: grad
352 : LOGICAL, OPTIONAL, INTENT(IN) :: mt_call
353 :
354 : ! kinED from Kohn-Sham equations:
355 : ! tau = sum[phi_i(r)^dag nabla phi_i(r)]
356 : ! see eq (2) in https://doi.org/10.1063/1.1565316
357 : ! (-0.5 is applied below)
358 : REAL, INTENT(IN), OPTIONAL :: kinEnergyDen_KS(:, :)
359 :
360 : #ifdef CPP_LIBXC
361 : TYPE(xc_f90_func_info_t) :: xc_info
362 14 : REAL :: excc(SIZE(exc))
363 : REAL :: cut_ratio = 0.1
364 : INTEGER :: cut_idx
365 : LOGICAL :: is_mt
366 :
367 : ! tau = 0.5 * sum[|grad phi_i(r)|²]
368 : ! see eq (3) in https://doi.org/10.1063/1.1565316
369 14 : REAL, ALLOCATABLE :: kinEnergyDen_libXC(:, :), pkzb_ratio(:, :), pkzb_zaehler(:, :), pkzb_nenner(:, :)
370 :
371 14 : is_mt = merge(mt_call, .False., present(mt_call))
372 14 : IF (xcpot%exc_is_gga()) THEN
373 11 : IF (.NOT. PRESENT(grad)) CALL judft_error("Bug: You called get_exc for a GGA potential without providing derivatives")
374 11 : CALL xc_f90_gga_exc(xcpot%exc_func_x, SIZE(rh, 1, kind=c_size_t), TRANSPOSE(rh), grad%sigma, exc)
375 11 : IF (xcpot%func_exc_id_c > 0) THEN
376 11 : CALL xc_f90_gga_exc(xcpot%exc_func_c, SIZE(rh, 1, kind=c_size_t), TRANSPOSE(rh), grad%sigma, excc)
377 891345 : exc = exc + excc
378 : END IF
379 3 : ELSEIF (xcpot%exc_is_LDA()) THEN !LDA potentials
380 3 : CALL xc_f90_lda_exc(xcpot%exc_func_x, SIZE(rh, 1, kind=c_size_t), TRANSPOSE(rh), exc)
381 3 : IF (xcpot%func_exc_id_c > 0) THEN
382 3 : CALL xc_f90_lda_exc(xcpot%exc_func_c, SIZE(rh, 1, kind=c_size_t), TRANSPOSE(rh), excc)
383 716711 : exc = exc + excc
384 : END IF
385 0 : ELSEIF (xcpot%exc_is_MetaGGA()) THEN
386 0 : IF (PRESENT(kinEnergyDen_KS)) THEN
387 : ! apply correction in eq (4) in https://doi.org/10.1063/1.1565316
388 0 : kinEnergyDen_libXC = transpose(kinEnergyDen_KS + 0.25*grad%laplace)
389 :
390 : !only cut core of muffin tin
391 0 : cut_idx = MERGE(NINT(size(rh, 1)*cut_ratio), 0, is_mt)
392 :
393 0 : exc = 0.0
394 0 : excc = 0.0
395 : call xc_f90_mgga_exc(xcpot%exc_func_x, SIZE(rh(cut_idx + 1:, :), 1, kind=c_size_t), &
396 : TRANSPOSE(rh(cut_idx + 1:, :)), &
397 : grad%sigma(:, cut_idx + 1:), &
398 : transpose(grad%laplace(cut_idx + 1:, :)), &
399 : kinEnergyDen_libXC(:, cut_idx + 1:), &
400 0 : exc(cut_idx + 1:))
401 :
402 : call xc_f90_gga_exc(xcpot%vxc_func_x, SIZE(rh(:cut_idx, :), 1, kind=c_size_t), &
403 : TRANSPOSE(rh(:cut_idx, :)), &
404 : grad%sigma(:, :cut_idx), &
405 0 : exc(:cut_idx))
406 :
407 0 : IF (xcpot%func_exc_id_c > 0) THEN
408 : call xc_f90_mgga_exc(xcpot%exc_func_c, SIZE(rh(cut_idx + 1:, :), 1, kind=c_size_t), &
409 : TRANSPOSE(rh(cut_idx + 1:, :)), &
410 : grad%sigma(:, cut_idx + 1:), &
411 : transpose(grad%laplace(cut_idx + 1:, :)), &
412 : kinEnergyDen_libXC(:, cut_idx + 1:), &
413 0 : excc(cut_idx + 1:))
414 :
415 : call xc_f90_gga_exc(xcpot%vxc_func_c, SIZE(rh(:cut_idx, :), 1, kind=c_size_t), &
416 : TRANSPOSE(rh(:cut_idx, :)), &
417 : grad%sigma(:, :cut_idx), &
418 0 : excc(:cut_idx))
419 0 : exc = exc + excc
420 : END IF
421 :
422 : ELSE ! first iteration is GGA
423 0 : IF (.NOT. PRESENT(grad)) CALL judft_error("Bug: You called get_exc for a MetaGGA potential without providing derivatives")
424 0 : CALL xc_f90_gga_exc(xcpot%vxc_func_x, SIZE(rh, 1, kind=c_size_t), TRANSPOSE(rh), grad%sigma, exc)
425 0 : IF (xcpot%func_exc_id_c > 0) THEN
426 0 : CALL xc_f90_gga_exc(xcpot%vxc_func_c, SIZE(rh, 1, kind=c_size_t), TRANSPOSE(rh), grad%sigma, excc)
427 0 : exc = exc + excc
428 : END IF
429 : ENDIF
430 :
431 : ELSE
432 0 : call juDFT_error("exc is part of a known Family", calledby="xcpot_get_exc@libxc")
433 : ENDIF
434 :
435 : #endif
436 14 : END SUBROUTINE xcpot_get_exc
437 :
438 0 : SUBROUTINE xcpot_get_fxc(xcpot, jspins, rh, fxc)
439 : USE, INTRINSIC :: IEEE_ARITHMETIC
440 : use iso_c_binding
441 :
442 : IMPLICIT NONE
443 :
444 : CLASS(t_xcpot_libxc), INTENT(IN) :: xcpot
445 : INTEGER, INTENT(IN) :: jspins
446 : REAL, INTENT(IN) :: rh(:, :)
447 : REAL, INTENT(OUT) :: fxc(:, :)
448 :
449 : #ifdef CPP_LIBXC
450 0 : REAL,ALLOCATABLE :: fxc_tmp(:,:), fx_tmp(:,:)
451 :
452 : integer(kind=c_size_t) :: idx
453 :
454 : !libxc uses the spin as a first index, hence we have to transpose....
455 0 : ALLOCATE (fxc_tmp(SIZE(fxc, 2), SIZE(fxc, 1))); fxc_tmp = 0.0
456 0 : ALLOCATE (fx_tmp(SIZE(fxc, 2), SIZE(fxc, 1))); fx_tmp = 0.0
457 :
458 0 : IF (xcpot%needs_grad().OR.xcpot%exc_is_MetaGGA()) THEN
459 0 : CALL judft_error("Bug: You called get_fxc for a (meta)GGA potential. This is not implemented (yet?).")
460 : ELSE !LDA potentials
461 0 : CALL xc_f90_lda_fxc(xcpot%vxc_func_x, SIZE(rh, 1, kind=c_size_t), TRANSPOSE(rh), fx_tmp)
462 0 : IF (xcpot%func_vxc_id_c > 0) THEN
463 0 : CALL xc_f90_lda_fxc(xcpot%vxc_func_c, SIZE(rh, 1, kind=c_size_t), TRANSPOSE(rh), fxc_tmp)
464 0 : fxc_tmp = fxc_tmp + fx_tmp
465 : ENDIF
466 : ENDIF
467 0 : fxc = TRANSPOSE(fxc_tmp)
468 :
469 : #endif
470 0 : END SUBROUTINE xcpot_get_fxc
471 :
472 44 : SUBROUTINE xcpot_alloc_gradients(ngrid, jspins, grad)
473 : INTEGER, INTENT(IN) :: jspins, ngrid
474 : TYPE(t_gradients), INTENT(INOUT):: grad
475 : !For libxc we only need the sigma array...
476 44 : IF (ALLOCATED(grad%sigma)) DEALLOCATE (grad%sigma, grad%gr, grad%laplace, grad%vsigma)
477 216 : ALLOCATE (grad%sigma(MERGE(1, 3, jspins == 1), ngrid))
478 176 : ALLOCATE (grad%gr(3, ngrid, jspins))
479 176 : ALLOCATE (grad%laplace(ngrid, jspins))
480 88 : ALLOCATE (grad%vsigma(MERGE(1, 3, jspins == 1), ngrid))
481 :
482 0 : END SUBROUTINE xcpot_alloc_gradients
483 :
484 0 : subroutine mpi_bc_xcpot_libxc(This, Mpi_comm, Irank)
485 : Use M_mpi_bc_tool
486 : Class(t_xcpot_libxc), Intent(Inout)::This
487 : Integer, Intent(In):: Mpi_comm
488 : Integer, Intent(In), Optional::Irank
489 : Integer ::Rank
490 0 : If (Present(Irank)) Then
491 0 : Rank = Irank
492 : Else
493 0 : Rank = 0
494 : End If
495 :
496 : ! Bcasts for abstract base class t_xcpot
497 0 : CALL mpi_bc(this%l_libxc, rank, mpi_comm)
498 0 : CALL mpi_bc(this%func_vxc_id_c, rank, mpi_comm)
499 0 : CALL mpi_bc(this%func_vxc_id_x, rank, mpi_comm)
500 0 : CALL mpi_bc(this%func_exc_id_c, rank, mpi_comm)
501 0 : CALL mpi_bc(this%func_exc_id_x, rank, mpi_comm)
502 0 : CALL mpi_bc(this%l_inbuild, rank, mpi_comm)
503 0 : CALL mpi_bc(rank, mpi_comm, this%inbuild_name)
504 0 : CALL mpi_bc(this%l_relativistic, rank, mpi_comm)
505 :
506 0 : END SUBROUTINE mpi_bc_xcpot_libxc
507 : #ifdef CPP_LIBXC
508 12 : SUBROUTINE write_xc_info(xc_func, is_E_func)
509 : IMPLICIT NONE
510 : LOGICAL, INTENT(IN), OPTIONAL :: is_E_func
511 : INTEGER :: i
512 : CHARACTER(len=120) :: kind, family
513 : LOGICAL :: is_energy_func
514 :
515 : TYPE(xc_f90_func_t),INTENT(IN) :: xc_func
516 : TYPE(xc_f90_func_info_t) :: xc_info
517 :
518 12 : xc_info = xc_f90_func_get_info(xc_func)
519 12 : is_energy_func = merge(is_E_func, .False., PRESENT(is_E_func))
520 :
521 6 : SELECT CASE(xc_f90_func_info_get_kind(xc_info))
522 : CASE (XC_EXCHANGE)
523 6 : WRITE (kind, '(a)') 'an exchange functional'
524 : CASE (XC_CORRELATION)
525 6 : WRITE (kind, '(a)') 'a correlation functional'
526 : CASE (XC_EXCHANGE_CORRELATION)
527 0 : WRITE (kind, '(a)') 'an exchange-correlation functional'
528 : CASE (XC_KINETIC)
529 0 : WRITE (kind, '(a)') 'a kinetic energy functional'
530 : CASE default
531 12 : WRITE (kind, '(a)') 'of unknown kind'
532 : END SELECT
533 4 : SELECT CASE (xc_f90_func_info_get_family(xc_info))
534 : CASE (XC_FAMILY_LDA);
535 4 : WRITE (family, '(a)') "LDA"
536 : CASE (XC_FAMILY_GGA);
537 8 : WRITE (family, '(a)') "GGA"
538 : CASE (XC_FAMILY_HYB_GGA);
539 0 : WRITE (family, '(a)') "hybrid GGA"
540 : CASE (XC_FAMILY_MGGA);
541 0 : WRITE (family, '(a)') "MGGA"
542 : CASE (XC_FAMILY_HYB_MGGA);
543 0 : WRITE (family, '(a)') "hybrid MGGA"
544 : CASE default;
545 12 : WRITE (family, '(a)') "unknown"
546 : END SELECT
547 :
548 12 : IF(.not. is_energy_func) THEN
549 : WRITE(*,'("The functional ''", a, "'' is ", a, ", it belongs to the ''", a, "'' family and is defined in the reference(s):")') &
550 12 : TRIM(xc_f90_func_info_get_name(xc_info)), TRIM(kind), TRIM(family)
551 : ELSE
552 : WRITE(*,'("The functional used for TotalE ''", a, "'' is ", a, ", it belongs to the ''", a, "'' family and is defined in the reference(s):")') &
553 0 : TRIM(xc_f90_func_info_get_name(xc_info)), TRIM(kind), TRIM(family)
554 : ENDIF
555 :
556 12 : i = 0
557 34 : DO WHILE(i >= 0)
558 34 : WRITE(*, '(a,i1,2a)') '[', i+1, '] ', TRIM(xc_f90_func_reference_get_ref(xc_f90_func_info_get_references(xc_info, i)))
559 : END DO
560 12 : END SUBROUTINE write_xc_info
561 :
562 12 : FUNCTION xc_get_family(xc_func) result(family)
563 : IMPLICIT NONE
564 : TYPE(xc_f90_func_t) :: xc_func
565 : integer :: family
566 12 : family = xc_f90_func_info_get_family(xc_f90_func_get_info(xc_func))
567 0 : END FUNCTION xc_get_family
568 : #endif
569 :
570 6 : END MODULE m_types_xcpot_libxc
|
