TABLE OF CONTENTS
ABINIT/drivexc [ Functions ]
NAME
drivexc
FUNCTION
Driver of XC functionals. Treat spin-polarized as well as non-spin-polarized. Treat local approximations or GGAs. Optionally, deliver the XC kernel, or even the derivative of the XC kernel (the third derivative of the XC energy)
COPYRIGHT
Copyright (C) 2002-2018 ABINIT group (XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
ixc=number of the XC functional (to be described) ndvxc= size of dvxc(npts,ndvxc) ngr2= size of grho2_updn(npts,ngr2) nvxcgrho= size of vxcgrho(npts,nvxcgrho) npts=number of real space points on which the density (and its gradients, if needed) is provided nspden=number of spin-density components (1 or 2) order=gives the maximal derivative of Exc computed. 1=usual value (return exc and vxc) 2=also computes the kernel (return exc,vxc,kxc) -2=like 2, except (to be described) 3=also computes the derivative of the kernel (return exc,vxc,kxc,k3xc) rho_updn(npts,nspden)=the spin-up and spin-down densities If nspden=1, only the spin-up density must be given. In the calling routine, the spin-down density must be equal to the spin-up density, and both are half the total density. If nspden=2, the spin-up and spin-down density must be given Optional inputs: [el_temp]= electronic temperature (to be used for finite temperature XC functionals) [exexch]= choice of <<<local>>> exact exchange. Active if exexch=3 (only for GGA, and NOT for libxc) [grho2_updn(npts,ngr2)]=the square of the gradients of spin-up, spin-down, and total density If nspden=1, only the square of the gradient of the spin-up density must be given. In the calling routine, the square of the gradient of the spin-down density must be equal to the square of the gradient of the spin-up density, and both must be equal to one-quarter of the square of the gradient of the total density. If nspden=2, the square of the gradients of spin-up, spin-down, and total density must be given. Note that the square of the gradient of the total density is usually NOT related to the square of the gradient of the spin-up and spin-down densities, because the gradients are not usually aligned. This is not the case when nspden=1. [hyb_mixing]= mixing parameter for the native PBEx functionals (ixc=41 and 42) [lrho_updn(npts,nspden)]=the Laplacian of spin-up and spin-down densities If nspden=1, only the spin-up Laplacian density must be given. In the calling routine, the spin-down Laplacian density must be equal to the spin-up Laplacian density, and both are half the total Laplacian density. If nspden=2, the Laplacian of spin-up and spin-down densities must be given [tau_updn(npts,nspden)]=the spin-up and spin-down kinetic energy densities If nspden=1, only the spin-up kinetic energy density must be given. In the calling routine, the spin-down kinetic energy density must be equal to the spin-up kinetic energy density, and both are half the total kinetic energy density. If nspden=2, the spin-up and spin-down kinetic energy densities must be given [xc_funcs(2)]= <type(libxc_functional_type)>, optional : libxc XC functionals. If not specified, the underlying xc_global(2) is used by libxc. [xc_tb09_c]= c parameter for the Tran-Blaha mGGA functional, within libxc
OUTPUT
exc(npts)=exchange-correlation energy density (hartree) vxcrho(npts,nspden)= (d($\rho$*exc)/d($\rho_up$)) (hartree) and (d($\rho$*exc)/d($\rho_down$)) (hartree) vxcgrho(npts,3)= 1/$|grad \rho_up|$ (d($\rho$*exc)/d($|grad \rho_up|$)) (hartree) 1/$|grad \rho_dn|$ (d($\rho$*exc)/d($|grad \rho_dn|$)) (hartree) and 1/$|grad \rho|$ (d($\rho$*exc)/d($|grad \rho|$)) (hartree) (will be zero if a LDA functional is used) vxclrho(npts,nspden)=(only for meta-GGA, i.e. optional output)= (d($\rho$*exc)/d($\lrho_up$)) (hartree) and (d($\rho$*exc)/d($\lrho_down$)) (hartree) vxctau(npts,nspden)=(only for meta-GGA, i.e. optional output)= derivative of XC energy density with respect to kinetic energy density (depsxcdtau). (d($\rho$*exc)/d($\tau_up$)) (hartree) and (d($\rho$*exc)/d($\tau_down$)) (hartree) Optional output: if(abs(order)>1) [dvxc(npts,ndvxc)]=partial second derivatives of the xc energy (This is a mess, to be rationalized !!) In case of local energy functional (option=1,-1 or 3): dvxc(npts,1+nspden)= (Hartree*bohr^3) if(nspden=1 .and. order==2): dvxci(:,1)=dvxc/d$\rho$ , dvxc(:,2) empty if(nspden=1 .and. order==-2): also compute dvxci(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$ if(nspden=2): dvxc(:,1)=dvxc($\uparrow$)/d$\rho(\downarrow)$, dvxc(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$, dvxc(:,3)=dvxc($\downarrow$)/d$\rho(\downarrow)$ In case of gradient corrected functional (option=2,-2, 4, -4, 5, 6, 7): dvxc(npts,15)= dvxc(:,1)= d2Ex/drho_up drho_up dvxc(:,2)= d2Ex/drho_dn drho_dn dvxc(:,3)= dEx/d(abs(grad(rho_up))) / abs(grad(rho_up)) dvxc(:,4)= dEx/d(abs(grad(rho_dn))) / abs(grad(rho_dn)) dvxc(:,5)= d2Ex/d(abs(grad(rho_up))) drho_up / abs(grad(rho_up)) dvxc(:,6)= d2Ex/d(abs(grad(rho_dn))) drho_dn / abs(grad(rho_dn)) dvxc(:,7)= 1/abs(grad(rho_up)) * d/d(abs(grad(rho_up)) (dEx/d(abs(grad(rho_up))) /abs(grad(rho_up))) dvxc(:,8)= 1/abs(grad(rho_dn)) * d/d(abs(grad(rho_dn)) (dEx/d(abs(grad(rho_dn))) /abs(grad(rho_dn))) dvxc(:,9)= d2Ec/drho_up drho_up dvxc(:,10)=d2Ec/drho_up drho_dn dvxc(:,11)=d2Ec/drho_dn drho_dn dvxc(:,12)=dEc/d(abs(grad(rho))) / abs(grad(rho)) dvxc(:,13)=d2Ec/d(abs(grad(rho))) drho_up / abs(grad(rho)) dvxc(:,14)=d2Ec/d(abs(grad(rho))) drho_dn / abs(grad(rho)) dvxc(:,15)=1/abs(grad(rho)) * d/d(abs(grad(rho)) (dEc/d(abs(grad(rho))) /abs(grad(rho))) if(abs(order)>2) (only available for LDA and nspden=1) [d2vxc(npts,nd2vxc)]=partial third derivatives of the xc energy if nspden=1 d2vxc(npts,1)=second derivative of the XC potential=3rd order derivative of energy if nspden=2 d2vxc(npts,1), d2vxc(npts,2), d2vxc(npts,3), d2vxc(npts,4) (3rd derivative of energy) [fxcT(npts)]=XC free energy of the electron gaz at finite temperature (to be used for plasma systems)
PARENTS
drivexc_main
CHILDREN
invcb,libxc_functionals_end,libxc_functionals_getvxc libxc_functionals_init,xchcth,xchelu,xciit,xclb,xcpbe,xcpzca,xcspol xctetr,xcwign,xcxalp
SOURCE
134 #if defined HAVE_CONFIG_H 135 #include "config.h" 136 #endif 137 138 #include "abi_common.h" 139 140 subroutine drivexc(exc,ixc,npts,nspden,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 141 & dvxc,d2vxc,grho2_updn,vxcgrho,el_temp,exexch,fxcT,& !Optional arguments 142 & hyb_mixing,lrho_updn,vxclrho,tau_updn,vxctau,xc_funcs,xc_tb09_c) !Optional arguments 143 144 145 use defs_basis 146 use m_profiling_abi 147 use m_errors 148 use libxc_functionals 149 150 !This section has been created automatically by the script Abilint (TD). 151 !Do not modify the following lines by hand. 152 #undef ABI_FUNC 153 #define ABI_FUNC 'drivexc' 154 use interfaces_41_xc_lowlevel, except_this_one => drivexc 155 !End of the abilint section 156 157 implicit none 158 159 !Arguments ------------------------------------ 160 !scalars 161 integer,intent(in) :: ixc,ndvxc,ngr2,nd2vxc,npts,nspden,nvxcgrho,order 162 integer,intent(in),optional :: exexch 163 real(dp),intent(in),optional :: el_temp,hyb_mixing,xc_tb09_c 164 !arrays 165 real(dp),intent(in) :: rho_updn(npts,nspden) 166 real(dp),intent(in),optional :: grho2_updn(npts,ngr2) 167 real(dp),intent(in),optional :: lrho_updn(npts,nspden), tau_updn(npts,nspden) 168 real(dp),intent(out) :: exc(npts),vxcrho(npts,nspden) 169 real(dp),intent(out),optional :: d2vxc(npts,nd2vxc),dvxc(npts,ndvxc),fxcT(:) 170 real(dp),intent(out),optional :: vxcgrho(npts,nvxcgrho) 171 real(dp),intent(out),optional :: vxclrho(npts,nspden),vxctau(npts,nspden) 172 type(libxc_functional_type),intent(inout),optional :: xc_funcs(2) 173 174 !Local variables------------------------------- 175 !scalars 176 integer :: exexch_,ixc_from_lib,ixc1,ixc2,ndvxc_x,optpbe,ispden 177 logical :: libxc_test,xc_err_ndvxc,xc_err_nvxcgrho1,xc_err_nvxcgrho2 178 logical :: is_gga,is_mgga 179 real(dp) :: alpha 180 real(dp),parameter :: rsfac=0.6203504908994000e0_dp 181 character(len=500) :: message 182 !arrays 183 real(dp),allocatable :: exci_rpa(:) 184 real(dp),allocatable :: rhotot(:),rspts(:),vxci_rpa(:,:),zeta(:) 185 real(dp),allocatable :: exc_c(:),exc_x(:),vxcrho_c(:,:),vxcrho_x(:,:) 186 real(dp),allocatable :: d2vxc_c(:,:),d2vxc_x(:,:),dvxc_c(:,:),dvxc_x(:,:) 187 real(dp),allocatable :: vxcgrho_x(:,:) 188 type(libxc_functional_type) :: xc_funcs_vwn3(2),xc_funcs_lyp(2) 189 190 ! ************************************************************************* 191 192 ! ================================================= 193 ! == Compatibility tests == 194 ! ================================================= 195 196 !Checks the values of order and other size parameters 197 if( (order<1 .and. order/=-2) .or. order>4)then 198 write(message, '(a,i0)' )& 199 & 'The only allowed values for order are 1,2,-2 or 3, while it is found to be ',order 200 MSG_BUG(message) 201 end if 202 xc_err_ndvxc=.false. 203 xc_err_nvxcgrho1=.false. 204 xc_err_nvxcgrho2=.false. 205 if (ixc==16.or.ixc==17.or.ixc==26.or.ixc==27) xc_err_nvxcgrho1=(nvxcgrho/=2) 206 if (order**2>1) then 207 if ((ixc>=2.and.ixc<=6).or.(ixc==50)) xc_err_ndvxc=(ndvxc/=1) 208 if (ixc==1.or.ixc==21.or.ixc==22) xc_err_ndvxc=(ndvxc/=nspden+1) 209 end if 210 if (order==2) then 211 if ((ixc>=7.and.ixc<=10).or.(ixc==13)) xc_err_nvxcgrho1=(ndvxc/=1+nspden.or.nvxcgrho/=0) 212 if (ixc==12.or.ixc==24) xc_err_nvxcgrho1=(ndvxc/=8.or.nvxcgrho/=3) 213 if (ixc==11.or.ixc==14.or.ixc==15.or.ixc==23.or.ixc==41.or.ixc==42) xc_err_nvxcgrho1=(ndvxc/=15.or.nvxcgrho/=3) 214 end if 215 if (order==3) then 216 if (ixc==3) xc_err_ndvxc=(ndvxc/=1) 217 if (ixc==10.or.ixc==13) xc_err_nvxcgrho1=(ndvxc/=1+nspden.or.nvxcgrho/=0) 218 if (ixc==12.or.ixc==24) xc_err_nvxcgrho1=(ndvxc/=8.or.nvxcgrho/=3) 219 if (ixc==11.or.ixc==14.or.ixc==15.or.ixc==23.or.ixc==41.or.ixc==42) xc_err_nvxcgrho1=(ndvxc/=15.or.nvxcgrho/=3) 220 if (ixc>=7.and.ixc<=9) xc_err_nvxcgrho2=(ndvxc/=1+nspden.or.nvxcgrho/=0.or.nd2vxc/=(3*nspden-2)) 221 if (ixc==50) xc_err_ndvxc=(nd2vxc/=0) 222 end if 223 if (xc_err_ndvxc) then 224 write(message, '(7a,i0,a,i0,a)' )& 225 & 'Wrong value of ndvxc:',ch10,& 226 & 'the value of the spin polarization',ch10,& 227 & 'is not compatible with the value of ixc:',ch10,& 228 & 'ixc=',ixc,' (nspden=',nspden,')' 229 MSG_BUG(message) 230 end if 231 if (xc_err_nvxcgrho1) then 232 write(message,'(3a,i0,a,i0,a,i0)' )& 233 & 'Wrong value of ndvxc or nvxcgrho:',ch10,& 234 & 'ixc=',ixc,'ndvxc=',ndvxc,'nvxcgrho=',nvxcgrho 235 MSG_BUG(message) 236 end if 237 if (xc_err_nvxcgrho2) then 238 write(message, '(3a,i0,a,i0,a,i0,a,i0)' )& 239 & 'Wrong value of ndvxc or nvxcgrho or nd2vxc:',ch10,& 240 & 'ixc=',ixc,'ndvxc=',ndvxc,'nvxcgrho=',nvxcgrho,'nd2vxc=',nd2vxc 241 MSG_BUG(message) 242 end if 243 244 !Check libXC 245 if (ixc<0 .or. ixc==1402) then 246 libxc_test=libxc_functionals_check(stop_if_error=.true.) 247 end if 248 249 if (ixc<0) then 250 ! Prepare the tests 251 if(present(xc_funcs))then 252 is_gga=libxc_functionals_isgga(xc_functionals=xc_funcs) 253 is_mgga=libxc_functionals_ismgga(xc_functionals=xc_funcs) 254 ixc_from_lib=libxc_functionals_ixc(xc_functionals=xc_funcs) 255 else 256 is_gga=libxc_functionals_isgga() 257 is_mgga=libxc_functionals_ismgga() 258 ixc_from_lib=libxc_functionals_ixc() 259 end if 260 ! Check consistency between ixc passed in input and the one used to initialize the library. 261 if (ixc /= ixc_from_lib) then 262 write(message, '(a,i0,2a,i0,2a)')& 263 & 'The value of ixc specified in input, ixc = ',ixc,ch10,& 264 & 'differs from the one used to initialize the functional ',ixc_from_lib,ch10,& 265 & 'Action: reinitialize the global structure funcs, see NOTES in m_libxc_functionals' 266 MSG_BUG(message) 267 end if 268 else if (ixc==1420)then 269 is_gga=.true. 270 is_mgga=.false. 271 end if 272 273 if (ixc<0 .or. ixc==1402) then 274 ! Check whether all the necessary arrays are present and have the correct dimensions 275 if (is_gga .or. is_mgga) then 276 if ( (.not. present(grho2_updn)) .or. (.not. present(vxcgrho))) then 277 write(message, '(5a,2L2,2a,2L2,2a,i10,a,i5,a,i5)' )& 278 & 'At least one of the functionals is a GGA or a MGGA,',ch10,& 279 & 'but not all the necessary arrays are present.',ch10,& 280 & 'is_gga, is_mgga=',is_gga,is_mgga,ch10,& 281 & 'present(grho2_updn),present(vxcgrho)=',present(grho2_updn),present(vxcgrho),ch10,& 282 & 'ixc=',ixc,' nvxcgrho=',nvxcgrho,' ngr2=',ngr2 283 MSG_BUG(message) 284 end if 285 if (ngr2==0.or.nvxcgrho/=3) then 286 write(message, '(5a,i7,a,i6,a,i6)' )& 287 & 'The value of the number of the XC functional ixc',ch10,& 288 & 'is not compatible with the value of nvxcgrho or ngr2',ch10,& 289 & 'ixc=',ixc,' nvxcgrho=',nvxcgrho,' ngr2=',ngr2 290 MSG_BUG(message) 291 end if 292 end if 293 if (is_mgga) then 294 if ( (.not. present(lrho_updn)) .or. (.not. present(vxclrho)) .or. & 295 (.not. present(tau_updn)) .or. (.not. present(vxctau)) ) then 296 write(message, '(5a,i7)' )& 297 & 'At least one of the functionals is a MGGA,',ch10,& 298 & 'but not all the necessary arrays are present.',ch10,& 299 & 'ixc=',ixc 300 MSG_BUG(message) 301 end if 302 end if 303 end if 304 305 !Checks the compatibility between the inputs and the presence of the optional arguments 306 if(present(dvxc))then 307 if(order**2<=1.or.ixc==16.or.ixc==17.or.ixc==26.or.ixc==27)then 308 write(message, '(5a,i0,a,i0)' )& 309 & 'The value of the number of the XC functional ixc',ch10,& 310 & 'or the value of order is not compatible with the presence of the array dvxc',ch10,& 311 & 'ixc=',ixc,'order=',order 312 MSG_BUG(message) 313 end if 314 end if 315 if(present(d2vxc))then 316 if(order/=3.or.(ixc/=3.and.(((ixc>15).and.(ixc/=23)) .or.& 317 & (ixc>=0.and.ixc<7).or.(ixc>=40 .and. ixc/=1402))) )then 318 write(message, '(5a,i6,a,i0)' )& 319 & 'The value of the number of the XC functional ixc',ch10,& 320 & 'or the value of order is not compatible with the presence of the array d2vxc',ch10,& 321 & 'ixc=',ixc,'order=',order 322 MSG_BUG(message) 323 end if 324 end if 325 if(present(vxcgrho))then 326 if(nvxcgrho==0.or. & 327 & ((((ixc>17.and.ixc/=23.and.ixc/=24.and.ixc/=26.and.ixc/=27) .or.& 328 & (ixc>= 0.and.ixc<7).or.(ixc>=40 .and. ixc/=1402)) .and. nvxcgrho /=3 ))) then 329 if(ixc<31.and.ixc>34)then !! additional if to include ixc 31 to 34 in the list (these ixc are used for mgga test see below) 330 write(message, '(5a,i0,a,i0)' )& 331 & 'The value of the number of the XC functional ixc',ch10,& 332 & 'or the value of nvxcgrho is not compatible with the presence of the array vxcgrho',ch10,& 333 & 'ixc=',ixc,' nvxcgrho=',nvxcgrho 334 MSG_BUG(message) 335 end if 336 end if 337 end if 338 if(present(grho2_updn))then 339 if (ngr2/=2*nspden-1 ) then 340 write(message, '(4a)' ) ch10,& 341 & 'drivexc : BUG -',ch10,& 342 & 'ngr2 must be 2*nspden-1 !' 343 ! MSG_BUG(message) 344 end if 345 if ((ixc>0.and.ixc<11).or.(ixc>17.and.ixc<23).or.ixc==25.or.(ixc>27.and.ixc<31).or. & 346 & (ixc>34.and.ixc<41).or.ixc==50) then 347 write(message, '(5a,i0)' )& 348 & 'The value of the number of the XC functional ixc',ch10,& 349 & 'is not compatible with the presence of the array grho2_updn',ch10,& 350 & 'ixc=',ixc 351 MSG_BUG(message) 352 end if 353 end if 354 if (present(exexch)) then 355 if(exexch/=0.and.(.not.present(grho2_updn))) then 356 message='exexch argument only valid for GGA!' 357 MSG_BUG(message) 358 end if 359 end if 360 if(ixc==31.or.ixc==32.or.ixc==33.or.ixc==34) then 361 if((.not.(present(vxcgrho))) .or. (.not.(present(vxclrho))) .or. (.not.(present(vxctau))))then 362 message = 'vxcgrho or vxclrho or vxctau is not present but they are all needed for MGGA XC tests.' 363 MSG_BUG(message) 364 end if 365 end if 366 if(ixc==50) then 367 if(.not.(present(el_temp)).or.(.not.present(fxcT)))then 368 message = 'el_temp or fxcT is not present but are needed for IIT XC functional.' 369 MSG_BUG(message) 370 end if 371 if (size(fxcT)/=npts) then 372 MSG_BUG('fxcT size must be npts!') 373 end if 374 end if 375 376 ! ================================================= 377 ! == Intermediate quantities computation == 378 ! ================================================= 379 380 !If needed, compute rhotot and rs 381 if (ixc==1.or.ixc==2.or.ixc==3.or.ixc==4.or.ixc==5.or.ixc==6.or.ixc==21.or.ixc==22.or.ixc==50) then 382 ABI_ALLOCATE(rhotot,(npts)) 383 ABI_ALLOCATE(rspts,(npts)) 384 if(nspden==1)then 385 rhotot(:)=two*rho_updn(:,1) 386 else 387 rhotot(:)=rho_updn(:,1)+rho_updn(:,2) 388 end if 389 call invcb(rhotot,rspts,npts) 390 rspts(:)=rsfac*rspts(:) 391 end if 392 393 !If needed, compute zeta 394 if (ixc==1.or.ixc==21.or.ixc==22) then 395 ABI_ALLOCATE(zeta,(npts)) 396 if(nspden==1)then 397 zeta(:)=zero 398 else 399 zeta(:)=two*rho_updn(:,1)/rhotot(:)-one 400 end if 401 end if 402 403 ! ================================================= 404 ! == XC energy, potentiel, ... computation == 405 ! ================================================= 406 407 exexch_=0;if(present(exexch)) exexch_=exexch 408 409 !>>>>> No exchange-correlation 410 if (ixc==0.or.ixc==40) then 411 exc=zero ; vxcrho=zero 412 if(present(d2vxc)) d2vxc(:,:)=zero 413 if(present(dvxc)) dvxc(:,:)=zero 414 if(present(vxcgrho)) vxcgrho(:,:)=zero 415 if(present(vxclrho)) vxclrho(:,:)=zero 416 if(present(vxctau)) vxctau(:,:)=zero 417 418 !>>>>> New Teter fit (4/93) to Ceperley-Alder data, with spin-pol option 419 else if (ixc==1 .or. ixc==21 .or. ixc==22) then 420 ! new Teter fit (4/93) to Ceperley-Alder data, with spin-pol option 421 if (order**2 <= 1) then 422 call xcspol(exc,npts,nspden,order,rspts,vxcrho,zeta,ndvxc) 423 else 424 call xcspol(exc,npts,nspden,order,rspts,vxcrho,zeta,ndvxc,dvxc) 425 end if 426 427 !>>>>> Perdew-Zunger fit to Ceperly-Alder data (no spin-pol) 428 else if (ixc==2) then 429 if (order**2 <= 1) then 430 call xcpzca(exc,npts,order,rhotot,rspts,vxcrho(:,1)) 431 else 432 call xcpzca(exc,npts,order,rhotot,rspts,vxcrho(:,1),dvxc) 433 end if 434 435 !>>>>> Teter fit (4/91) to Ceperley-Alder values (no spin-pol) 436 else if (ixc==3) then 437 if (order**2 <= 1) then 438 call xctetr(exc,npts,order,rhotot,rspts,vxcrho(:,1)) 439 else if (order == 2) then 440 call xctetr(exc,npts,order,rhotot,rspts,vxcrho(:,1),dvxc=dvxc) 441 else if (order == 3) then 442 call xctetr(exc,npts,order,rhotot,rspts,vxcrho(:,1),d2vxc=d2vxc,dvxc=dvxc) 443 end if 444 445 !>>>>> Wigner xc (no spin-pol) 446 else if (ixc==4) then 447 if (order**2 <= 1) then 448 call xcwign(exc,npts,order,rspts,vxcrho(:,1)) 449 else 450 call xcwign(exc,npts,order,rspts,vxcrho(:,1),dvxc) 451 end if 452 453 !>>>>> Hedin-Lundqvist xc (no spin-pol) 454 else if (ixc==5) then 455 if (order**2 <= 1) then 456 call xchelu(exc,npts,order,rspts,vxcrho(:,1)) 457 else 458 call xchelu(exc,npts,order,rspts,vxcrho(:,1),dvxc) 459 end if 460 461 !>>>>> X-alpha (no spin-pol) 462 else if (ixc==6) then 463 if (order**2 <= 1) then 464 call xcxalp(exc,npts,order,rspts,vxcrho(:,1)) 465 else 466 call xcxalp(exc,npts,order,rspts,vxcrho(:,1),dvxc) 467 end if 468 469 !>>>>> PBE and alternatives 470 else if (((ixc>=7.and.ixc<=15).or.(ixc>=23.and.ixc<=24)).and.ixc/=10.and.ixc/=13) then 471 ! Perdew-Wang LSD is coded in Perdew-Burke-Ernzerhof GGA, with optpbe=1 472 if(ixc==7)optpbe=1 473 ! x-only part of Perdew-Wang 474 if(ixc==8)optpbe=-1 475 ! Exchange + RPA correlation from Perdew-Wang 476 if(ixc==9)optpbe=3 477 ! Perdew-Burke-Ernzerhof GGA 478 if(ixc==11)optpbe=2 479 ! x-only part of PBE 480 if(ixc==12)optpbe=-2 481 ! C09x exchange of V. R. Cooper 482 if(ixc==24)optpbe=-4 483 ! revPBE of Zhang and Yang 484 if(ixc==14)optpbe=5 485 ! RPBE of Hammer, Hansen and Norskov 486 if(ixc==15)optpbe=6 487 ! Wu and Cohen 488 if(ixc==23)optpbe=7 489 if (ixc >=7.and.ixc<=9) then 490 if (order**2 <= 1) then 491 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 492 else if (order /=3) then 493 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,dvxci=dvxc) 494 else if (order ==3) then 495 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,d2vxci=d2vxc,dvxci=dvxc) 496 end if 497 else if ((ixc >= 11 .and. ixc <= 15) .or. (ixc>=23 .and. ixc<=24)) then 498 if (order**2 <= 1) then 499 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 500 & dvxcdgr=vxcgrho,exexch=exexch_,grho2_updn=grho2_updn) 501 else if (order /=3) then 502 if(ixc==12 .or. ixc==24)then 503 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 504 & dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 505 else if(ixc/=12 .or. ixc/=24) then 506 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 507 & dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 508 end if 509 else if (order ==3) then 510 if(ixc==12 .or. ixc==24)then 511 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 512 & d2vxci=d2vxc,dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 513 else if(ixc/=12 .or. ixc/=24) then 514 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 515 & d2vxci=d2vxc,dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 516 end if 517 end if 518 end if 519 520 !>>>>> RPA correlation from Perdew-Wang 521 else if (ixc==10) then 522 if (order**2 <= 1) then 523 ABI_ALLOCATE(exci_rpa,(npts)) 524 ABI_ALLOCATE(vxci_rpa,(npts,2)) 525 optpbe=3 526 call xcpbe(exci_rpa,npts,nspden,optpbe,order,rho_updn,vxci_rpa,ndvxc,ngr2,nd2vxc) 527 optpbe=1 528 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 529 exc(:)=exc(:)-exci_rpa(:) 530 ! PMA: second index of vxcrho is nspden while that of rpa is 2 they can mismatch 531 vxcrho(:,1:min(nspden,2))=vxcrho(:,1:min(nspden,2))-vxci_rpa(:,1:min(nspden,2)) 532 ABI_DEALLOCATE(exci_rpa) 533 ABI_DEALLOCATE(vxci_rpa) 534 else if (order /=3) then 535 ABI_ALLOCATE(exci_rpa,(npts)) 536 ABI_ALLOCATE(vxci_rpa,(npts,2)) 537 optpbe=3 538 call xcpbe(exci_rpa,npts,nspden,optpbe,order,rho_updn,vxci_rpa,ndvxc,ngr2,nd2vxc,dvxci=dvxc) 539 optpbe=1 540 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,dvxci=dvxc) 541 exc(:)=exc(:)-exci_rpa(:) 542 vxcrho(:,:)=vxcrho(:,:)-vxci_rpa(:,:) 543 ABI_DEALLOCATE(exci_rpa) 544 ABI_DEALLOCATE(vxci_rpa) 545 else if (order ==3) then 546 ABI_ALLOCATE(exci_rpa,(npts)) 547 ABI_ALLOCATE(vxci_rpa,(npts,2)) 548 optpbe=3 549 call xcpbe(exci_rpa,npts,nspden,optpbe,order,rho_updn,vxci_rpa,ndvxc,ngr2,nd2vxc,& 550 & d2vxci=d2vxc,dvxci=dvxc) 551 optpbe=1 552 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 553 & d2vxci=d2vxc,dvxci=dvxc) 554 exc(:)=exc(:)-exci_rpa(:) 555 vxcrho(:,:)=vxcrho(:,:)-vxci_rpa(:,:) 556 ABI_DEALLOCATE(exci_rpa) 557 ABI_DEALLOCATE(vxci_rpa) 558 end if 559 560 !>>>>> LDA xc energy like ixc==7, and Leeuwen-Baerends GGA xc potential 561 else if(ixc==13) then 562 if (order**2 <= 1) then 563 optpbe=1 564 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 565 call xclb(grho2_updn,npts,nspden,rho_updn,vxcrho) 566 else if (order /=3) then 567 optpbe=1 568 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,dvxci=dvxc) 569 call xclb(grho2_updn,npts,nspden,rho_updn,vxcrho) 570 else if (order ==3) then 571 optpbe=1 572 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,d2vxci=d2vxc,dvxci=dvxc) 573 call xclb(grho2_updn,npts,nspden,rho_updn,vxcrho) 574 end if 575 576 !>>>>> HTCH93, HTCH120, HTCH107, HTCH147 577 else if(ixc==16 .or. ixc==17 .or. ixc==26 .or. ixc==27) then 578 call xchcth(vxcgrho,exc,grho2_updn,ixc,npts,nspden,order,rho_updn,vxcrho) 579 580 !>>>>> Only for test purpose (test various part of MGGA implementation) 581 else if(ixc==31 .or. ixc==32 .or. ixc==33 .or. ixc==34) then 582 exc(:)=zero 583 vxcrho(:,:)=zero 584 vxcgrho(:,:)=zero 585 vxclrho(:,:)=zero 586 vxctau(:,:)=zero 587 588 !>>>>> Perdew-Wang LSD is coded in Perdew-Burke-Ernzerhof GGA, with optpbe=1 589 optpbe=1 590 select case(ixc) 591 case (31) 592 alpha=1.00d0-(1.00d0/1.01d0) 593 ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc) 594 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 595 if (nspden==1) then 596 ! it should be : exc_tot= exc_spin up + exc_spin down = 2*exc_spin up but this applies to tau and rho (so it cancels) 597 exc(:)=exc(:)+alpha*tau_updn(:,1)/rho_updn(:,1) 598 else 599 do ispden=1,nspden 600 exc(:)=exc(:)+alpha*tau_updn(:,ispden)/(rho_updn(:,1)+rho_updn(:,2)) 601 end do 602 end if 603 vxctau(:,:)=alpha 604 case (32) 605 alpha=0.01d0 606 ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc) 607 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 608 if (nspden==1) then 609 exc(:)=exc(:)+2.0d0*alpha*lrho_updn(:,1) 610 vxcrho(:,1) =vxcrho(:,1)+2.0d0*alpha*lrho_updn(:,1) 611 vxclrho(:,1)=alpha*2.0d0*rho_updn(:,1) 612 else 613 do ispden=1,nspden 614 exc(:)=exc(:)+alpha*lrho_updn(:,ispden) 615 vxcrho(:,ispden) =vxcrho(:,ispden)+alpha*(lrho_updn(:,1)+lrho_updn(:,2)) 616 vxclrho(:,ispden)=alpha*(rho_updn(:,1)+rho_updn(:,2)) 617 end do 618 end if 619 case (33) 620 alpha=-0.010d0 621 ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc) 622 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 623 if (nspden==1) then 624 ! it should be : exc_tot= exc_spin up + exc_spin down = 2*exc_spin up but this applies to grho2 and rho 625 ! (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) 626 exc(:)=exc(:)+alpha*2.0d0*grho2_updn(:,1)/rho_updn(:,1) 627 if(nvxcgrho==2)vxcgrho(:,1:2)=2.0d0*alpha 628 if(nvxcgrho==3)vxcgrho(:,3)=2.0d0*alpha 629 else 630 exc(:)=exc(:)+alpha*grho2_updn(:,3)/(rho_updn(:,1)+rho_updn(:,2)) 631 if(nvxcgrho==2)vxcgrho(:,1:2)=2.0d0*alpha 632 if(nvxcgrho==3)vxcgrho(:,3)=2.0d0*alpha 633 end if 634 case (34) 635 alpha=-0.010d0 636 ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc) 637 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 638 if (nspden==1) then 639 exc(:)=exc(:)+16.0d0*alpha*tau_updn(:,1) 640 vxcrho(:,1)=vxcrho(:,1)+16.0d0*alpha*tau_updn(:,1) 641 vxctau(:,1)=16.0d0*alpha*rho_updn(:,1) 642 else 643 do ispden=1,nspden 644 exc(:)=exc(:)+8.0d0*alpha*tau_updn(:,ispden) 645 vxcrho(:,ispden)=vxcrho(:,ispden)+8.0d0*alpha*(tau_updn(:,1)+tau_updn(:,2)) 646 vxctau(:,ispden)=8.0d0*alpha*(rho_updn(:,1)+rho_updn(:,2)) 647 end do 648 end if 649 end select 650 651 !>>>>> Hybrid PBE0 (1/4 and 1/3) 652 else if(ixc>=41.and.ixc<=42) then 653 ! Requires to evaluate exchange-correlation with PBE (optpbe=2) 654 ! minus hyb_mixing*exchange with PBE (optpbe=-2) 655 ndvxc_x=8 656 ABI_ALLOCATE(exc_x,(npts)) 657 ABI_ALLOCATE(vxcrho_x,(npts,nspden)) 658 ABI_ALLOCATE(vxcgrho_x,(npts,nvxcgrho)) 659 exc_x=zero;vxcrho_x=zero;vxcgrho_x=zero 660 if (order**2 <= 1) then 661 optpbe=2 !PBE exchange correlation 662 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 663 & dvxcdgr=vxcgrho,exexch=exexch_,grho2_updn=grho2_updn) 664 optpbe=-2 !PBE exchange-only 665 call xcpbe(exc_x,npts,nspden,optpbe,order,rho_updn,vxcrho_x,ndvxc,ngr2,nd2vxc,& 666 & dvxcdgr=vxcgrho_x,exexch=exexch_,grho2_updn=grho2_updn) 667 exc=exc-exc_x*hyb_mixing 668 vxcrho=vxcrho-vxcrho_x*hyb_mixing 669 vxcgrho=vxcgrho-vxcgrho_x*hyb_mixing 670 else if (order /=3) then 671 ABI_ALLOCATE(dvxc_x,(npts,ndvxc_x)) 672 optpbe=2 !PBE exchange correlation 673 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 674 dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 675 optpbe=-2 !PBE exchange-only 676 call xcpbe(exc_x,npts,nspden,optpbe,order,rho_updn,vxcrho_x,ndvxc_x,ngr2,nd2vxc,& 677 & dvxcdgr=vxcgrho_x,dvxci=dvxc_x,grho2_updn=grho2_updn) 678 exc=exc-exc_x*hyb_mixing 679 vxcrho=vxcrho-vxcrho_x*hyb_mixing 680 vxcgrho=vxcgrho-vxcgrho_x*hyb_mixing 681 dvxc(:,1:ndvxc_x)=dvxc(:,1:ndvxc_x)-dvxc_x(:,1:ndvxc_x)*hyb_mixing 682 ABI_DEALLOCATE(dvxc_x) 683 else if (order ==3) then 684 ! The size of exchange-correlation with PBE (optpbe=2) 685 ! is the one which defines the size for ndvxc. 686 ABI_ALLOCATE(dvxc_x,(npts,ndvxc_x)) 687 ABI_ALLOCATE(d2vxc_x,(npts,nd2vxc)) 688 optpbe=2 !PBE exchange correlation 689 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 690 & d2vxci=d2vxc,dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 691 optpbe=-2 !PBE exchange-only 692 call xcpbe(exc_x,npts,nspden,optpbe,order,rho_updn,vxcrho_x,ndvxc_x,ngr2,nd2vxc,& 693 & d2vxci=d2vxc_x,dvxcdgr=vxcgrho_x,dvxci=dvxc_x,grho2_updn=grho2_updn) 694 exc=exc-exc_x*hyb_mixing 695 vxcrho=vxcrho-vxcrho_x*hyb_mixing 696 vxcgrho=vxcgrho-vxcgrho_x*hyb_mixing 697 d2vxc=d2vxc-d2vxc_x*hyb_mixing 698 dvxc(:,1:ndvxc_x)=dvxc(:,1:ndvxc_x)-dvxc_x(:,1:ndvxc_x)*hyb_mixing 699 ABI_DEALLOCATE(dvxc_x) 700 ABI_DEALLOCATE(d2vxc_x) 701 end if 702 ABI_DEALLOCATE(exc_x) 703 ABI_DEALLOCATE(vxcrho_x) 704 ABI_DEALLOCATE(vxcgrho_x) 705 706 !>>>>> Ichimaru,Iyetomi,Tanaka, XC at finite temp (e- gaz) 707 else if (ixc==50) then 708 if (order**2 <= 1) then 709 call xciit(exc,fxcT,npts,order,rspts,el_temp,vxcrho(:,1)) 710 else 711 call xciit(exc,fxcT,npts,order,rspts,el_temp,vxcrho(:,1),dvxc) 712 end if 713 714 !>>>>> GGA counterpart of the B3LYP functional 715 else if(ixc==1402000) then 716 ! Requires to evaluate exchange-correlation 717 ! with 5/4 B3LYP - 1/4 B3LYPc, where 718 ! B3LYPc = (0.19 Ec VWN3 + 0.81 Ec LYP) 719 720 ! First evaluate B3LYP. 721 if(present(xc_funcs))then 722 if (abs(order)==1) then 723 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 724 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,xc_functionals=xc_funcs) 725 elseif (abs(order)==2) then 726 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 727 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,xc_functionals=xc_funcs) 728 else 729 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 730 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,d2vxc=d2vxc,xc_functionals=xc_funcs) 731 end if 732 else 733 if (abs(order)==1) then 734 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 735 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho) 736 elseif (abs(order)==2) then 737 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 738 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc) 739 else 740 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 741 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,d2vxc=d2vxc) 742 end if 743 end if 744 745 ! Then renormalize B3LYP and subtract VWN3 contribution 746 ABI_ALLOCATE(exc_c,(npts)) 747 ABI_ALLOCATE(vxcrho_c,(npts,nspden)) 748 if(order**2>1)then 749 ABI_ALLOCATE(dvxc_c,(npts,ndvxc)) 750 end if 751 if(order**2>4)then 752 ABI_ALLOCATE(d2vxc_c,(npts,nd2vxc)) 753 end if 754 exc_c=zero;vxcrho_c=zero 755 call libxc_functionals_init(-30,nspden,xc_functionals=xc_funcs_vwn3) 756 if (order**2 <= 1) then 757 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 758 & vxcrho_c,xc_functionals=xc_funcs_vwn3) 759 elseif (order**2 <= 4) then 760 dvxc_c=zero 761 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 762 & vxcrho_c,dvxc=dvxc_c,xc_functionals=xc_funcs_vwn3) 763 else 764 dvxc_c=zero 765 d2vxc_c=zero 766 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 767 & vxcrho_c,dvxc=dvxc_c,d2vxc=d2vxc,xc_functionals=xc_funcs_vwn3) 768 end if 769 exc=1.25d0*exc-quarter*0.19d0*exc_c 770 vxcrho=1.25d0*vxcrho-quarter*0.19d0*vxcrho_c 771 if(order**2>1)dvxc=1.25d0*dvxc-quarter*0.19d0*dvxc_c 772 if(order**2>4)d2vxc=1.25d0*d2vxc-quarter*0.19d0*d2vxc_c 773 call libxc_functionals_end(xc_functionals=xc_funcs_vwn3) 774 775 ! Then subtract LYP contribution 776 call libxc_functionals_init(-131,nspden,xc_functionals=xc_funcs_lyp) 777 if (order**2 <= 1) then 778 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 779 & vxcrho_c,grho2=grho2_updn,vxcgr=vxcgrho,xc_functionals=xc_funcs_lyp) 780 elseif (order**2 <= 4) then 781 dvxc_c=zero 782 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 783 & vxcrho_c,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc_c,xc_functionals=xc_funcs_lyp) 784 else 785 dvxc_c=zero 786 d2vxc_c=zero 787 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 788 & vxcrho_c,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc_c,d2vxc=d2vxc,xc_functionals=xc_funcs_lyp) 789 end if 790 exc=exc-quarter*0.81d0*exc_c 791 vxcrho=vxcrho-quarter*0.81d0*vxcrho_c 792 if(order**2>1)dvxc=dvxc-quarter*0.81d0*dvxc_c 793 if(order**2>4)d2vxc=d2vxc-quarter*0.81d0*d2vxc_c 794 call libxc_functionals_end(xc_functionals=xc_funcs_lyp) 795 796 ABI_DEALLOCATE(exc_c) 797 ABI_DEALLOCATE(vxcrho_c) 798 if(allocated(dvxc_c))then 799 ABI_DEALLOCATE(dvxc_c) 800 end if 801 if(allocated(d2vxc_c))then 802 ABI_DEALLOCATE(d2vxc_c) 803 end if 804 805 !>>>>> All libXC functionals 806 else if( ixc<0 ) then 807 if (is_mgga) then 808 if(present(xc_funcs))then 809 if (PRESENT(xc_tb09_c)) then 810 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 811 & vxcrho,grho2_updn,vxcgrho,lrho_updn,vxclrho,tau_updn,vxctau,xc_functionals=xc_funcs,xc_tb09_c=xc_tb09_c) 812 else 813 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 814 & vxcrho,grho2_updn,vxcgrho,lrho_updn,vxclrho,tau_updn,vxctau,xc_functionals=xc_funcs) 815 end if 816 else 817 if (PRESENT(xc_tb09_c)) then 818 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 819 & vxcrho,grho2_updn,vxcgrho,lrho_updn,vxclrho,tau_updn,vxctau,xc_tb09_c=xc_tb09_c) 820 else 821 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 822 & vxcrho,grho2_updn,vxcgrho,lrho_updn,vxclrho,tau_updn,vxctau) 823 end if 824 end if 825 ixc1 = (-ixc)/1000 826 ixc2 = (-ixc) - ixc1*1000 827 if(ixc1==206 .or. ixc1==207 .or. ixc1==208 .or. ixc1==209 .or. & 828 & ixc2==206 .or. ixc2==207 .or. ixc2==208 .or. ixc2==209 )then 829 ! Assume that that type of mGGA can only be used with a LDA correlation (see doc) 830 vxcgrho(:,:)=zero 831 vxclrho(:,:)=zero 832 vxctau(:,:)=zero 833 end if 834 elseif (is_gga) then 835 if(present(xc_funcs))then 836 if (order**2 <= 1) then 837 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 838 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,xc_functionals=xc_funcs) 839 else 840 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 841 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,xc_functionals=xc_funcs) 842 end if 843 else 844 if (order**2 <= 1) then 845 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 846 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho) 847 else 848 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 849 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc) 850 end if 851 end if 852 else 853 if(present(xc_funcs))then 854 if (order**2 <= 1) then 855 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 856 & vxcrho,xc_functionals=xc_funcs) 857 elseif (order**2 <= 4) then 858 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 859 & vxcrho,dvxc=dvxc,xc_functionals=xc_funcs) 860 else 861 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 862 & vxcrho,dvxc=dvxc,d2vxc=d2vxc,xc_functionals=xc_funcs) 863 end if 864 else 865 if (order**2 <= 1) then 866 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 867 & vxcrho) 868 elseif (order**2 <= 4) then 869 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 870 & vxcrho,dvxc=dvxc) 871 else 872 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 873 & vxcrho,dvxc=dvxc,d2vxc=d2vxc) 874 end if 875 end if 876 end if 877 878 end if 879 880 ! ================================================= 881 ! == Finalization == 882 ! ================================================= 883 !Deallocate arrays 884 if(allocated(rhotot)) then 885 ABI_DEALLOCATE(rhotot) 886 end if 887 if(allocated(rspts)) then 888 ABI_DEALLOCATE(rspts) 889 end if 890 if(allocated(zeta)) then 891 ABI_DEALLOCATE(zeta) 892 end if 893 894 end subroutine drivexc