Branch data 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 : 5340 : 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 : 5340 : message =''
78 : 5340 : citation =''
79 : :
80 : : !normal case (not libxc)
81 [ + + ]: 5340 : 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 : 2181 : message = 'LDA: new Teter (4/93) with spin-polarized option - ixc=1'
90 : 2181 : 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 : 673 : message = 'LDA: Perdew-Wang 92 LSD fit to Ceperley-Alder data - ixc=7'
108 : 673 : 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 : 618 : message = 'GGA: Perdew-Burke-Ernzerhof functional - ixc=11'
121 : 618 : 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 [ + + + + : 4148 : ABI_WARNING(message)
+ + + + +
- - + + -
+ + + + +
+ + + + +
+ + + + +
- + + + +
+ + - ]
201 : : end select
202 : :
203 : : message = " Exchange-correlation functional for the present dataset will be:" // ch10 &
204 : 4148 : & // " " // trim(message)
205 : :
206 : 4148 : l_citation=len_trim(citation)
207 : 4148 : citation = " Citation for XC functional:" // ch10 // " " // trim(citation)
208 : :
209 : 4148 : call wrtout(ab_out,message,'COLL')
210 : 4148 : call wrtout(std_out,message,'COLL')
211 : :
212 [ + + ]: 4148 : if(l_citation/=0)then
213 : 3728 : call wrtout(ab_out,citation,'COLL')
214 : 3728 : call wrtout(std_out,citation,'COLL')
215 : : end if
216 : :
217 : 4148 : message =' '
218 : 4148 : call wrtout(ab_out,message,'COLL')
219 : 4148 : call wrtout(std_out,message,'COLL')
220 : :
221 : : end if ! end libxc if
222 : :
223 : 5340 : 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 : 1268150 : 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 : 1268150 : xc_need_kden=.false.
247 : :
248 [ + + ]: 1268150 : if (ixc>=0) then
249 [ + + - + ]: 752941 : xc_need_kden=(ixc==31.or.ixc==34.or.ixc==35)
250 : : else
251 [ + + ]: 515209 : if(present(xc_funcs)) xc_need_kden=libxc_functionals_needs_tau(xc_functionals=xc_funcs)
252 : 513664 : if(.not.present(xc_funcs)) xc_need_kden=libxc_functionals_needs_tau()
253 : : end if
254 : :
255 : 1268150 : 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 : 9138 : 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 : 9138 : has_kxc=.false.
281 : :
282 [ + + ]: 9138 : if (ixc>=0) then
283 [ + + + + ]: 7322 : has_kxc=(ixc/=16.and.ixc/=17.and.ixc/=26.and.ixc/=27)
284 [ + + ]: 1816 : 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 [ + + ]: 1746 : if (present(xc_funcs)) then
288 : 1151 : has_kxc=libxc_functionals_has_kxc(xc_funcs)
289 : : else
290 : 595 : has_kxc=libxc_functionals_has_kxc()
291 : : end if
292 : : end if
293 : :
294 : 9138 : 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 : 2124 : 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 : 2124 : has_k3xc=.false.
321 : :
322 [ + + ]: 2124 : if (ixc>=0) then
323 : : has_k3xc=(ixc==0.or.ixc==3.or.(ixc>=7.and.ixc<=15).or. &
324 [ + + ]: 1520 : & ixc==23.or.ixc==24.or.ixc==41.or.ixc==42.or.ixc==1402000)
325 [ + + ]: 604 : 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 [ - + ]: 595 : if (present(xc_funcs)) then
329 : 0 : has_k3xc=libxc_functionals_has_k3xc(xc_funcs)
330 : : else
331 : 595 : has_k3xc=libxc_functionals_has_k3xc()
332 : : end if
333 : : end if
334 : :
335 : 2124 : 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 : 2124 : 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 [ + + ]: 2124 : check_k3xc_=.false. ; if (present(check_k3xc)) check_k3xc_=check_k3xc
366 : :
367 : 2124 : kxc_available=has_kxc(ixc)
368 : 2124 : k3xc_available=has_k3xc(ixc)
369 : :
370 [ + + ]: 2124 : if (ixc>=0) then
371 [ - + ]: 1520 : 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 [ - + ]: 1520 : 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 [ - + ]: 604 : 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 [ - + ]: 604 : 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 [ - + ]: 2124 : 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 [ - + ]: 2124 : 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 : 2124 : 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 : 1260422 : 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 [ + + ]: 1260422 : my_add_tfw=.false.;if (present(add_tfw)) my_add_tfw=add_tfw
476 : 1260422 : libxc_isgga=.false. ; libxc_ismgga=.false. ; libxc_ishybrid=.false.
477 [ + + ]: 1260422 : if(ixc<0)then
478 [ + + ]: 513720 : 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 : 512175 : libxc_has_kxc=libxc_functionals_has_kxc()
486 : 512175 : libxc_has_k3xc=libxc_functionals_has_k3xc()
487 : 512175 : libxc_isgga=libxc_functionals_isgga()
488 : 512175 : libxc_ismgga=libxc_functionals_ismgga()
489 : 512175 : 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 [ + + ]: 1260422 : & (ixc==41.or.ixc==42).or.ixc==1402000.or.ixc==60)
497 [ + + + + : 1260422 : if (ixc<0.and.(libxc_isgga.or.libxc_ismgga.or.libxc_ishybrid)) need_gradient=.true.
+ - ]
498 [ + + ]: 1260422 : if (my_add_tfw) need_gradient=.true.
499 [ + - + + ]: 2042907 : if (present(usegradient)) usegradient=merge(1,0,need_gradient)
500 : :
501 : : !Do we use the laplacian?
502 : 1260422 : need_laplacian=(ixc==32.or.ixc==35)
503 [ + + ]: 1260422 : if (ixc<0) then
504 [ + + ]: 513720 : if(present(xc_funcs)) need_laplacian=libxc_functionals_needs_laplacian(xc_functionals=xc_funcs)
505 : 512175 : if(.not.present(xc_funcs)) need_laplacian=libxc_functionals_needs_laplacian()
506 : : end if
507 [ + - + + ]: 2485069 : if (present(uselaplacian)) uselaplacian=merge(1,0,need_laplacian)
508 : :
509 : : !Do we use the kinetic energy density?
510 [ + + ]: 1260422 : if(present(xc_funcs)) need_kden=xc_need_kden(ixc,xc_funcs)
511 : 1253026 : if(.not.present(xc_funcs)) need_kden=xc_need_kden(ixc)
512 [ + + ]: 1260422 : usekden=merge(1,0,need_kden)
513 : :
514 : : !First derivative(s) of XC functional wrt gradient of density
515 [ + + ]: 1260422 : if (present(nvxcgrho)) then
516 : 1203258 : nvxcgrho=0
517 [ + - ]: 1203258 : if (abs(order)>=1) then
518 [ + + ]: 1203258 : if (need_gradient) nvxcgrho=3
519 [ - + ]: 1203258 : if (ixc==13) nvxcgrho=0
520 [ + + + + ]: 1203258 : 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 [ + + ]: 1260422 : if (present(nvxclrho)) then
526 : 1203258 : nvxclrho=0
527 [ + - ]: 1203258 : if (abs(order)>=1) then
528 [ + + ]: 1203258 : 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 [ + + ]: 1260422 : if (present(nvxctau)) then
534 : 1203258 : nvxctau=0
535 [ + - ]: 1203258 : if (abs(order)>=1) then
536 [ + + ]: 1203258 : if (need_kden) nvxctau=min(nspden,2)
537 : : end if
538 : : end if
539 : :
540 : : !Second derivative(s) of XC functional wrt density
541 [ + + ]: 1260422 : if (present(ndvxc)) then
542 : 1203258 : ndvxc=0
543 [ + + ]: 1203258 : 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 [ + + + - : 164509 : & ixc==21.or.ixc==22) then
+ - + - +
- + + +
+ ]
546 : 48547 : ndvxc=min(nspden,2)+1
547 [ + + + + ]: 115962 : 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 [ + + ]: 62144 : 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 [ + + + + : 62140 : & ixc==23.or.ixc==41.or.ixc==42.or.ixc==1402000) then
+ + + + +
- - + ]
554 : 13949 : ndvxc=15
555 [ + + ]: 48191 : else if (ixc<0) then
556 [ + + ]: 48161 : if (libxc_has_kxc.or.ixc==-406.or.ixc==-427.or.ixc==-428.or.ixc==-456) then
557 [ - + ]: 48161 : ndvxc=2*min(nspden,2)+1 ; if (order==-2) ndvxc=2
558 [ + + ]: 48161 : 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 [ + + ]: 1260422 : if (present(nd2vxc)) then
567 : 1203258 : nd2vxc=0
568 [ + + ]: 1203258 : 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 : 1260422 : 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 : 203663 : 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 [ + + ]: 814652 : do idir=1,3
646 : :
647 [ + + ]: 814652 : if(nspden==1)then
648 : : !$OMP PARALLEL DO PRIVATE(ifft) SHARED(depsxc,idir,nfft,rhonow)
649 [ + + ]: 1787340303 : do ifft=1,nfft
650 : 1787340303 : 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 [ + + ]: 254950836 : do ifft=1,nfft
658 : 254883420 : rho_tot=rhonow(ifft,1,1+idir)
659 : 254883420 : rho_up =rhonow(ifft,2,1+idir)
660 : 254883420 : rhonow(ifft,1,1+idir)=rho_up *depsxc(ifft,3) + rho_tot*depsxc(ifft,5)
661 : 254950836 : 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 : 203663 : 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 : 618045 : 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 : 618045 : numneg=0
734 : 618045 : worst=zero
735 : :
736 [ + + ]: 618045 : 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 [ + + ]: 1835742152 : do ifft=1,nfft
741 : 1835258176 : rhotmp=rhonow(ifft,1)
742 [ + + ]: 1835742152 : if(rhotmp<xc_denpos)then
743 [ + + ]: 60040945 : if(rhotmp<-xc_denpos)then
744 : : ! This case is probably beyond machine precision considerations
745 : 58380589 : worst=min(worst,rhotmp)
746 : 58380589 : numneg=numneg+1
747 : : end if
748 : 60040945 : rhonow(ifft,1)=xc_denpos
749 : : end if
750 : : end do
751 [ + - ]: 134069 : else if (nspden==2) then
752 : :
753 : : ! Spin-polarized
754 : :
755 : : ! rhonow is stored as (up,dn)
756 [ + + ]: 134069 : 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 [ + + ]: 122981378 : do ifft=1,nfft
761 : : ! For polarized case, rho(1) is spin-up density, rho(2) is spin-down density
762 : 122856745 : rho(1)=rhonow(ifft,1)
763 : 122856745 : rho(2)=rhonow(ifft,2)
764 [ + + ]: 368694868 : do ispden=1,nspden
765 [ + + ]: 368570235 : 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 [ + - ]: 9436 : 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 [ + + ]: 212178002 : do ifft=1,nfft
783 : : ! For polarized case, rho(1) is spin-up density, rho(2) is spin-down density
784 : 212168566 : rho(1)=rhonow(ifft,2)
785 : 212168566 : rho(2)=rhonow(ifft,1)-rho(1)
786 [ + + ]: 636515134 : do ispden=1,nspden
787 [ + + ]: 636505698 : if (rho(ispden)<xc_denpos) then
788 [ + + ]: 17365681 : if (rho(ispden)<-xc_denpos) then
789 : : ! This case is probably beyond machine precision considerations
790 : 1546853 : worst=min(worst,rho(ispden))
791 : 1546853 : numneg=numneg+1
792 : : end if
793 : 17365681 : rho(ispden)=xc_denpos
794 : 17365681 : rhonow(ifft,1)=rho(1)+rho(2)
795 : 17365681 : 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 [ + + ]: 618045 : if (numneg>0) then
807 [ + + ]: 1841 : if (iwarn==0) then
808 : : write(message,'(a,i0,a,a,a,es10.2,a,e10.2,11a)')&
809 : 1436 : & 'Density went too small (lower than xc_denpos) at ',numneg,' points',ch10,&
810 : 1436 : & 'and was set to xc_denpos = ',xc_denpos,'. Lowest was ',worst,'.',ch10,&
811 : 1436 : & 'This might be due to (1) too low boxcut or (2) too low ecut for',ch10,&
812 : 1436 : & ' pseudopotential core charge, or (3) too low ecut for estimated initial density.',ch10,&
813 : 1436 : & ' Possible workarounds : increase ecut, or define the input variable densty,',ch10,&
814 : 1436 : & ' with a value larger than the guess for the decay length, or initialize your,',ch10,&
815 : 2872 : & ' density with a preliminary LDA or GGA-PBE if you are using a more exotic xc functional.'
816 : 1436 : ABI_WARNING(message)
817 : : end if
818 : 1841 : iwarn=iwarn+1
819 : : end if
820 : :
821 : 618045 : 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 : 901188 : subroutine drivexc(ixc,order,npts,nspden,usegradient,uselaplacian,usekden,&
932 : 901188 : & rho_updn,exc,tsxc,vxcrho,nvxcgrho,nvxclrho,nvxctau,ndvxc,nd2vxc,el_temp, & ! mandatory arguments
933 : 4851689 : & 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 : 1802376 : real(dp),allocatable :: exci_rpa(:),rhotot(:),rspts(:),vxci_rpa(:,:),zeta(:)
964 : 901188 : real(dp),allocatable :: exc_c(:),exc_x(:),vxcrho_c(:,:),vxcrho_x(:,:)
965 : 901188 : real(dp),allocatable :: d2vxc_c(:,:),d2vxc_x(:,:),dvxc_c(:,:),dvxc_x(:,:)
966 : 901188 : real(dp),allocatable :: vxcgrho_x(:,:)
967 [ + + + + ]: 4505940 : type(libxc_functional_type) :: xc_funcs_vwn3(2),xc_funcs_lyp(2)
968 : : ! *************************************************************************
969 : :
970 : : !optional arguments
971 [ + + ]: 901188 : my_exexch=0;if(present(exexch)) my_exexch=exexch
972 : 901188 : my_hyb_mixing=0
973 [ + + ]: 901188 : if (ixc==41) my_hyb_mixing=quarter
974 [ + + ]: 901188 : if (ixc==42) my_hyb_mixing=third
975 [ + - ]: 901188 : if (present(hyb_mixing)) my_hyb_mixing=hyb_mixing
976 : :
977 : : ! =================================================
978 : : ! == Compatibility tests ==
979 : : ! =================================================
980 : :
981 : : !Check libXC initialization
982 [ + + ]: 901188 : if (ixc<0 .or. ixc==1402) then
983 : 383834 : 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 [ + + ]: 901188 : if (ixc<0) then
989 : : ixc_from_lib=libxc_functionals_ixc()
990 [ + + ]: 383834 : if (present(xc_funcs)) then
991 : 1024 : ixc_from_lib=libxc_functionals_ixc(xc_functionals=xc_funcs)
992 : : else
993 : 382810 : ixc_from_lib=libxc_functionals_ixc()
994 : : end if
995 [ - + ]: 383834 : 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 [ + - - + ]: 901188 : 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 : 901188 : has_gradient=.false.;has_laplacian=.false.;has_kden=.false.
1012 [ + + ]: 901188 : if (usegradient==1) then
1013 [ - + ]: 374305 : if (.not.present(grho2_updn)) then
1014 : 0 : ABI_BUG('missing grho2_updn argument!')
1015 : : end if
1016 [ + - ]: 374305 : if (nvxcgrho>0) then
1017 [ - + ]: 374305 : if (.not.present(vxcgrho)) then
1018 : 0 : ABI_BUG('missing vxcgrho argument!')
1019 : : end if
1020 : : has_gradient=.true.
1021 : : end if
1022 [ + - ]: 526883 : else if (nvxcgrho>0) then
1023 : 0 : ABI_BUG('nvxcgrho>0 and usegradient=0!')
1024 : : end if
1025 [ + + ]: 901188 : 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 [ + - ]: 865909 : else if (nvxclrho>0) then
1036 : 0 : ABI_BUG('nvxclrho>0 and uselaplacian=0!')
1037 : : end if
1038 [ + + ]: 901188 : 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 [ + - ]: 749500 : else if (nvxctau>0) then
1049 : 0 : ABI_BUG('nvxctau>0 and usekden=0!')
1050 : : end if
1051 [ + + ]: 901188 : if (abs(order)>=2) then
1052 [ - + ]: 149095 : if (.not.present(dvxc)) then
1053 : 0 : message='order>=2 needs argument dvxc!'
1054 : 0 : ABI_BUG(message)
1055 [ - + ]: 149095 : 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 [ + + ]: 901188 : 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 [ + + ]: 901188 : 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 : 894313 : & nvxctau=need_nvxctau,ndvxc=need_ndvxc,nd2vxc=need_nd2vxc)
1082 : : end if
1083 : : if ((has_gradient.and.need_gradient>usegradient).or.&
1084 [ + + + - : 901188 : & (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 [ + + + - : 901188 : & 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 [ + + - + ]: 901188 : 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 [ + + ]: 2044231740 : 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 [ + + ]: 901188 : & ixc==21.or.ixc==22.or.ixc==50.or.ixc==51.or.ixc==60) then
1123 [ - + + - : 599259 : ABI_MALLOC(rhotot,(npts))
+ - + - -
+ - + ]
1124 [ - + + - : 399506 : ABI_MALLOC(rspts,(npts))
- + - + ]
1125 [ + + ]: 199753 : if(nspden==1)then
1126 [ + + ]: 518638316 : rhotot(:)=two*rho_updn(:,1)
1127 : : else
1128 [ + + ]: 79535808 : rhotot(:)=rho_updn(:,1)+rho_updn(:,2)
1129 : : end if
1130 : 199753 : call invcb(rhotot,rspts,npts)
1131 [ + + ]: 598174124 : rspts(:)=rsfac*rspts(:)
1132 : : end if
1133 : :
1134 : : !If needed, compute zeta
1135 [ + + + + ]: 901188 : if (ixc==1.or.ixc==21.or.ixc==22) then
1136 [ - + + - : 314844 : ABI_MALLOC(zeta,(npts))
+ - + - -
+ - + ]
1137 [ + + ]: 104948 : if(nspden==1)then
1138 [ + + ]: 270032157 : zeta(:)=zero
1139 : : else
1140 [ + + ]: 80332048 : 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 [ - + + + : 901188 : 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 [ + + ]: 104948 : if (order**2 <= 1) then
1161 : 99284 : call xcspol(exc,npts,nspden,order,rspts,vxcrho,zeta,ndvxc)
1162 : : else
1163 : 5664 : 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 [ + + + - : 701800 : 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 [ + + ]: 310955 : if(ixc==7)optpbe=1
1212 : : ! x-only part of Perdew-Wang
1213 [ + + ]: 310955 : if(ixc==8)optpbe=-1
1214 : : ! Exchange + RPA correlation from Perdew-Wang
1215 [ - + ]: 310955 : if(ixc==9)optpbe=3
1216 : : ! Perdew-Burke-Ernzerhof GGA
1217 [ + + ]: 310955 : if(ixc==11)optpbe=2
1218 : : ! x-only part of PBE
1219 [ + + ]: 310955 : if(ixc==12)optpbe=-2
1220 : : ! C09x exchange of V. R. Cooper
1221 [ + + ]: 310955 : if(ixc==24)optpbe=-4
1222 : : ! revPBE of Zhang and Yang
1223 [ + + ]: 310955 : if(ixc==14)optpbe=5
1224 : : ! RPBE of Hammer, Hansen and Norskov
1225 [ + + ]: 310955 : if(ixc==15)optpbe=6
1226 : : ! Wu and Cohen
1227 [ + + ]: 310955 : if(ixc==23)optpbe=7
1228 [ + + ]: 310955 : if (ixc >=7.and.ixc<=9) then
1229 [ + + ]: 233563 : if (order**2 <= 1) then
1230 : 196683 : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc)
1231 [ + + ]: 36880 : else if (order /=3) then
1232 : 34364 : 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 [ + - ]: 77392 : else if ((ixc >= 11 .and. ixc <= 15) .or. (ixc>=23 .and. ixc<=24)) then
1237 [ + + ]: 77392 : if (order**2 <= 1) then
1238 : : call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,&
1239 : 66402 : & dvxcdgr=vxcgrho,exexch=my_exexch,grho2_updn=grho2_updn)
1240 [ + - ]: 10990 : else if (order /=3) then
1241 [ + + ]: 10990 : 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 : 10987 : 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 [ + - ]: 383834 : else if( ixc<0 ) then
1594 : :
1595 : : ! ===== meta-GGA =====
1596 [ + + + + ]: 383834 : 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 [ + + ]: 232195 : else if (need_gradient==1) then
1686 [ + + ]: 138382 : if (abs(order)<=1) then
1687 [ + + ]: 111460 : 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 : 110436 : & grho2=grho2_updn,vxcgr=vxcgrho)
1694 : : end if
1695 : : else
1696 [ - + ]: 26922 : 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 : 26922 : & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc)
1703 : : end if
1704 : : end if
1705 : :
1706 : : ! ===== LDA =====
1707 : : else
1708 [ + + ]: 93813 : if (abs(order)<=1) then
1709 [ - + ]: 79477 : 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 : 79477 : call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho)
1714 : : end if
1715 [ + + ]: 14336 : else if (abs(order)<=2) then
1716 [ - + ]: 14334 : 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 : 14334 : & 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 [ + + ]: 901188 : if(allocated(rhotot)) then
1741 : 199753 : ABI_FREE(rhotot)
1742 : : end if
1743 [ + + ]: 901188 : if(allocated(rspts)) then
1744 : 199753 : ABI_FREE(rspts)
1745 : : end if
1746 [ + + ]: 901188 : if(allocated(zeta)) then
1747 : 104948 : ABI_FREE(zeta)
1748 : : end if
1749 : :
1750 [ + + + + : 4857540 : end subroutine drivexc
+ + + + +
+ + + + -
+ + ]
1751 : : !!***
1752 : :
1753 : : end module m_drivexc
1754 : : !!***
|
