Line data Source code
1 : !!****m* ABINIT/m_drivexc
2 : !! NAME
3 : !! m_drivexc
4 : !!
5 : !! FUNCTION
6 : !! Driver of XC functionals. Optionally, deliver the XC kernel, or even the derivative
7 : !! of the XC kernel (the third derivative of the XC energy)
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2012-2026 ABINIT group (MT, MJV, CE, TD, XG, AB)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_drivexc
24 :
25 : use defs_basis
26 : use m_abicore
27 : use m_errors
28 : use libxc_functionals
29 : use m_numeric_tools, only: invcb
30 : use m_xciit, only: xciit
31 : use m_xcpbe, only: xcpbe,xckdt16
32 : use m_xchcth, only: xchcth
33 : use m_xclda, only: xcpzca,xcspol,xctetr,xcwign,xchelu,xcxalp,xclb,xcksdt
34 :
35 : implicit none
36 :
37 : private
38 : !!***
39 :
40 : public :: drivexc ! Driver of XC functionals. Optionally, deliver the XC kernel, or even the derivative
41 : public :: echo_xc_name ! Write to log and output the xc functional which will be used for this dataset
42 : public :: xc_need_kden ! Given a XC functional (defined by ixc), return TRUE if it needs kinetic energy density.
43 : public :: has_kxc ! Given a XC functional (defined by ixc), return TRUE if Kxc (dVxc/drho) is avalaible.
44 : public :: has_k3xc ! Given a XC functional (defined by ixc), return TRUE if K3xc (d2Vxc/drho2) is avalaible.
45 : public :: check_kxc ! Given a XC functional (defined by ixc), check if Kxc and/or K3xc is avalaible.
46 : public :: size_dvxc ! Give the size of the array dvxc(npts,ndvxc) and the second dimension of the d2vxc(npts,nd2vxc)
47 : public :: xcmult ! (GGA) Multiply the different gradient of spin-density by the derivative of the XC functional
48 : ! with respect to the norm of the gradient, then divide it by the norm of the gradient
49 : public :: mkdenpos ! Make a ground-state density positive everywhere.
50 : !!***
51 :
52 : contains
53 : !!***
54 :
55 : !!****f* m_drivexc/echo_xc_name
56 : !! NAME
57 : !! echo_xc_name
58 : !!
59 : !! FUNCTION
60 : !! Write to log and output the xc functional which will be used for this dataset
61 : !!
62 : !! INPUTS
63 : !! ixc = internal code for xc functional
64 : !!
65 : !! SOURCE
66 :
67 5281 : subroutine echo_xc_name (ixc)
68 :
69 : !Arguments -------------------------------
70 : integer, intent(in) :: ixc
71 :
72 : !Local variables -------------------------
73 : integer :: l_citation
74 : character(len=500) :: message, citation
75 : ! *********************************************************************
76 :
77 5281 : message =''
78 5281 : citation =''
79 :
80 : !normal case (not libxc)
81 5281 : if (ixc >= 0) then
82 :
83 39 : select case (ixc)
84 : case (0)
85 39 : message = 'No xc applied (usually for testing) - ixc=0'
86 39 : citation = ''
87 : ! LDA,LSD
88 : case (1)
89 2175 : message = 'LDA: new Teter (4/93) with spin-polarized option - ixc=1'
90 2175 : citation = 'S. Goedecker, M. Teter, J. Huetter, PRB 54, 1703 (1996)' ! [[cite:Goedecker1996]]
91 : case (2)
92 181 : message = 'LDA: Perdew-Zunger-Ceperley-Alder - ixc=2'
93 181 : citation = 'J.P.Perdew and A.Zunger, PRB 23, 5048 (1981) ' ! [[cite:Perdew1981]]
94 : case (3)
95 297 : message = 'LDA: old Teter (4/91) fit to Ceperley-Alder data - ixc=3'
96 297 : citation = ''
97 : case (4)
98 2 : message = 'LDA: Wigner - ixc=4'
99 2 : citation = 'E.P.Wigner, Trans. Faraday Soc. 34, 67 (1938)' ! [[cite:Wigner1938]]
100 : case (5)
101 2 : message = 'LDA: Hedin-Lundqvist - ixc=5'
102 2 : citation = 'L.Hedin and B.I.Lundqvist, J. Phys. C4, 2064 (1971)' ! [[cite:Hedin1971]]
103 : case (6)
104 1 : message = 'LDA: "X-alpha" xc - ixc=6'
105 1 : citation = 'Slater J. C., Phys. Rev. 81, 385 (1951)' ! [[cite:Slater1951]]
106 : case (7)
107 665 : message = 'LDA: Perdew-Wang 92 LSD fit to Ceperley-Alder data - ixc=7'
108 665 : citation = 'J.P.Perdew and Y.Wang, PRB 45, 13244 (1992)' ! [[cite:Perdew1992a]]
109 : case (8)
110 1 : message = 'LDA: Perdew-Wang 92 LSD , exchange-only - ixc=8'
111 1 : citation = 'J.P.Perdew and Y.Wang, PRB 45, 13244 (1992)' ! [[cite:Perdew1992a]]
112 : case (9)
113 0 : message = 'LDA: Perdew-Wang 92 Ex+Ec_RPA energy - ixc=9'
114 0 : citation = 'J.P.Perdew and Y.Wang, PRB 45, 13244 (1992)' ! [[cite:Perdew1992]]
115 : case (10)
116 0 : message = 'LDA: RPA LSD energy (only the energy !!) - ixc=10'
117 0 : citation = ''
118 : ! GGA
119 : case (11)
120 614 : message = 'GGA: Perdew-Burke-Ernzerhof functional - ixc=11'
121 614 : citation = 'J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996)' ! [[cite:Perdew1996]]
122 : case (12)
123 7 : message = 'GGA: x-only Perdew-Burke-Ernzerhof functional - ixc=12'
124 7 : citation = 'J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996)' ! [[cite:Perdew1996]]
125 : case (13)
126 0 : message = 'GGA: LDA (ixc==7) energy, and the xc _potential_ is given by van Leeuwen-Baerends GGA - ixc=13'
127 0 : citation = 'R. van Leeuwen and E. J. Baerends PRA 49, 2421 (1994)' ! [[cite:VanLeeuwen1994]]
128 : case (14)
129 10 : message = 'GGA: revPBE functional - ixc=14'
130 10 : citation = 'Zhang and Yang, PRL 80, 890 (1998)' ! [[cite:Zhang1998]]
131 : case (15)
132 8 : message = 'GGA: RPBE functional - ixc=15'
133 8 : citation = 'Hammer, L. B. Hansen, and J. K. Norskov, PRB 59, 7413 (1999)' ! [[cite:Hammer1999]]
134 : case (16)
135 10 : message = 'GGA: HCTH93 functional - ixc=16'
136 10 : citation = 'F.A. Hamprecht, A.J. Cohen, D.J. Tozer, N.C. Handy, JCP 109, 6264 (1998)' ! [[cite:Hamprecht1998]]
137 : case (17)
138 2 : message = 'GGA: HCTH120 functional - ixc=17'
139 2 : citation = 'A.D. Boese, N.L. Doltsinis, N.C. Handy, and M. Sprik, JCP 112, 1670 (2000)' ! [[cite:Boese2000]]
140 : case (23)
141 21 : message = 'GGA: Wu Cohen functional - ixc=23'
142 21 : citation = 'Z. Wu and R. E. Cohen, PRB 73, 235116 (2006)' ! [[cite:Wu2006]]
143 : case (24)
144 2 : message = 'GGA: C09x exchange functional - ixc=24'
145 2 : citation = 'Valentino R. Cooper, PRB 81, 161104(R) (2010)' ! [[cite:Cooper2010]]
146 : case (26)
147 2 : message = 'GGA: HCTH147 functional - ixc=26'
148 2 : citation = 'A.D. Boese, N.L. Doltsinis, N.C. Handy, and M. Sprik, JCP 112, 1670 (2000)' ! [[cite:Boese2000]]
149 : case (27)
150 2 : message = 'GGA: HCTH407 functional - ixc=27'
151 2 : citation = 'A.D. Boese, and N.C. Handy, JCP 114, 5497 (2001)' ! [[cite:Boese2001]]
152 : ! Fermi-Amaldi
153 : case (20)
154 7 : message = 'Fermi-Amaldi correction - ixc=20'
155 7 : citation = ''
156 : case (21)
157 1 : message = 'Fermi-Amaldi correction with LDA(ixc=1) kernel - ixc=21'
158 1 : citation = ''
159 : case (22)
160 1 : message = 'Fermi-Amaldi correction with hybrid BPG kernel - ixc=22'
161 1 : citation = ''
162 : case (31)
163 2 : message = 'Meta-GGA fake1 - ixc=31'
164 2 : citation = ''
165 : case (32)
166 2 : message = 'Meta-GGA fake2 - ixc=32'
167 2 : citation = ''
168 : case (33)
169 4 : message = 'Meta-GGA fake3 - ixc=33'
170 4 : citation = ''
171 : case (34)
172 2 : message = 'Meta-GGA fake4 - ixc=34'
173 2 : citation = ''
174 : case (35)
175 0 : message = 'Meta-GGA fake5 - ixc=35'
176 0 : citation = ''
177 : case (40)
178 8 : message = 'Hartree-Fock with mixing coefficient alpha=1'
179 8 : citation = ''
180 : case (41)
181 53 : message = 'PBE0 with alpha=0.25'
182 53 : citation = ''
183 : case (42)
184 4 : message = 'modified PBE0 with alpha=0.33'
185 4 : citation = ''
186 : ! Temperature-dependant XC
187 : ! LDA
188 : case (50)
189 1 : message = 'LDA at finite T Ichimaru-Iyetomy-Tanaka - ixc=50'
190 1 : citation = 'Ichimaru S., Iyetomi H., Tanaka S., Phys. Rep. 149, 91-205 (1987)' ! [[cite:Ichimaru1987]]
191 : case (51)
192 2 : message = 'TLDA: corrKSDT Karasiev-Sjostrom-Dufty-Trickey - ixc=51'
193 2 : citation = 'V.V. Karasiev, T. Sjostrom, J. Dufty, and S.B. Trickey, PRL 112, 076403 (2014)' ! [[cite:Karasiev2014]]
194 : ! GGA
195 : case (60)
196 2 : message = 'TGGA: KDT16 Karasiev-Dufty-Trickey - ixc=60'
197 2 : citation = 'V.V. Karasiev, J.W. Dufty, and S.B. Trickey, PRL 120(7), 076401 (2018)' ! [[cite:Karasiev2018]]
198 : case default
199 0 : write(message,'(a,i0)')" echo_xc_name does not know how to handle ixc = ",ixc
200 4130 : ABI_WARNING(message)
201 : end select
202 :
203 : message = " Exchange-correlation functional for the present dataset will be:" // ch10 &
204 4130 : & // " " // trim(message)
205 :
206 4130 : l_citation=len_trim(citation)
207 4130 : citation = " Citation for XC functional:" // ch10 // " " // trim(citation)
208 :
209 4130 : call wrtout(ab_out,message,'COLL')
210 4130 : call wrtout(std_out,message,'COLL')
211 :
212 4130 : if(l_citation/=0)then
213 3710 : call wrtout(ab_out,citation,'COLL')
214 3710 : call wrtout(std_out,citation,'COLL')
215 : end if
216 :
217 4130 : message =' '
218 4130 : call wrtout(ab_out,message,'COLL')
219 4130 : call wrtout(std_out,message,'COLL')
220 :
221 : end if ! end libxc if
222 :
223 5281 : end subroutine echo_xc_name
224 : !!***
225 :
226 : !!****f* m_drivexc/xc_need_kden
227 : !! NAME
228 : !! xc_need_kden
229 : !!
230 : !! FUNCTION
231 : !! Check if kinetic energy density is used in XC functional
232 : !!
233 : !! INPUTS
234 : !! ixc= choice of exchange-correlation scheme
235 : !! [xc_funcs(2)]= <type(libxc_functional_type)>= optional - libXC functional(s)
236 : !!
237 : !! SOURCE
238 :
239 1227376 : logical function xc_need_kden(ixc,xc_funcs)
240 :
241 : !Arguments ------------------------------------
242 : integer,intent(in) :: ixc
243 : type(libxc_functional_type),intent(in),optional :: xc_funcs(2)
244 : ! *************************************************************************
245 :
246 1227376 : xc_need_kden=.false.
247 :
248 1227376 : if (ixc>=0) then
249 748077 : xc_need_kden=(ixc==31.or.ixc==34.or.ixc==35)
250 : else
251 479299 : if(present(xc_funcs)) xc_need_kden=libxc_functionals_needs_tau(xc_functionals=xc_funcs)
252 477754 : if(.not.present(xc_funcs)) xc_need_kden=libxc_functionals_needs_tau()
253 : end if
254 :
255 1227376 : end function xc_need_kden
256 : !!***
257 :
258 : !!****f* m_drivexc/has_kxc
259 : !! NAME
260 : !! has_kxc
261 : !!
262 : !! FUNCTION
263 : !! Given a XC functional (defined by ixc), return TRUE if Kxc (dVxc/drho) is avalaible.
264 : !!
265 : !! INPUTS
266 : !! ixc = internal code for xc functional
267 : !! [xc_funcs(2)]= <type(libxc_functional_type)> = optional - libXC set of functionals
268 : !!
269 : !! OUTPUT
270 : !!
271 : !! SOURCE
272 :
273 9053 : logical function has_kxc(ixc,xc_funcs)
274 :
275 : !Arguments -------------------------------
276 : integer, intent(in) :: ixc
277 : type(libxc_functional_type),intent(in),optional :: xc_funcs(2)
278 : ! *********************************************************************
279 :
280 9053 : has_kxc=.false.
281 :
282 9053 : if (ixc>=0) then
283 7302 : has_kxc=(ixc/=16.and.ixc/=17.and.ixc/=26.and.ixc/=27)
284 1751 : else if (ixc==-406.or.ixc==-427.or.ixc==-428.or.ixc==-456)then
285 : has_kxc=.true.
286 : else ! ixc<0 and not one of the allowed hybrids
287 1681 : if (present(xc_funcs)) then
288 1110 : has_kxc=libxc_functionals_has_kxc(xc_funcs)
289 : else
290 571 : has_kxc=libxc_functionals_has_kxc()
291 : end if
292 : end if
293 :
294 9053 : end function has_kxc
295 : !!***
296 :
297 : !!****f* m_drivexc/has_k3xc
298 : !! NAME
299 : !! has_k3xc
300 : !!
301 : !! FUNCTION
302 : !! Given a XC functional (defined by ixc), return TRUE if K3xc (d2Vxc/drho2) is avalaible.
303 : !!
304 : !! INPUTS
305 : !! ixc = internal code for xc functional
306 : !! [xc_funcs(2)]= <type(libxc_functional_type)> = optional - libXC set of functionals
307 : !!
308 : !!
309 : !! OUTPUT
310 : !!
311 : !! SOURCE
312 :
313 2098 : logical function has_k3xc(ixc,xc_funcs)
314 :
315 : !Arguments -------------------------------
316 : integer, intent(in) :: ixc
317 : type(libxc_functional_type),intent(in),optional :: xc_funcs(2)
318 : ! *********************************************************************
319 :
320 2098 : has_k3xc=.false.
321 :
322 2098 : if (ixc>=0) then
323 : has_k3xc=(ixc==0.or.ixc==3.or.(ixc>=7.and.ixc<=15).or. &
324 1518 : & ixc==23.or.ixc==24.or.ixc==41.or.ixc==42.or.ixc==1402000)
325 580 : else if (ixc==-406.or.ixc==-427.or.ixc==-428.or.ixc==-456)then
326 : has_k3xc=.false.
327 : else ! ixc<0 and not one of the allowed hybrids
328 571 : if (present(xc_funcs)) then
329 0 : has_k3xc=libxc_functionals_has_k3xc(xc_funcs)
330 : else
331 571 : has_k3xc=libxc_functionals_has_k3xc()
332 : end if
333 : end if
334 :
335 2098 : end function has_k3xc
336 : !!***
337 :
338 : !!****f* m_drivexc/check_kxc
339 : !! NAME
340 : !! check_kxc
341 : !!
342 : !! FUNCTION
343 : !! Given a XC functional (defined by ixc), check if Kxc and/or K3xc is avalaible.
344 : !!
345 : !! INPUTS
346 : !! ixc = internal code for xc functional
347 : !! optdriver=type of calculation (ground-state, response function, GW, ...)
348 : !! [check_k3xc]= optional ; check also k3xc availability
349 : !!
350 : !! OUTPUT
351 : !!
352 : !! SOURCE
353 :
354 2098 : subroutine check_kxc(ixc,optdriver,check_k3xc)
355 :
356 : !Arguments -------------------------------
357 : integer, intent(in) :: ixc,optdriver
358 : logical,intent(in),optional :: check_k3xc
359 :
360 : !Local variables -------------------------
361 : logical :: check_k3xc_,kxc_available,k3xc_available
362 : character(len=500) :: msg
363 : ! *********************************************************************
364 :
365 2098 : check_k3xc_=.false. ; if (present(check_k3xc)) check_k3xc_=check_k3xc
366 :
367 2098 : kxc_available=has_kxc(ixc)
368 2098 : k3xc_available=has_k3xc(ixc)
369 :
370 2098 : if (ixc>=0) then
371 1518 : if (.not.kxc_available) then
372 : write(msg,'(a,i0,3a)') &
373 0 : & 'The selected XC functional (ixc=',ixc,')',ch10,&
374 0 : & 'does not provide Kxc (dVxc/drho) !'
375 : end if
376 1518 : if (check_k3xc_.and.(.not.k3xc_available)) then
377 : write(msg,'(a,i0,3a)') &
378 0 : & 'The selected XC functional (ixc=',ixc,')',ch10,&
379 0 : & 'does not provide K3xc (d^2Vxc/drho^2) !'
380 : end if
381 : else ! ixc<0
382 580 : if (.not.kxc_available) then
383 : write(msg,'(a,i0,7a)') &
384 0 : & 'The selected XC functional (ixc=',ixc,'):',ch10,&
385 0 : & ' <<',trim(libxc_functionals_fullname()),'>>',ch10,&
386 0 : & 'does not provide Kxc (dVxc/drho) !'
387 : end if
388 580 : if (check_k3xc_.and.(.not.k3xc_available)) then
389 : write(msg,'(a,i0,7a)') &
390 0 : & 'The selected XC functional (ixc=',ixc,'):',ch10,&
391 0 : & ' <<',trim(libxc_functionals_fullname()),'>>',ch10,&
392 0 : & 'does not provide K3xc (d^2Vxc/d^2rho) !'
393 : end if
394 : end if
395 :
396 2098 : if (.not.kxc_available) then
397 0 : write(msg,'(7a)') trim(msg),ch10,&
398 0 : & 'However, with the current input options, ABINIT needs Kxc.',ch10,&
399 0 : & '>Possible action:',ch10,&
400 0 : & 'Change the XC functional in psp file or input file.'
401 0 : if (optdriver==0) then
402 0 : write(msg,'(13a)') trim(msg),ch10,&
403 0 : & '>Possible action (2):',ch10,&
404 0 : & 'If you are using density mixing for the SCF cycle',ch10,&
405 0 : & '(iscf>=10, which is the default for PAW),',ch10,&
406 0 : & 'change to potential mixing (iscf=7, for instance).',ch10,&
407 0 : & '>Possible action (3):',ch10,&
408 0 : & 'Switch to another value of densfor_pred (=5, for instance).'
409 : end if
410 0 : ABI_ERROR(msg)
411 2098 : else if (check_k3xc_.and.(.not.k3xc_available)) then
412 0 : write(msg,'(13a)') trim(msg),ch10,&
413 0 : & 'However, with the current input options, ABINIT needs K3xc.',ch10,&
414 0 : & '>Possible actions:',ch10,&
415 0 : & '- Recompile libXC using --enable-kxc.',ch10,&
416 0 : & ' or',ch10,&
417 0 : & '- Change the XC functional in psp file or input file:',ch10,&
418 0 : & ' use one of the internal LDA (ixc=3, 7 to 15, 23, 24).'
419 0 : ABI_ERROR(msg)
420 : end if
421 :
422 2098 : end subroutine check_kxc
423 : !!***
424 :
425 : !!****f* m_drivexc/size_dvxc
426 : !! NAME
427 : !! size_dvxc
428 : !!
429 : !! FUNCTION
430 : !! Give the sizes of the several arrays involved in exchange-correlation calculation
431 : !! needed to allocated them for the drivexc routine
432 : !!
433 : !! INPUTS
434 : !! ixc= choice of exchange-correlation scheme
435 : !! order= gives the maximal derivative of Exc computed.
436 : !! 1=usual value (return exc and vxc)
437 : !! 2=also computes the kernel (return exc,vxc,kxc)
438 : !! -2=like 2, except (to be described)
439 : !! 3=also computes the derivative of the kernel (return exc,vxc,kxc,k3xc)
440 : !! nspden= number of spin components
441 : !! [xc_funcs(2)]= <type(libxc_functional_type)>= optional - libXC functional(s)
442 : !! [add_tfw]= optional flag controling the addition of Weiszacker gradient correction to Thomas-Fermi XC energy
443 : !!
444 : !! OUTPUT
445 : !! --- All optionals
446 : !! [usegradient]= [flag] 1 if the XC functional needs the gradient of the density (grho2_updn)
447 : !! [uselaplacian]= [flag] 1 if the XC functional needs the laplacian of the density (lrho_updn)
448 : !! [usekden]= [flag] 1 if the XC functional needs the kinetic energy density (lrho_updn)
449 : !! [nvxcgrho]= size of the array dvxcdgr(npts,nvxcgrho) (derivative of Exc wrt to gradient)
450 : !! [nvxclrho]= size of the array dvxclpl(npts,nvxclrho) (derivative of Exc wrt to laplacian)
451 : !! [nvxctau]= size of the array dvxctau(npts,nvxctau) (derivative of Exc wrt to kin. ener. density)
452 : !! [ndvxc]= size of the array dvxc(npts,ndvxc) (second derivatives of Exc wrt to density and gradient)
453 : !! [nd2vxc]= size of the array d2vxc(npts,nd2vxc) (third derivatives of Exc wrt density)
454 : !!
455 : !! SOURCE
456 :
457 1219696 : subroutine size_dvxc(ixc,order,nspden,&
458 : & usegradient,uselaplacian,usekden,&
459 : & nvxcgrho,nvxclrho,nvxctau,ndvxc,nd2vxc,&
460 : & add_tfw,xc_funcs) ! Optional
461 :
462 : !Arguments----------------------
463 : integer,intent(in) :: ixc,nspden,order
464 : integer,intent(out),optional :: nvxcgrho,nvxclrho,nvxctau,ndvxc,nd2vxc
465 : integer,intent(out),optional :: usegradient,uselaplacian,usekden
466 : logical, intent(in),optional :: add_tfw
467 : type(libxc_functional_type),intent(in),optional :: xc_funcs(2)
468 :
469 : !Local variables----------------
470 : logical :: libxc_has_kxc,libxc_has_k3xc,libxc_isgga,libxc_ismgga,libxc_ishybrid,my_add_tfw
471 : logical :: need_gradient,need_laplacian,need_kden
472 : ! *************************************************************************
473 :
474 : !Several flags
475 1219696 : my_add_tfw=.false.;if (present(add_tfw)) my_add_tfw=add_tfw
476 1219696 : libxc_isgga=.false. ; libxc_ismgga=.false. ; libxc_ishybrid=.false.
477 1219696 : if(ixc<0)then
478 477844 : if(present(xc_funcs))then
479 1545 : libxc_has_kxc=libxc_functionals_has_kxc(xc_funcs)
480 1545 : libxc_has_k3xc=libxc_functionals_has_k3xc(xc_funcs)
481 1545 : libxc_isgga=libxc_functionals_isgga(xc_functionals=xc_funcs)
482 1545 : libxc_ismgga=libxc_functionals_ismgga(xc_functionals=xc_funcs)
483 1545 : libxc_ishybrid=libxc_functionals_is_hybrid(xc_functionals=xc_funcs)
484 : else
485 476299 : libxc_has_kxc=libxc_functionals_has_kxc()
486 476299 : libxc_has_k3xc=libxc_functionals_has_k3xc()
487 476299 : libxc_isgga=libxc_functionals_isgga()
488 476299 : libxc_ismgga=libxc_functionals_ismgga()
489 476299 : libxc_ishybrid=libxc_functionals_is_hybrid()
490 : end if
491 : end if
492 :
493 : !Do we use the gradient?
494 : need_gradient=((ixc>=11.and.ixc<=17).or.(ixc==23.or.ixc==24).or. &
495 : & (ixc==26.or.ixc==27).or.(ixc>=31.and.ixc<=35).or. &
496 1219696 : & (ixc==41.or.ixc==42).or.ixc==1402000.or.ixc==60)
497 1219696 : if (ixc<0.and.(libxc_isgga.or.libxc_ismgga.or.libxc_ishybrid)) need_gradient=.true.
498 1219696 : if (my_add_tfw) need_gradient=.true.
499 1983951 : if (present(usegradient)) usegradient=merge(1,0,need_gradient)
500 :
501 : !Do we use the laplacian?
502 1219696 : need_laplacian=(ixc==32.or.ixc==35)
503 1219696 : if (ixc<0) then
504 477844 : if(present(xc_funcs)) need_laplacian=libxc_functionals_needs_laplacian(xc_functionals=xc_funcs)
505 476299 : if(.not.present(xc_funcs)) need_laplacian=libxc_functionals_needs_laplacian()
506 : end if
507 2403617 : if (present(uselaplacian)) uselaplacian=merge(1,0,need_laplacian)
508 :
509 : !Do we use the kinetic energy density?
510 1219696 : if(present(xc_funcs)) need_kden=xc_need_kden(ixc,xc_funcs)
511 1212300 : if(.not.present(xc_funcs)) need_kden=xc_need_kden(ixc)
512 1219696 : usekden=merge(1,0,need_kden)
513 :
514 : !First derivative(s) of XC functional wrt gradient of density
515 1219696 : if (present(nvxcgrho)) then
516 1163469 : nvxcgrho=0
517 1163469 : if (abs(order)>=1) then
518 1163469 : if (need_gradient) nvxcgrho=3
519 1163469 : if (ixc==13) nvxcgrho=0
520 1163469 : if (ixc==16.or.ixc==17.or.ixc==26.or.ixc==27) nvxcgrho=2
521 : end if
522 : end if
523 :
524 : !First derivative(s) of XC functional wrt laplacian of density
525 1219696 : if (present(nvxclrho)) then
526 1163469 : nvxclrho=0
527 1163469 : if (abs(order)>=1) then
528 1163469 : if (need_laplacian) nvxclrho=min(nspden,2)
529 : end if
530 : end if
531 :
532 : !First derivative(s) of XC functional wrt kinetic energy density
533 1219696 : if (present(nvxctau)) then
534 1163469 : nvxctau=0
535 1163469 : if (abs(order)>=1) then
536 1163469 : if (need_kden) nvxctau=min(nspden,2)
537 : end if
538 : end if
539 :
540 : !Second derivative(s) of XC functional wrt density
541 1219696 : if (present(ndvxc)) then
542 1163469 : ndvxc=0
543 1163469 : if (abs(order)>=2) then
544 : if (ixc==1.or.ixc==7.or.ixc==8.or.ixc==9.or.ixc==10.or.ixc==13.or. &
545 160784 : & ixc==21.or.ixc==22) then
546 47262 : ndvxc=min(nspden,2)+1
547 113522 : else if ((ixc>=2.and.ixc<=6).or.(ixc>=31.and.ixc<=35).or.&
548 : & (ixc==50.or.ixc==51.or.ixc==60)) then
549 53818 : ndvxc=1
550 59704 : else if (ixc==12.or.ixc==24) then
551 4 : ndvxc=8
552 : else if (ixc==11.or.ixc==12.or.ixc==14.or.ixc==15.or. &
553 59700 : & ixc==23.or.ixc==41.or.ixc==42.or.ixc==1402000) then
554 13958 : ndvxc=15
555 45742 : else if (ixc<0) then
556 45712 : if (libxc_has_kxc.or.ixc==-406.or.ixc==-427.or.ixc==-428.or.ixc==-456) then
557 45712 : ndvxc=2*min(nspden,2)+1 ; if (order==-2) ndvxc=2
558 45712 : if (need_gradient) ndvxc=15 ! This is for GGA, but also for mGGA
559 : ! (we dont consider derivatives wrt Tau or Laplacian)
560 : end if
561 : end if
562 : end if
563 : end if
564 :
565 : !Third derivative(s) of XC functional wrt density
566 1219696 : if (present(nd2vxc)) then
567 1163469 : nd2vxc=0
568 1163469 : if (abs(order)>=3) then
569 : if (ixc==3.or.(ixc>=11.and.ixc<=15.and.ixc/=13).or. &
570 2588 : & ixc==23.or.ixc==24.or.ixc==41.or.ixc==42) then
571 39 : nd2vxc=1
572 2549 : else if ((ixc>=7.and.ixc<=10).or.ixc==13.or.ixc==1402000) then
573 2542 : nd2vxc=3*min(nspden,2)-2
574 7 : else if (ixc<0) then
575 4 : if (libxc_has_k3xc) then
576 4 : if (.not.need_gradient) nd2vxc=3*min(nspden,2)-2
577 : end if
578 : end if
579 : end if
580 : end if
581 :
582 1219696 : end subroutine size_dvxc
583 : !!***
584 :
585 : !!****f* m_drivexc/xcmult
586 : !! NAME
587 : !! xcmult
588 : !!
589 : !! FUNCTION
590 : !! In the case of GGA, multiply the different gradient of spin-density
591 : !! by the derivative of the XC functional with respect
592 : !! to the norm of the gradient, then divide it by the norm of the gradient
593 : !!
594 : !! INPUTS
595 : !! depsxc(nfft,nspgrad)=derivative of Exc with respect to the (spin-)density,
596 : !! or to the norm of the gradient of the (spin-)density,
597 : !! further divided by the norm of the gradient of the (spin-)density
598 : !! The different components of depsxc will be
599 : !! for nspden=1, depsxc(:,1)=d(rho.exc)/d(rho)
600 : !! and if ngrad=2, depsxc(:,2)=1/2*1/|grad rho_up|*d(rho.exc)/d(|grad rho_up|)
601 : !! + 1/|grad rho|*d(rho.exc)/d(|grad rho|)
602 : !! (do not forget : |grad rho| /= |grad rho_up| + |grad rho_down|
603 : !! for nspden=2, depsxc(:,1)=d(rho.exc)/d(rho_up)
604 : !! depsxc(:,2)=d(rho.exc)/d(rho_down)
605 : !! and if ngrad=2, depsxc(:,3)=1/|grad rho_up|*d(rho.exc)/d(|grad rho_up|)
606 : !! depsxc(:,4)=1/|grad rho_down|*d(rho.exc)/d(|grad rho_down|)
607 : !! depsxc(:,5)=1/|grad rho|*d(rho.exc)/d(|grad rho|)
608 : !! nfft=(effective) number of FFT grid points (for this processor)
609 : !! ngrad = must be 2
610 : !! nspden=number of spin-density components
611 : !! nspgrad=number of spin-density and spin-density-gradient components
612 : !!
613 : !! OUTPUT
614 : !! (see side effects)
615 : !!
616 : !! SIDE EFFECTS
617 : !! rhonow(nfft,nspden,ngrad*ngrad)=
618 : !! at input :
619 : !! electron (spin)-density in real space and its gradient,
620 : !! either on the unshifted grid (if ishift==0,
621 : !! then equal to rhor), or on the shifted grid
622 : !! rhonow(:,:,1)=electron density in electrons/bohr**3
623 : !! rhonow(:,:,2:4)=gradient of electron density in el./bohr**4
624 : !! at output :
625 : !! rhonow(:,:,2:4) has been multiplied by the proper factor,
626 : !! described above.
627 : !!
628 : !! SOURCE
629 :
630 203251 : subroutine xcmult (depsxc,nfft,ngrad,nspden,nspgrad,rhonow)
631 :
632 : !Arguments ------------------------------------
633 : !scalars
634 : integer,intent(in) :: nfft,ngrad,nspden,nspgrad
635 : !arrays
636 : real(dp),intent(in) :: depsxc(nfft,nspgrad)
637 : real(dp),intent(inout) :: rhonow(nfft,nspden,ngrad*ngrad)
638 :
639 : !Local variables-------------------------------
640 : !scalars
641 : integer :: idir,ifft
642 : real(dp) :: rho_tot,rho_up
643 : ! *************************************************************************
644 :
645 813004 : do idir=1,3
646 :
647 813004 : if(nspden==1)then
648 : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(depsxc,idir,nfft,rhonow)
649 1614874776 : do ifft=1,nfft
650 1614874776 : rhonow(ifft,1,1+idir)=rhonow(ifft,1,1+idir)*depsxc(ifft,2)
651 : end do
652 :
653 : else
654 :
655 : ! In the spin-polarized case, there are more factors to take into account
656 : !$OMP PARALLEL DO PRIVATE(ifft,rho_tot,rho_up) SHARED(depsxc,idir,nfft,rhonow)
657 253116483 : do ifft=1,nfft
658 253049166 : rho_tot=rhonow(ifft,1,1+idir)
659 253049166 : rho_up =rhonow(ifft,2,1+idir)
660 253049166 : rhonow(ifft,1,1+idir)=rho_up *depsxc(ifft,3) + rho_tot*depsxc(ifft,5)
661 253116483 : rhonow(ifft,2,1+idir)=(rho_tot-rho_up)*depsxc(ifft,4)+ rho_tot*depsxc(ifft,5)
662 : end do
663 :
664 : end if ! nspden==1
665 :
666 : end do ! End loop on directions
667 :
668 203251 : end subroutine xcmult
669 : !!***
670 :
671 : !!****f* m_drivexc/mkdenpos
672 : !! NAME
673 : !! mkdenpos
674 : !!
675 : !! FUNCTION
676 : !! Make a density positive everywhere:
677 : !! when the density (or spin-density) is smaller than xc_denpos,
678 : !! set it to the value of xc_denpos
679 : !!
680 : !! INPUTS
681 : !! nfft=(effective) number of FFT grid points (for this processor)
682 : !! nspden=number of spin-density components (max. 2)
683 : !! option=0 if density rhonow is stored as (up,dn)
684 : !! 1 if density rhonow is stored as (up+dn,up)
685 : !! Active only when nspden=2
686 : !! xc_denpos= lowest allowed density (usually for the computation of the XC functionals)
687 : !!
688 : !! OUTPUT
689 : !! (see side effects)
690 : !!
691 : !! SIDE EFFECTS
692 : !! Input/output
693 : !! iwarn=At input: iwarn=0 a warning will be printed when rho is negative
694 : !! iwarn>0 no warning will be printed out
695 : !! At output: iwarn is increased by 1
696 : !! rhonow(nfft,nspden)=electron (spin)-density in real space,
697 : !! either on the unshifted grid (if ishift==0,
698 : !! then equal to rhor),or on the shifted grid
699 : !!
700 : !! NOTES
701 : !! At this stage, rhonow(:,1:nspden) contains the density in real space,
702 : !! on the unshifted or shifted grid. Now test for negative densities
703 : !! Note that, ignoring model core charge, as long as boxcut>=2
704 : !! the shifted density is derivable from the square of a Fourier
705 : !! interpolated charge density => CANNOT go < 0.
706 : !! However, actually can go < 0 to within machine precision;
707 : !! do not print useless warnings in this case, just fix it.
708 : !! Fourier interpolated core charge can go < 0 due to Gibbs
709 : !! oscillations; could avoid this by recomputing the model core
710 : !! charge at the new real space grid points (future work).
711 : !!
712 : !! SOURCE
713 :
714 613857 : subroutine mkdenpos(iwarn,nfft,nspden,option,rhonow,xc_denpos)
715 :
716 : !Arguments ------------------------------------
717 : !scalars
718 : integer,intent(in) :: nfft,nspden,option
719 : integer,intent(inout) :: iwarn
720 : real(dp),intent(in) :: xc_denpos
721 : !arrays
722 : real(dp),intent(inout) :: rhonow(nfft,nspden)
723 :
724 : !Local variables-------------------------------
725 : !scalars
726 : integer :: ifft,ispden,numneg
727 : real(dp) :: rhotmp,worst
728 : character(len=600) :: message
729 : !arrays
730 : real(dp) :: rho(2)
731 : ! *************************************************************************
732 :
733 613857 : numneg=0
734 613857 : worst=zero
735 :
736 613857 : if(nspden==1)then
737 :
738 : ! Non spin-polarized
739 : !$OMP PARALLEL DO PRIVATE(ifft,rhotmp) REDUCTION(MIN:worst) REDUCTION(+:numneg) SHARED(nfft,rhonow)
740 1765858080 : do ifft=1,nfft
741 1765377882 : rhotmp=rhonow(ifft,1)
742 1765858080 : if(rhotmp<xc_denpos)then
743 45426772 : if(rhotmp<-xc_denpos)then
744 : ! This case is probably beyond machine precision considerations
745 43772719 : worst=min(worst,rhotmp)
746 43772719 : numneg=numneg+1
747 : end if
748 45426772 : rhonow(ifft,1)=xc_denpos
749 : end if
750 : end do
751 133659 : else if (nspden==2) then
752 :
753 : ! Spin-polarized
754 :
755 : ! rhonow is stored as (up,dn)
756 133659 : if (option==0) then
757 :
758 : !$OMP PARALLEL DO PRIVATE(ifft,ispden,rho,rhotmp) REDUCTION(MIN:worst) REDUCTION(+:numneg) &
759 : !$OMP&SHARED(nfft,nspden,rhonow)
760 122881566 : do ifft=1,nfft
761 : ! For polarized case, rho(1) is spin-up density, rho(2) is spin-down density
762 122757235 : rho(1)=rhonow(ifft,1)
763 122757235 : rho(2)=rhonow(ifft,2)
764 368396036 : do ispden=1,nspden
765 368271705 : if (rho(ispden)<xc_denpos) then
766 41302 : if (rho(ispden)<-xc_denpos) then
767 : ! This case is probably beyond machine precision considerations
768 28 : worst=min(worst,rho(ispden))
769 28 : numneg=numneg+1
770 : end if
771 41302 : rhonow(ifft,ispden)=xc_denpos
772 : end if
773 : end do
774 : end do
775 :
776 : ! rhonow is stored as (up+dn,up)
777 9328 : else if (option==1) then
778 :
779 : !$OMP PARALLEL DO PRIVATE(ifft,ispden,rho,rhotmp) &
780 : !$OMP&REDUCTION(MIN:worst) REDUCTION(+:numneg) &
781 : !$OMP&SHARED(nfft,nspden,rhonow)
782 206058516 : do ifft=1,nfft
783 : ! For polarized case, rho(1) is spin-up density, rho(2) is spin-down density
784 206049188 : rho(1)=rhonow(ifft,2)
785 206049188 : rho(2)=rhonow(ifft,1)-rho(1)
786 618156892 : do ispden=1,nspden
787 618147564 : if (rho(ispden)<xc_denpos) then
788 17365657 : if (rho(ispden)<-xc_denpos) then
789 : ! This case is probably beyond machine precision considerations
790 1546829 : worst=min(worst,rho(ispden))
791 1546829 : numneg=numneg+1
792 : end if
793 17365657 : rho(ispden)=xc_denpos
794 17365657 : rhonow(ifft,1)=rho(1)+rho(2)
795 17365657 : rhonow(ifft,2)=rho(1)
796 : end if
797 : end do
798 : end do
799 :
800 : end if ! option
801 :
802 : else
803 0 : ABI_BUG('nspden>2 not allowed !')
804 : end if ! End choice between non-spin polarized and spin-polarized.
805 :
806 613857 : if (numneg>0) then
807 1741 : if (iwarn==0) then
808 : write(message,'(a,i0,a,a,a,es10.2,a,e10.2,11a)')&
809 1336 : & 'Density went too small (lower than xc_denpos) at ',numneg,' points',ch10,&
810 1336 : & 'and was set to xc_denpos = ',xc_denpos,'. Lowest was ',worst,'.',ch10,&
811 1336 : & 'This might be due to (1) too low boxcut or (2) too low ecut for',ch10,&
812 1336 : & ' pseudopotential core charge, or (3) too low ecut for estimated initial density.',ch10,&
813 1336 : & ' Possible workarounds : increase ecut, or define the input variable densty,',ch10,&
814 1336 : & ' with a value larger than the guess for the decay length, or initialize your,',ch10,&
815 2672 : & ' density with a preliminary LDA or GGA-PBE if you are using a more exotic xc functional.'
816 1336 : ABI_WARNING(message)
817 : end if
818 1741 : iwarn=iwarn+1
819 : end if
820 :
821 613857 : end subroutine mkdenpos
822 : !!***
823 :
824 : !!****f* m_drivexc/drivexc
825 : !! NAME
826 : !! drivexc
827 : !!
828 : !! FUNCTION
829 : !! Driver of XC functionals. Treat spin-polarized as well as non-spin-polarized.
830 : !! Treat local approximations, GGAs, meta-GGAs or hybrid functionals.
831 : !! Optionally, deliver the XC kernel, or even the derivative
832 : !! of the XC kernel (the third derivative of the XC energy)
833 : !!
834 : !! INPUTS
835 : !! ixc=index of the XC functional
836 : !! xclevel=XC functional level (lda, gga, etc...)
837 : !! usegradient=[flag] 1 if the XC functional depends on density gradient (grho2_updn)
838 : !! uselaplacian=[flag] 1 if the XC functional depends on density laplacian (lrho_updn)
839 : !! usekden=[flag] 1 if the XC functional depends on kinetic energy density (tau_updn)
840 : !! order=gives the maximal derivative of Exc computed.
841 : !! 1=usual value (return exc and vxc)
842 : !! 2=also computes the kernel (return exc,vxc,kxc)
843 : !! -2=like 2, except (to be described)
844 : !! 3=also computes the derivative of the kernel (return exc,vxc,kxc,k3xc)
845 : !! npts=number of real space points on which the density is provided
846 : !! nspden=number of spin-density components (1 or 2)
847 : !! nvxcgrho=number of components of 1st-derivative of Exc wrt density gradient (nvxcgrho)
848 : !! nvxclrho=number of components of 1st-derivative of Exc wrt density laplacian (nvxclrho)
849 : !! nvxctau=number of components of 1st-derivative of Exc wrt kinetic energy density (nvxctau)
850 : !! ndvxc=number of components of 1st-derivative of Vxc (dvxc)
851 : !! nd2vxc=number of components of 2nd-derivative of Vxc (d2vxc)
852 : !! el_temp=electronic temperature (hartree)
853 : !! rho_updn(npts,nspden)=spin-up and spin-down densities
854 : !! In the calling routine, spin-down density must be equal to spin-up density.
855 : !! If nspden=1, only spin-up density must be given (half the total density).
856 : !! If nspden=2, spin-up and spin-down densities must be given.
857 : !! === Optional input arguments ===
858 : !! [grho2_updn(npts,(2*nspden-1)*usegradient)]=the square of the gradients
859 : !! of spin-up, spin-down, and total density.
860 : !! If nspden=1, only the square of the gradient of the spin-up density must be given.
861 : !! In the calling routine, the square of the gradient of the spin-down density must be equal
862 : !! to the square of the gradient of the spin-up density, and both must be equal to
863 : !! one-quarter of the square of the gradient of the total density.
864 : !! If nspden=2, the square of the gradients of spin-up, spin-down, and total density must be given.
865 : !! Note that the square of the gradient of the total density is usually NOT related to
866 : !! the square of the gradient of the spin-up and spin-down densities, because the gradients
867 : !! are not usually aligned. This is not the case when nspden=1.
868 : !! [lrho_updn(npts,nspden*uselaplacian)]=the Laplacian of spin-up and spin-down densities.
869 : !! If nspden=1, only the spin-up Laplacian density must be given and must
870 : !! be equal to the spin-up Laplacian density.
871 : !! If nspden=2, the Laplacian of spin-up and spin-down densities must be given.
872 : !! [tau_updn(npts,nspden*usekden)]=the spin-up and spin-down kinetic energy densities.
873 : !! If nspden=1, only the spin-up kinetic energy density must be given and must
874 : !! be equal to the half the total kinetic energy density.
875 : !! If nspden=2, the spin-up and spin-down kinetic energy densities must be given.
876 : !! [exexch]=choice of <<<local>>> exact exchange. Active if exexch=3 (only for GGA, and NOT for libxc)
877 : !! [hyb_mixing]= mixing parameter for the native PBEx functionals (ixc=41 and 42)
878 : !! [xc_funcs(2)]= <type(libxc_functional_type)>: libxc XC functionals.
879 : !!
880 : !! OUTPUT
881 : !! exc(npts)=exchange-correlation energy density (hartree)
882 : !! tsxc(npts)=exchange-correlation entropy energy density (hartree)
883 : !! vxcrho(npts,nspden)= (d($\rho$*exc)/d($\rho_up$)) (hartree)
884 : !! and (d($\rho$*exc)/d($\rho_down$)) (hartree)
885 : !! === Optional output arguments ===
886 : !! [vxcgrho(npts,nvxcgrho)]=1st-derivative of the xc energy wrt density gradient>
887 : !! = 1/$|grad \rho_up|$ (d($\rho$*exc)/d($|grad \rho_up|$))
888 : !! 1/$|grad \rho_dn|$ (d($\rho$*exc)/d($|grad \rho_dn|$))
889 : !! 1/$|grad \rho|$ (d($\rho$*exc)/d($|grad \rho|$))
890 : !! [vxclrho(npts,nvxclrho)]=1st-derivative of the xc energy wrt density laplacian.
891 : !! = d($\rho$*exc)/d($\lrho_up$)
892 : !! d($\rho$*exc)/d($\lrho_down$)
893 : !! [vxctau(npts,nvxctau)]=1st-derivative of the xc energy wrt kinetic energy density.
894 : !! = d($\rho$*exc)/d($\tau_up$)
895 : !! d($\rho$*exc)/d($\tau_down$)
896 : !! [dvxc(npts,ndvxc)]=partial second derivatives of the XC energy
897 : !! === Only if abs(order)>1 ===
898 : !! In case of local energy functional (option=1,-1 or 3):
899 : !! dvxc(npts,1+nspden)=
900 : !! if(nspden=1 .and. order==2): dvxci(:,1)=dvxc/d$\rho$ , dvxc(:,2) empty
901 : !! if(nspden=1 .and. order==-2): also compute dvxci(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$
902 : !! if(nspden=2): dvxc(:,1)=dvxc($\uparrow$)/d$\rho(\downarrow)$,
903 : !! dvxc(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$,
904 : !! dvxc(:,3)=dvxc($\downarrow$)/d$\rho(\downarrow)$
905 : !! In case of gradient corrected functional (option=2,-2, 4, -4, 5, 6, 7):
906 : !! dvxc(npts,15)=
907 : !! dvxc(:,1)= d2Ex/drho_up drho_up
908 : !! dvxc(:,2)= d2Ex/drho_dn drho_dn
909 : !! dvxc(:,3)= dEx/d(abs(grad(rho_up))) / abs(grad(rho_up))
910 : !! dvxc(:,4)= dEx/d(abs(grad(rho_dn))) / abs(grad(rho_dn))
911 : !! dvxc(:,5)= d2Ex/d(abs(grad(rho_up))) drho_up / abs(grad(rho_up))
912 : !! dvxc(:,6)= d2Ex/d(abs(grad(rho_dn))) drho_dn / abs(grad(rho_dn))
913 : !! dvxc(:,7)= 1/abs(grad(rho_up)) * d/d(abs(grad(rho_up)) (dEx/d(abs(grad(rho_up))) /abs(grad(rho_up)))
914 : !! dvxc(:,8)= 1/abs(grad(rho_dn)) * d/d(abs(grad(rho_dn)) (dEx/d(abs(grad(rho_dn))) /abs(grad(rho_dn)))
915 : !! dvxc(:,9)= d2Ec/drho_up drho_up
916 : !! dvxc(:,10)=d2Ec/drho_up drho_dn
917 : !! dvxc(:,11)=d2Ec/drho_dn drho_dn
918 : !! dvxc(:,12)=dEc/d(abs(grad(rho))) / abs(grad(rho))
919 : !! dvxc(:,13)=d2Ec/d(abs(grad(rho))) drho_up / abs(grad(rho))
920 : !! dvxc(:,14)=d2Ec/d(abs(grad(rho))) drho_dn / abs(grad(rho))
921 : !! dvxc(:,15)=1/abs(grad(rho)) * d/d(abs(grad(rho)) (dEc/d(abs(grad(rho))) /abs(grad(rho)))
922 : !! Note about mGGA: 2nd derivatives involving Tau or Laplacian are not output
923 : !! [d2vxc(npts,nd2vxc)]=second derivative of the XC potential=3rd order derivative of XC energy
924 : !! === Only if abs(order)>1 ===
925 : !! === At present only available for LDA ===
926 : !! if nspden=1 d2vxc(npts,1)=second derivative of the XC potential=3rd order derivative of energy
927 : !! if nspden=2 d2vxc(npts,1), d2vxc(npts,2), d2vxc(npts,3), d2vxc(npts,4) (3rd derivative of energy)
928 : !!
929 : !! SOURCE
930 :
931 872659 : subroutine drivexc(ixc,order,npts,nspden,usegradient,uselaplacian,usekden,&
932 872659 : & rho_updn,exc,tsxc,vxcrho,nvxcgrho,nvxclrho,nvxctau,ndvxc,nd2vxc,el_temp, & ! mandatory arguments
933 4678387 : & grho2_updn,vxcgrho,lrho_updn,vxclrho,tau_updn,vxctau,dvxc,d2vxc, & ! optional arguments
934 : & exexch,hyb_mixing,xc_funcs) ! optional parameters
935 :
936 : !Arguments ------------------------------------
937 : !scalars
938 : integer,intent(in) :: ixc,npts,nspden
939 : integer,intent(in) :: ndvxc,nd2vxc,nvxcgrho,nvxclrho,nvxctau,order
940 : integer,intent(in) :: usegradient,uselaplacian,usekden
941 : real(dp),intent(in) :: el_temp
942 : integer,intent(in),optional :: exexch
943 : real(dp),intent(in),optional :: hyb_mixing
944 : !arrays
945 : real(dp),intent(in) :: rho_updn(npts,nspden)
946 : real(dp),intent(in),optional :: grho2_updn(npts,(2*nspden-1)*usegradient)
947 : real(dp),intent(in),optional :: lrho_updn(npts,nspden*uselaplacian),tau_updn(npts,nspden*usekden)
948 : real(dp),intent(out) :: exc(npts),tsxc(npts),vxcrho(npts,nspden)
949 : real(dp),intent(out),optional :: dvxc(npts,ndvxc),d2vxc(npts,nd2vxc)
950 : real(dp),intent(out),optional :: vxcgrho(npts,nvxcgrho),vxclrho(npts,nvxclrho),vxctau(npts,nvxctau)
951 : type(libxc_functional_type),intent(inout),optional :: xc_funcs(2)
952 :
953 : !Local variables-------------------------------
954 : !scalars
955 : integer :: ispden,ixc_from_lib,ixc1,ixc2,ndvxc_x
956 : integer :: my_exexch,need_ndvxc,need_nd2vxc,need_nvxcgrho,need_nvxclrho,need_nvxctau
957 : integer :: need_gradient,need_laplacian,need_kden,optpbe
958 : logical :: has_gradient,has_laplacian,has_kden,libxc_test
959 : real(dp) :: alpha,beta,my_hyb_mixing
960 : real(dp),parameter :: rsfac=0.6203504908994000e0_dp
961 : character(len=500) :: message
962 : !arrays
963 1745318 : real(dp),allocatable :: exci_rpa(:),rhotot(:),rspts(:),vxci_rpa(:,:),zeta(:)
964 872659 : real(dp),allocatable :: exc_c(:),exc_x(:),vxcrho_c(:,:),vxcrho_x(:,:)
965 872659 : real(dp),allocatable :: d2vxc_c(:,:),d2vxc_x(:,:),dvxc_c(:,:),dvxc_x(:,:)
966 872659 : real(dp),allocatable :: vxcgrho_x(:,:)
967 4363295 : type(libxc_functional_type) :: xc_funcs_vwn3(2),xc_funcs_lyp(2)
968 : ! *************************************************************************
969 :
970 : !optional arguments
971 872659 : my_exexch=0;if(present(exexch)) my_exexch=exexch
972 872659 : my_hyb_mixing=0
973 872659 : if (ixc==41) my_hyb_mixing=quarter
974 872659 : if (ixc==42) my_hyb_mixing=third
975 872659 : if (present(hyb_mixing)) my_hyb_mixing=hyb_mixing
976 :
977 : ! =================================================
978 : ! == Compatibility tests ==
979 : ! =================================================
980 :
981 : !Check libXC initialization
982 872659 : if (ixc<0 .or. ixc==1402) then
983 358389 : libxc_test=libxc_functionals_check(stop_if_error=.true.)
984 : end if
985 :
986 : ! Check libXC consistency between ixc passed in input
987 : ! and the one used to initialize the libXC library
988 872659 : if (ixc<0) then
989 : ixc_from_lib=libxc_functionals_ixc()
990 358389 : if (present(xc_funcs)) then
991 1024 : ixc_from_lib=libxc_functionals_ixc(xc_functionals=xc_funcs)
992 : else
993 357365 : ixc_from_lib=libxc_functionals_ixc()
994 : end if
995 358389 : if (ixc/=ixc_from_lib) then
996 : write(message, '(a,i0,2a,i0)')&
997 0 : & 'The value of ixc specified in input, ixc = ',ixc,ch10,&
998 0 : & 'differs from the one used to initialize the functional ',ixc_from_lib
999 0 : ABI_BUG(message)
1000 : end if
1001 : end if
1002 :
1003 : !Check value of order
1004 872659 : if( (order<1.and.order/=-2).or.order>4)then
1005 : write(message, '(a,i0)' )&
1006 0 : & 'The only allowed values for order are 1, 2, -2 or 3, while it is found to be ',order
1007 0 : ABI_BUG(message)
1008 : end if
1009 :
1010 : !Determine quantities available in input arguments
1011 872659 : has_gradient=.false.;has_laplacian=.false.;has_kden=.false.
1012 872659 : if (usegradient==1) then
1013 356137 : if (.not.present(grho2_updn)) then
1014 0 : ABI_BUG('missing grho2_updn argument!')
1015 : end if
1016 356137 : if (nvxcgrho>0) then
1017 356137 : if (.not.present(vxcgrho)) then
1018 0 : ABI_BUG('missing vxcgrho argument!')
1019 : end if
1020 : has_gradient=.true.
1021 : end if
1022 516522 : else if (nvxcgrho>0) then
1023 0 : ABI_BUG('nvxcgrho>0 and usegradient=0!')
1024 : end if
1025 872659 : if (uselaplacian==1) then
1026 35279 : if (.not.present(lrho_updn)) then
1027 0 : ABI_BUG('missing lrho_updn argument!')
1028 : end if
1029 35279 : if (nvxclrho>0) then
1030 35279 : if (.not.present(vxclrho)) then
1031 0 : ABI_BUG('missing vxclrho argument!')
1032 : end if
1033 : has_laplacian=.true.
1034 : end if
1035 837380 : else if (nvxclrho>0) then
1036 0 : ABI_BUG('nvxclrho>0 and uselaplacian=0!')
1037 : end if
1038 872659 : if (usekden==1) then
1039 151688 : if (.not.present(tau_updn)) then
1040 0 : ABI_BUG('missing tau_updn argument!')
1041 : end if
1042 151688 : if (nvxctau>0) then
1043 151688 : if (.not.present(vxctau)) then
1044 0 : ABI_BUG('missing vxctau argument!')
1045 : end if
1046 : has_kden=.true.
1047 : end if
1048 720971 : else if (nvxctau>0) then
1049 0 : ABI_BUG('nvxctau>0 and usekden=0!')
1050 : end if
1051 872659 : if (abs(order)>=2) then
1052 145757 : if (.not.present(dvxc)) then
1053 0 : message='order>=2 needs argument dvxc!'
1054 0 : ABI_BUG(message)
1055 145757 : else if (ndvxc==0) then
1056 0 : message='order>=2 needs ndvxc>0!'
1057 0 : ABI_BUG(message)
1058 : end if
1059 : end if
1060 872659 : if (abs(order)>=3) then
1061 2544 : if (.not.present(d2vxc)) then
1062 0 : message='order>=3 needs argument d2vxc!'
1063 0 : ABI_BUG(message)
1064 2544 : else if (nd2vxc==0) then
1065 0 : message='order>=3 needs nd2vxc>0!'
1066 0 : ABI_BUG(message)
1067 : end if
1068 : end if
1069 :
1070 : !Determine quantities needed by XC functional
1071 872659 : if (present(xc_funcs)) then
1072 : call size_dvxc(ixc,order,nspden,usegradient=need_gradient,&
1073 : & uselaplacian=need_laplacian,usekden=need_kden,&
1074 : & nvxcgrho=need_nvxcgrho,nvxclrho=need_nvxclrho,&
1075 : & nvxctau=need_nvxctau,ndvxc=need_ndvxc,nd2vxc=need_nd2vxc,&
1076 6875 : & xc_funcs=xc_funcs)
1077 : else
1078 : call size_dvxc(ixc,order,nspden,usegradient=need_gradient,&
1079 : & uselaplacian=need_laplacian,usekden=need_kden,&
1080 : & nvxcgrho=need_nvxcgrho,nvxclrho=need_nvxclrho,&
1081 865784 : & nvxctau=need_nvxctau,ndvxc=need_ndvxc,nd2vxc=need_nd2vxc)
1082 : end if
1083 : if ((has_gradient.and.need_gradient>usegradient).or.&
1084 872659 : & (has_laplacian.and.need_laplacian>uselaplacian).or.&
1085 : & (has_kden.and.need_kden>usekden)) then
1086 : write(message, '(3a)' )&
1087 0 : & 'one of the arguments usegradient/uselaplacian/usesekden',ch10,&
1088 0 : & 'doesnt match the requirements of the XC functional!'
1089 0 : ABI_BUG(message)
1090 : end if
1091 : if ((has_gradient.and.need_nvxcgrho>nvxcgrho).or.&
1092 : & (has_laplacian.and.need_nvxclrho>nvxclrho).or.&
1093 : & (has_kden.and.need_nvxctau>nvxctau).or.&
1094 872659 : & need_ndvxc>ndvxc.or.need_nd2vxc>nd2vxc) then
1095 : write(message, '(3a)' )&
1096 0 : & 'one of the arguments nvxcgrho/nvxclrho/nvxctau/ndvxc/nd2vxc',ch10,&
1097 0 : & 'doesnt match the requirements of the XC functional!'
1098 0 : ABI_BUG(message)
1099 : end if
1100 : !Deactivate this test because, in case of mGGA, we can output derivatives involving
1101 : ! the density and its gradient. Derivatives involving tau or Laplacian will not be output.
1102 : !if (abs(order)>1.and.ixc<0.and.(need_laplacian==1.or.need_kden==1)) then
1103 : ! message='Derivatives of XC potential are not available in mGGA!'
1104 : ! ABI_BUG(message)
1105 : !end if
1106 :
1107 : !Check other optional arguments
1108 872659 : if (my_exexch/=0.and.usegradient==0) then
1109 0 : message='exexch argument only valid for GGA!'
1110 0 : ABI_BUG(message)
1111 : end if
1112 :
1113 : !Initialize exchange-correlation entropy energy density to zero
1114 1959819889 : tsxc(:)=zero
1115 :
1116 : ! =================================================
1117 : ! == Intermediate quantities computation ==
1118 : ! =================================================
1119 :
1120 : !If needed, compute rhotot and rs
1121 : if (ixc==1.or.ixc==2.or.ixc==3.or.ixc==4.or.ixc==5.or.ixc==6.or.&
1122 872659 : & ixc==21.or.ixc==22.or.ixc==50.or.ixc==51.or.ixc==60) then
1123 599205 : ABI_MALLOC(rhotot,(npts))
1124 399470 : ABI_MALLOC(rspts,(npts))
1125 199735 : if(nspden==1)then
1126 518577548 : rhotot(:)=two*rho_updn(:,1)
1127 : else
1128 79535808 : rhotot(:)=rho_updn(:,1)+rho_updn(:,2)
1129 : end if
1130 199735 : call invcb(rhotot,rspts,npts)
1131 598113356 : rspts(:)=rsfac*rspts(:)
1132 : end if
1133 :
1134 : !If needed, compute zeta
1135 872659 : if (ixc==1.or.ixc==21.or.ixc==22) then
1136 314790 : ABI_MALLOC(zeta,(npts))
1137 104930 : if(nspden==1)then
1138 269971389 : zeta(:)=zero
1139 : else
1140 80303537 : zeta(:)=two*rho_updn(:,1)/rhotot(:)-one
1141 : end if
1142 : end if
1143 :
1144 : ! =================================================
1145 : ! == XC energy, potentiel, ... computation ==
1146 : ! =================================================
1147 :
1148 : !>>>>> No exchange-correlation
1149 872659 : if (ixc==0.or.ixc==40) then
1150 0 : exc=zero ; vxcrho=zero
1151 0 : if (present(dvxc).and.ndvxc>0) dvxc(:,:)=zero
1152 0 : if (present(d2vxc).and.nd2vxc>0) d2vxc(:,:)=zero
1153 0 : if (present(vxcgrho).and.nvxcgrho>0) vxcgrho(:,:)=zero
1154 0 : if (present(vxclrho).and.nvxclrho>0) vxclrho(:,:)=zero
1155 0 : if (present(vxctau).and.nvxctau>0) vxctau(:,:)=zero
1156 :
1157 : !>>>>> New Teter fit (4/93) to Ceperley-Alder data, with spin-pol option
1158 : else if (ixc==1 .or. ixc==21 .or. ixc==22) then
1159 : ! new Teter fit (4/93) to Ceperley-Alder data, with spin-pol option
1160 104930 : if (order**2 <= 1) then
1161 99270 : call xcspol(exc,npts,nspden,order,rspts,vxcrho,zeta,ndvxc)
1162 : else
1163 5660 : call xcspol(exc,npts,nspden,order,rspts,vxcrho,zeta,ndvxc,dvxc)
1164 : end if
1165 :
1166 : !>>>>> Perdew-Zunger fit to Ceperly-Alder data (no spin-pol)
1167 : else if (ixc==2) then
1168 88406 : if (order**2 <= 1) then
1169 36632 : call xcpzca(exc,npts,order,rhotot,rspts,vxcrho(:,1))
1170 : else
1171 51774 : call xcpzca(exc,npts,order,rhotot,rspts,vxcrho(:,1),dvxc)
1172 : end if
1173 :
1174 : !>>>>> Teter fit (4/91) to Ceperley-Alder values (no spin-pol)
1175 : else if (ixc==3) then
1176 5044 : if (order**2 <= 1) then
1177 4881 : call xctetr(exc,npts,order,rhotot,rspts,vxcrho(:,1))
1178 163 : else if (order == 2) then
1179 137 : call xctetr(exc,npts,order,rhotot,rspts,vxcrho(:,1),dvxc=dvxc)
1180 26 : else if (order == 3) then
1181 26 : call xctetr(exc,npts,order,rhotot,rspts,vxcrho(:,1),d2vxc=d2vxc,dvxc=dvxc)
1182 : end if
1183 :
1184 : !>>>>> Wigner xc (no spin-pol)
1185 : else if (ixc==4) then
1186 396 : if (order**2 <= 1) then
1187 396 : call xcwign(exc,npts,order,rspts,vxcrho(:,1))
1188 : else
1189 0 : call xcwign(exc,npts,order,rspts,vxcrho(:,1),dvxc)
1190 : end if
1191 :
1192 : !>>>>> Hedin-Lundqvist xc (no spin-pol)
1193 : else if (ixc==5) then
1194 396 : if (order**2 <= 1) then
1195 396 : call xchelu(exc,npts,order,rspts,vxcrho(:,1))
1196 : else
1197 0 : call xchelu(exc,npts,order,rspts,vxcrho(:,1),dvxc)
1198 : end if
1199 :
1200 : !>>>>> X-alpha (no spin-pol)
1201 : else if (ixc==6) then
1202 198 : if (order**2 <= 1) then
1203 198 : call xcxalp(exc,npts,order,rspts,vxcrho(:,1))
1204 : else
1205 0 : call xcxalp(exc,npts,order,rspts,vxcrho(:,1),dvxc)
1206 : end if
1207 :
1208 : !>>>>> PBE and alternatives
1209 673289 : else if (((ixc>=7.and.ixc<=15).or.(ixc>=23.and.ixc<=24)).and.ixc/=10.and.ixc/=13) then
1210 : ! Perdew-Wang LSD is coded in Perdew-Burke-Ernzerhof GGA, with optpbe=1
1211 307889 : if(ixc==7)optpbe=1
1212 : ! x-only part of Perdew-Wang
1213 307889 : if(ixc==8)optpbe=-1
1214 : ! Exchange + RPA correlation from Perdew-Wang
1215 307889 : if(ixc==9)optpbe=3
1216 : ! Perdew-Burke-Ernzerhof GGA
1217 307889 : if(ixc==11)optpbe=2
1218 : ! x-only part of PBE
1219 307889 : if(ixc==12)optpbe=-2
1220 : ! C09x exchange of V. R. Cooper
1221 307889 : if(ixc==24)optpbe=-4
1222 : ! revPBE of Zhang and Yang
1223 307889 : if(ixc==14)optpbe=5
1224 : ! RPBE of Hammer, Hansen and Norskov
1225 307889 : if(ixc==15)optpbe=6
1226 : ! Wu and Cohen
1227 307889 : if(ixc==23)optpbe=7
1228 307889 : if (ixc >=7.and.ixc<=9) then
1229 230221 : if (order**2 <= 1) then
1230 194531 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc)
1231 35690 : else if (order /=3) then
1232 33174 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,dvxci=dvxc)
1233 : else if (order ==3) then
1234 2516 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,d2vxci=d2vxc,dvxci=dvxc)
1235 : end if
1236 77668 : else if ((ixc >= 11 .and. ixc <= 15) .or. (ixc>=23 .and. ixc<=24)) then
1237 77668 : if (order**2 <= 1) then
1238 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,&
1239 66672 : & dvxcdgr=vxcgrho,exexch=my_exexch,grho2_updn=grho2_updn)
1240 10996 : else if (order /=3) then
1241 10996 : if(ixc==12 .or. ixc==24)then
1242 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,&
1243 3 : & dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn)
1244 10993 : else if(ixc/=12 .or. ixc/=24) then
1245 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,&
1246 : & dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn)
1247 : end if
1248 : else if (order ==3) then
1249 0 : if(ixc==12 .or. ixc==24)then
1250 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,&
1251 0 : & d2vxci=d2vxc,dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn)
1252 0 : else if(ixc/=12 .or. ixc/=24) then
1253 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,&
1254 : & d2vxci=d2vxc,dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn)
1255 : end if
1256 : end if
1257 : end if
1258 :
1259 : !>>>>> RPA correlation from Perdew-Wang
1260 : else if (ixc==10) then
1261 0 : if (order**2 <= 1) then
1262 0 : ABI_MALLOC(exci_rpa,(npts))
1263 0 : ABI_MALLOC(vxci_rpa,(npts,2))
1264 0 : optpbe=3
1265 0 : call xcpbe(exci_rpa,npts,nspden,optpbe,order,rho_updn,vxci_rpa,ndvxc,nd2vxc)
1266 0 : optpbe=1
1267 0 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc)
1268 0 : exc(:)=exc(:)-exci_rpa(:)
1269 : ! PMA: second index of vxcrho is nspden while that of rpa is 2 they can mismatch
1270 0 : vxcrho(:,1:min(nspden,2))=vxcrho(:,1:min(nspden,2))-vxci_rpa(:,1:min(nspden,2))
1271 0 : ABI_FREE(exci_rpa)
1272 0 : ABI_FREE(vxci_rpa)
1273 0 : else if (order /=3) then
1274 0 : ABI_MALLOC(exci_rpa,(npts))
1275 0 : ABI_MALLOC(vxci_rpa,(npts,2))
1276 0 : optpbe=3
1277 0 : call xcpbe(exci_rpa,npts,nspden,optpbe,order,rho_updn,vxci_rpa,ndvxc,nd2vxc,dvxci=dvxc)
1278 0 : optpbe=1
1279 0 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,dvxci=dvxc)
1280 0 : exc(:)=exc(:)-exci_rpa(:)
1281 0 : vxcrho(:,:)=vxcrho(:,:)-vxci_rpa(:,:)
1282 0 : ABI_FREE(exci_rpa)
1283 0 : ABI_FREE(vxci_rpa)
1284 : else if (order ==3) then
1285 0 : ABI_MALLOC(exci_rpa,(npts))
1286 0 : ABI_MALLOC(vxci_rpa,(npts,2))
1287 0 : optpbe=3
1288 : call xcpbe(exci_rpa,npts,nspden,optpbe,order,rho_updn,vxci_rpa,ndvxc,nd2vxc,&
1289 0 : & d2vxci=d2vxc,dvxci=dvxc)
1290 0 : optpbe=1
1291 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,&
1292 0 : & d2vxci=d2vxc,dvxci=dvxc)
1293 0 : exc(:)=exc(:)-exci_rpa(:)
1294 0 : vxcrho(:,:)=vxcrho(:,:)-vxci_rpa(:,:)
1295 0 : ABI_FREE(exci_rpa)
1296 0 : ABI_FREE(vxci_rpa)
1297 : end if
1298 :
1299 : !>>>>> LDA xc energy like ixc==7, and Leeuwen-Baerends GGA xc potential
1300 : else if(ixc==13) then
1301 0 : if (order**2 <= 1) then
1302 0 : optpbe=1
1303 0 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc)
1304 0 : call xclb(grho2_updn,npts,nspden,rho_updn,vxcrho)
1305 0 : else if (order /=3) then
1306 0 : optpbe=1
1307 0 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,dvxci=dvxc)
1308 0 : call xclb(grho2_updn,npts,nspden,rho_updn,vxcrho)
1309 : else if (order ==3) then
1310 0 : optpbe=1
1311 0 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,d2vxci=d2vxc,dvxci=dvxc)
1312 0 : call xclb(grho2_updn,npts,nspden,rho_updn,vxcrho)
1313 : end if
1314 :
1315 : !>>>>> HTCH93, HTCH120, HTCH107, HTCH147
1316 : else if(ixc==16 .or. ixc==17 .or. ixc==26 .or. ixc==27) then
1317 1265 : call xchcth(vxcgrho,exc,grho2_updn,ixc,npts,nspden,order,rho_updn,vxcrho)
1318 :
1319 : !>>>>> Only for test purpose (test various part of MGGA implementation)
1320 : else if(ixc==31 .or. ixc==32 .or. ixc==33 .or. ixc==34 .or. ixc==35) then
1321 13437411 : exc(:)=zero ; vxcrho(:,:)=zero
1322 16048836 : if (present(vxcgrho).and.nvxcgrho>0) vxcgrho(:,:)=zero
1323 1260789 : if (present(vxclrho).and.nvxclrho>0) vxclrho(:,:)=zero
1324 3904315 : if (present(vxctau).and.nvxctau>0) vxctau(:,:)=zero
1325 753 : if (present(dvxc).and.ndvxc>0) dvxc(:,:)=zero
1326 753 : if (present(d2vxc).and.nd2vxc>0) d2vxc(:,:)=zero
1327 :
1328 : !>>>>> Perdew-Wang LSD is coded in Perdew-Burke-Ernzerhof GGA, with optpbe=1
1329 753 : optpbe=1
1330 32 : select case(ixc)
1331 : case (31)
1332 32 : alpha=1.00d0-(1.00d0/1.01d0)
1333 : ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc)
1334 32 : call xcpbe(exc,npts,nspden,optpbe,1,rho_updn,vxcrho,0,0)
1335 32 : if (nspden==1) then
1336 746512 : exc(:)=exc(:)+alpha*tau_updn(:,1)/rho_updn(:,1)
1337 : else
1338 : ! It should be : exc_tot= exc_spin up + exc_spin down = 2*exc_spin up
1339 : ! but this applies to tau and rho (so it cancels)
1340 48 : do ispden=1,nspden
1341 1493040 : exc(:)=exc(:)+alpha*tau_updn(:,ispden)/(rho_updn(:,1)+rho_updn(:,2))
1342 : end do
1343 : end if
1344 2239568 : vxctau(:,:)=alpha
1345 : case (32)
1346 216 : alpha=0.01d0
1347 : ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc)
1348 216 : call xcpbe(exc,npts,nspden,optpbe,1,rho_updn,vxcrho,0,0)
1349 216 : if (nspden==1) then
1350 420012 : exc(:)=exc(:)+2.0d0*alpha*lrho_updn(:,1)
1351 420012 : vxcrho(:,1) =vxcrho(:,1)+2.0d0*alpha*lrho_updn(:,1)
1352 420012 : vxclrho(:,1)=alpha*2.0d0*rho_updn(:,1)
1353 : else
1354 324 : do ispden=1,nspden
1355 840024 : exc(:)=exc(:)+alpha*lrho_updn(:,ispden)
1356 840024 : vxcrho(:,ispden) =vxcrho(:,ispden)+alpha*(lrho_updn(:,1)+lrho_updn(:,2))
1357 840132 : vxclrho(:,ispden)=alpha*(rho_updn(:,1)+rho_updn(:,2))
1358 : end do
1359 : end if
1360 : case (33)
1361 488 : alpha=-0.010d0
1362 : ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc)
1363 488 : call xcpbe(exc,npts,nspden,optpbe,1,rho_updn,vxcrho,0,0)
1364 488 : if (nspden==1) then
1365 : ! it should be : exc_tot= exc_spin up + exc_spin down = 2*exc_spin up but this applies to grho2 and rho
1366 : ! (for grho2 it is a factor 4 to have total energy and for rho it is just a factor 2. So we end with factor 2 only)
1367 932140 : exc(:)=exc(:)+alpha*2.0d0*grho2_updn(:,1)/rho_updn(:,1)
1368 236 : if(nvxcgrho==2)vxcgrho(:,1:2)=2.0d0*alpha
1369 932140 : if(nvxcgrho==3)vxcgrho(:,3)=2.0d0*alpha
1370 : else
1371 996156 : exc(:)=exc(:)+alpha*grho2_updn(:,3)/(rho_updn(:,1)+rho_updn(:,2))
1372 252 : if(nvxcgrho==2)vxcgrho(:,1:2)=2.0d0*alpha
1373 996156 : if(nvxcgrho==3)vxcgrho(:,3)=2.0d0*alpha
1374 : end if
1375 : case (34)
1376 17 : alpha=-0.010d0
1377 : ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc)
1378 17 : call xcpbe(exc,npts,nspden,optpbe,1,rho_updn,vxcrho,0,0)
1379 17 : if (nspden==1) then
1380 512008 : exc(:)=exc(:)+16.0d0*alpha*tau_updn(:,1)
1381 512008 : vxcrho(:,1)=vxcrho(:,1)+16.0d0*alpha*tau_updn(:,1)
1382 512008 : vxctau(:,1)=16.0d0*alpha*rho_updn(:,1)
1383 : else
1384 27 : do ispden=1,nspden
1385 1152018 : exc(:)=exc(:)+8.0d0*alpha*tau_updn(:,ispden)
1386 1152018 : vxcrho(:,ispden)=vxcrho(:,ispden)+8.0d0*alpha*(tau_updn(:,1)+tau_updn(:,2))
1387 1152027 : vxctau(:,ispden)=8.0d0*alpha*(rho_updn(:,1)+rho_updn(:,2))
1388 : end do
1389 : end if
1390 : case (35)
1391 0 : alpha=0.01d0 ; beta=1.00d0-(1.00d0/1.01d0)
1392 : ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc)
1393 0 : call xcpbe(exc,npts,nspden,optpbe,1,rho_updn,vxcrho,0,0)
1394 0 : if (nspden==1) then
1395 0 : exc(:)=exc(:)+2.0d0*alpha*lrho_updn(:,1)+beta*tau_updn(:,1)/rho_updn(:,1)
1396 0 : vxcrho(:,1) =vxcrho(:,1)+2.0d0*alpha*lrho_updn(:,1)
1397 0 : vxclrho(:,1)=alpha*2.0d0*rho_updn(:,1)
1398 : else
1399 0 : do ispden=1,nspden
1400 : exc(:)=exc(:)+alpha*lrho_updn(:,ispden) &
1401 0 : & +beta*tau_updn(:,ispden)/(rho_updn(:,1)+rho_updn(:,2))
1402 0 : vxcrho(:,ispden) =vxcrho(:,ispden)+alpha*(lrho_updn(:,1)+lrho_updn(:,2))
1403 0 : vxclrho(:,ispden)=alpha*(rho_updn(:,1)+rho_updn(:,2))
1404 : end do
1405 : end if
1406 753 : vxctau(:,:)=beta
1407 : end select
1408 :
1409 : !>>>>> Hybrid PBE0 (1/4 and 1/3)
1410 : else if(ixc>=41.and.ixc<=42) then
1411 : ! Requires to evaluate exchange-correlation with PBE (optpbe=2)
1412 : ! minus hyb_mixing*exchange with PBE (optpbe=-2)
1413 4610 : ndvxc_x=8
1414 13830 : ABI_MALLOC(exc_x,(npts))
1415 18440 : ABI_MALLOC(vxcrho_x,(npts,nspden))
1416 18440 : ABI_MALLOC(vxcgrho_x,(npts,nvxcgrho))
1417 64164326 : exc_x=zero;vxcrho_x=zero;vxcgrho_x=zero
1418 4610 : if (order**2 <= 1) then
1419 2734 : optpbe=2 !PBE exchange correlation
1420 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,&
1421 2734 : & dvxcdgr=vxcgrho,exexch=my_exexch,grho2_updn=grho2_updn)
1422 2734 : optpbe=-2 !PBE exchange-only
1423 : call xcpbe(exc_x,npts,nspden,optpbe,order,rho_updn,vxcrho_x,ndvxc,nd2vxc,&
1424 2734 : & dvxcdgr=vxcgrho_x,exexch=my_exexch,grho2_updn=grho2_updn)
1425 4749610 : exc=exc-exc_x*my_hyb_mixing
1426 6046276 : vxcrho=vxcrho-vxcrho_x*my_hyb_mixing
1427 14251564 : vxcgrho=vxcgrho-vxcgrho_x*my_hyb_mixing
1428 1876 : else if (order /=3) then
1429 7504 : ABI_MALLOC(dvxc_x,(npts,ndvxc_x))
1430 1876 : optpbe=2 !PBE exchange correlation
1431 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,&
1432 1876 : dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn)
1433 1876 : optpbe=-2 !PBE exchange-only
1434 : call xcpbe(exc_x,npts,nspden,optpbe,order,rho_updn,vxcrho_x,ndvxc_x,nd2vxc,&
1435 1876 : & dvxcdgr=vxcgrho_x,dvxci=dvxc_x,grho2_updn=grho2_updn)
1436 7211952 : exc=exc-exc_x*my_hyb_mixing
1437 10276412 : vxcrho=vxcrho-vxcrho_x*my_hyb_mixing
1438 21637732 : vxcgrho=vxcgrho-vxcgrho_x*my_hyb_mixing
1439 57697492 : dvxc(:,1:ndvxc_x)=dvxc(:,1:ndvxc_x)-dvxc_x(:,1:ndvxc_x)*my_hyb_mixing
1440 1876 : ABI_FREE(dvxc_x)
1441 : else if (order ==3) then
1442 : ! The size of exchange-correlation with PBE (optpbe=2)
1443 : ! is the one which defines the size for ndvxc.
1444 0 : ABI_MALLOC(dvxc_x,(npts,ndvxc_x))
1445 0 : ABI_MALLOC(d2vxc_x,(npts,nd2vxc))
1446 0 : optpbe=2 !PBE exchange correlation
1447 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,&
1448 0 : & d2vxci=d2vxc,dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn)
1449 0 : optpbe=-2 !PBE exchange-only
1450 : call xcpbe(exc_x,npts,nspden,optpbe,order,rho_updn,vxcrho_x,ndvxc_x,nd2vxc,&
1451 0 : & d2vxci=d2vxc_x,dvxcdgr=vxcgrho_x,dvxci=dvxc_x,grho2_updn=grho2_updn)
1452 0 : exc=exc-exc_x*my_hyb_mixing
1453 0 : vxcrho=vxcrho-vxcrho_x*my_hyb_mixing
1454 0 : vxcgrho=vxcgrho-vxcgrho_x*my_hyb_mixing
1455 0 : d2vxc=d2vxc-d2vxc_x*my_hyb_mixing
1456 0 : dvxc(:,1:ndvxc_x)=dvxc(:,1:ndvxc_x)-dvxc_x(:,1:ndvxc_x)*my_hyb_mixing
1457 0 : ABI_FREE(dvxc_x)
1458 0 : ABI_FREE(d2vxc_x)
1459 : end if
1460 4610 : ABI_FREE(exc_x)
1461 4610 : ABI_FREE(vxcrho_x)
1462 4610 : ABI_FREE(vxcgrho_x)
1463 :
1464 : !>>>>> Finite-temperature XC functionals.
1465 : !>>>>> exc is the xc free energy density.
1466 : !>>>>> XC entropy energy density 'tsxc' is needed to retrieve
1467 : !>>>>> the proper internal energy E_xc
1468 :
1469 : !>>>>> Ichimaru-Iyetomi-Tanaka, TLDA (no spin-pol) (IIT)
1470 : else if (ixc==50) then
1471 252 : if (order**2 <= 1) then
1472 252 : call xciit(exc,tsxc,npts,order,rspts,el_temp,vxcrho(:,1))
1473 : else
1474 0 : call xciit(exc,tsxc,npts,order,rspts,el_temp,vxcrho(:,1),dvxc=dvxc)
1475 : end if
1476 :
1477 : !>>>>> Karasiev-Sjostrom-Dufty-Trickey, TLDA (no spin-pol) (KSDT)
1478 : else if (ixc==51) then
1479 56 : if (order**2 <= 1) then
1480 49 : call xcksdt(exc,tsxc,npts,order,rhotot,rspts,el_temp,vxcrho(:,1))
1481 : else
1482 7 : call xcksdt(exc,tsxc,npts,order,rhotot,rspts,el_temp,vxcrho(:,1),dvxc=dvxc)
1483 : end if
1484 :
1485 : !>>>>> Karasiev-Dufty-Trickey, TGGA (no spin-pol) (KDT16)
1486 : else if(ixc==60) then
1487 57 : if(nvxcgrho /= 3 )then
1488 : write(message, '(3a,i0,a,i0)')&
1489 0 : & 'Wrong value of nvxcgrho:',ch10,&
1490 0 : & 'ixc=',ixc,'ndvxcdgr=',nvxcgrho
1491 0 : ABI_BUG(message)
1492 : end if
1493 57 : if (order**2 <= 1) then
1494 : call xckdt16(vxcgrho,exc,tsxc,grho2_updn,ixc,npts,nspden,order,&
1495 50 : & rhotot,rspts,el_temp,vxcrho)
1496 : else
1497 : call xckdt16(vxcgrho,exc,tsxc,grho2_updn,ixc,npts,nspden,order,&
1498 7 : & rhotot,rspts,el_temp,vxcrho,dvxci=dvxc)
1499 : end if
1500 :
1501 : !>>>>> GGA counterpart of the B3LYP functional
1502 : else if(ixc==1402000) then
1503 : ! Requires to evaluate exchange-correlation
1504 : ! with 5/4 B3LYP - 1/4 B3LYPc, where
1505 : ! B3LYPc = (0.19 Ec VWN3 + 0.81 Ec LYP)
1506 :
1507 : ! First evaluate B3LYP.
1508 18 : if(present(xc_funcs))then
1509 0 : if (abs(order)==1) then
1510 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,&
1511 0 : & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,xc_functionals=xc_funcs)
1512 0 : else if (abs(order)==2) then
1513 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,&
1514 0 : & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,xc_functionals=xc_funcs)
1515 : else
1516 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,&
1517 0 : & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,d2vxc=d2vxc,xc_functionals=xc_funcs)
1518 : end if
1519 : else
1520 18 : if (abs(order)==1) then
1521 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,&
1522 18 : & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho)
1523 0 : else if (abs(order)==2) then
1524 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,&
1525 0 : & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc)
1526 : else
1527 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,&
1528 0 : & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,d2vxc=d2vxc)
1529 : end if
1530 : end if
1531 :
1532 : ! Then renormalize B3LYP and subtract VWN3 contribution
1533 54 : ABI_MALLOC(exc_c,(npts))
1534 72 : ABI_MALLOC(vxcrho_c,(npts,nspden))
1535 18 : if(order**2>1)then
1536 0 : ABI_MALLOC(dvxc_c,(npts,ndvxc))
1537 : end if
1538 18 : if(order**2>4)then
1539 0 : ABI_MALLOC(d2vxc_c,(npts,nd2vxc))
1540 : end if
1541 62244 : exc_c=zero;vxcrho_c=zero
1542 18 : call libxc_functionals_init(-30,nspden,xc_functionals=xc_funcs_vwn3)
1543 18 : if (order**2 <= 1) then
1544 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,&
1545 18 : & vxcrho_c,xc_functionals=xc_funcs_vwn3)
1546 0 : elseif (order**2 <= 4) then
1547 0 : dvxc_c=zero
1548 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,&
1549 0 : & vxcrho_c,dvxc=dvxc_c,xc_functionals=xc_funcs_vwn3)
1550 : else
1551 0 : dvxc_c=zero
1552 0 : d2vxc_c=zero
1553 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,&
1554 0 : & vxcrho_c,dvxc=dvxc_c,d2vxc=d2vxc,xc_functionals=xc_funcs_vwn3)
1555 : end if
1556 31122 : exc=1.25d0*exc-quarter*0.19d0*exc_c
1557 31140 : vxcrho=1.25d0*vxcrho-quarter*0.19d0*vxcrho_c
1558 18 : if(order**2>1)dvxc=1.25d0*dvxc-quarter*0.19d0*dvxc_c
1559 18 : if(order**2>4)d2vxc=1.25d0*d2vxc-quarter*0.19d0*d2vxc_c
1560 18 : call libxc_functionals_end(xc_functionals=xc_funcs_vwn3)
1561 :
1562 : ! Then subtract LYP contribution
1563 18 : call libxc_functionals_init(-131,nspden,xc_functionals=xc_funcs_lyp)
1564 18 : if (order**2 <= 1) then
1565 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,&
1566 18 : & vxcrho_c,grho2=grho2_updn,vxcgr=vxcgrho,xc_functionals=xc_funcs_lyp)
1567 0 : elseif (order**2 <= 4) then
1568 0 : dvxc_c=zero
1569 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,&
1570 0 : & vxcrho_c,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc_c,xc_functionals=xc_funcs_lyp)
1571 : else
1572 0 : dvxc_c=zero
1573 0 : d2vxc_c=zero
1574 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,&
1575 0 : & vxcrho_c,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc_c,d2vxc=d2vxc,xc_functionals=xc_funcs_lyp)
1576 : end if
1577 31122 : exc=exc-quarter*0.81d0*exc_c
1578 31140 : vxcrho=vxcrho-quarter*0.81d0*vxcrho_c
1579 18 : if(order**2>1)dvxc=dvxc-quarter*0.81d0*dvxc_c
1580 18 : if(order**2>4)d2vxc=d2vxc-quarter*0.81d0*d2vxc_c
1581 18 : call libxc_functionals_end(xc_functionals=xc_funcs_lyp)
1582 :
1583 18 : ABI_FREE(exc_c)
1584 18 : ABI_FREE(vxcrho_c)
1585 18 : if(allocated(dvxc_c))then
1586 0 : ABI_FREE(dvxc_c)
1587 : end if
1588 18 : if(allocated(d2vxc_c))then
1589 0 : ABI_FREE(d2vxc_c)
1590 : end if
1591 :
1592 : !>>>>> All libXC functionals
1593 358389 : else if( ixc<0 ) then
1594 :
1595 : ! ===== meta-GGA =====
1596 358389 : if (need_laplacian==1.or.need_kden==1) then
1597 151639 : if (need_laplacian==1.and.need_kden==1) then
1598 35063 : if (abs(order)<=1) then
1599 35063 : if (present(xc_funcs)) then
1600 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1601 : & grho2=grho2_updn,vxcgr=vxcgrho,&
1602 : & lrho=lrho_updn,vxclrho=vxclrho,&
1603 : & tau=tau_updn,vxctau=vxctau,&
1604 0 : & xc_functionals=xc_funcs)
1605 : else
1606 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1607 : & grho2=grho2_updn,vxcgr=vxcgrho,&
1608 : & lrho=lrho_updn,vxclrho=vxclrho,&
1609 35063 : & tau=tau_updn,vxctau=vxctau)
1610 : end if
1611 : else
1612 0 : if (present(xc_funcs)) then
1613 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1614 : & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,&
1615 : & lrho=lrho_updn,vxclrho=vxclrho,&
1616 : & tau=tau_updn,vxctau=vxctau,&
1617 0 : & xc_functionals=xc_funcs)
1618 : else
1619 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1620 : & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,&
1621 : & lrho=lrho_updn,vxclrho=vxclrho,&
1622 0 : & tau=tau_updn,vxctau=vxctau)
1623 : end if
1624 : end if
1625 116576 : else if (need_laplacian==1) then
1626 0 : if (abs(order)<=1) then
1627 0 : if (present(xc_funcs)) then
1628 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1629 : & grho2=grho2_updn,vxcgr=vxcgrho,&
1630 : & lrho=lrho_updn,vxclrho=vxclrho,&
1631 0 : & xc_functionals=xc_funcs)
1632 : else
1633 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1634 : & grho2=grho2_updn,vxcgr=vxcgrho,&
1635 0 : & lrho=lrho_updn,vxclrho=vxclrho)
1636 : end if
1637 : else
1638 0 : if (present(xc_funcs)) then
1639 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1640 : & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,&
1641 : & lrho=lrho_updn,vxclrho=vxclrho,&
1642 0 : & xc_functionals=xc_funcs)
1643 : else
1644 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1645 : & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,&
1646 0 : & lrho=lrho_updn,vxclrho=vxclrho)
1647 : end if
1648 : end if
1649 116576 : else if (need_kden==1) then
1650 116576 : if (abs(order)<=1) then
1651 116100 : if (present(xc_funcs)) then
1652 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1653 : & grho2=grho2_updn,vxcgr=vxcgrho,&
1654 : & tau=tau_updn,vxctau=vxctau,&
1655 0 : & xc_functionals=xc_funcs)
1656 : else
1657 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1658 : & grho2=grho2_updn,vxcgr=vxcgrho,&
1659 116100 : & tau=tau_updn,vxctau=vxctau)
1660 : end if
1661 : else
1662 476 : if (present(xc_funcs)) then
1663 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1664 : & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,&
1665 : & tau=tau_updn,vxctau=vxctau,&
1666 0 : & xc_functionals=xc_funcs)
1667 : else
1668 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1669 : & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,&
1670 476 : & tau=tau_updn,vxctau=vxctau)
1671 : end if
1672 : end if
1673 : end if
1674 : !Some meta-GGAs can only be used with a LDA correlation (see doc)
1675 151639 : ixc1=(-ixc)/1000;ixc2=(-ixc)-ixc1*1000
1676 : if (ixc1==206 .or. ixc1==207 .or. ixc1==208 .or. ixc1==209 .or. &
1677 151639 : & ixc2==206 .or. ixc2==207 .or. ixc2==208 .or. ixc2==209 )then
1678 162077498 : if (present(vxcgrho)) vxcgrho(:,:)=zero
1679 54049208 : if (present(vxclrho)) vxclrho(:,:)=zero
1680 54049208 : if (present(vxctau)) vxctau(:,:)=zero
1681 35063 : if (present(dvxc)) dvxc(:,:)=zero
1682 : end if
1683 :
1684 : ! ===== GGA =====
1685 206750 : else if (need_gradient==1) then
1686 119938 : if (abs(order)<=1) then
1687 95026 : if (present(xc_funcs)) then
1688 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1689 : & grho2=grho2_updn,vxcgr=vxcgrho,&
1690 1024 : & xc_functionals=xc_funcs)
1691 : else
1692 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1693 94002 : & grho2=grho2_updn,vxcgr=vxcgrho)
1694 : end if
1695 : else
1696 24912 : if (present(xc_funcs)) then
1697 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1698 : & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,&
1699 0 : & xc_functionals=xc_funcs)
1700 : else
1701 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1702 24912 : & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc)
1703 : end if
1704 : end if
1705 :
1706 : ! ===== LDA =====
1707 : else
1708 86812 : if (abs(order)<=1) then
1709 72616 : if (present(xc_funcs)) then
1710 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1711 0 : & xc_functionals=xc_funcs)
1712 : else
1713 72616 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho)
1714 : end if
1715 14196 : else if (abs(order)<=2) then
1716 14194 : if (present(xc_funcs)) then
1717 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1718 0 : & dvxc=dvxc,xc_functionals=xc_funcs)
1719 : else
1720 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1721 14194 : & dvxc=dvxc)
1722 : end if
1723 2 : else if (abs(order)<=3) then
1724 2 : if (present(xc_funcs)) then
1725 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1726 0 : & dvxc=dvxc,d2vxc=d2vxc,xc_functionals=xc_funcs)
1727 : else
1728 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,&
1729 2 : & dvxc=dvxc,d2vxc=d2vxc)
1730 : end if
1731 : end if
1732 :
1733 : end if ! mGGA, GGA, LDA
1734 : end if ! libXC
1735 :
1736 : ! =================================================
1737 : ! == Finalization ==
1738 : ! =================================================
1739 : !Deallocate arrays
1740 872659 : if(allocated(rhotot)) then
1741 199735 : ABI_FREE(rhotot)
1742 : end if
1743 872659 : if(allocated(rspts)) then
1744 199735 : ABI_FREE(rspts)
1745 : end if
1746 872659 : if(allocated(zeta)) then
1747 104930 : ABI_FREE(zeta)
1748 : end if
1749 :
1750 4684238 : end subroutine drivexc
1751 : !!***
1752 :
1753 : end module m_drivexc
1754 : !!***
|
