TABLE OF CONTENTS
ABINIT/drivexc_main [ Functions ]
NAME
drivexc_main
FUNCTION
Driver of XC functionals.Optionally, deliver the XC kernel, or even the derivative of the XC kernel (the third derivative of the XC energy)
COPYRIGHT
Copyright (C) 2012-2018 ABINIT group (MT) 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=index of the XC functional mgga=1 if functional is metaGGA ndvxc=size of dvxc(npts,ndvxc) nd2vxc=size of d2vxc(npts,nd2vxc) ngr2=size of grho2(npts,ngr2) npts=number of real space points on which the density (and its gradients) is provided nspden=number of spin-density components (1 or 2) nvxcgrho=size of vxcgrho(npts,nvxcgrho) 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(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 density5 and both are half the total density. If nspden=2, the spin-up and spin-down density must be given xclevel= XC functional level === Optional input arguments === [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(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(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(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. [xc_tb09_c]=c parameter for the mgga TB09 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) === Optional output arguments === [dvxc]=partial second derivatives of the xc energy, only if abs(order)>1 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))) [d2vxc]=second derivative of the XC potential=3rd order derivative of energy, only if abs(order)>2 (only available for LDA and nspden=1) 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) [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)
PARENTS
m_pawxc,rhotoxc
CHILDREN
drivexc
SOURCE
127 #if defined HAVE_CONFIG_H 128 #include "config.h" 129 #endif 130 131 #include "abi_common.h" 132 133 subroutine drivexc_main(exc,ixc,mgga,ndvxc,nd2vxc,ngr2,npts,nspden,nvxcgrho,order,rho,vxcrho,xclevel, & 134 & dvxc,d2vxc,el_temp,exexch,fxcT,grho2,& ! Optional arguments 135 & hyb_mixing,lrho,tau,vxcgrho,vxclrho,vxctau,xc_funcs,xc_tb09_c) ! Optional arguments 136 137 use defs_basis 138 use m_profiling_abi 139 use m_errors 140 use libxc_functionals 141 142 !This section has been created automatically by the script Abilint (TD). 143 !Do not modify the following lines by hand. 144 #undef ABI_FUNC 145 #define ABI_FUNC 'drivexc_main' 146 use interfaces_41_xc_lowlevel, except_this_one => drivexc_main 147 !End of the abilint section 148 149 implicit none 150 151 !Arguments ------------------------------------ 152 !scalars 153 integer,intent(in) :: ixc,mgga,ndvxc,nd2vxc,ngr2,npts,nspden,nvxcgrho,order,xclevel 154 integer,intent(in),optional :: exexch 155 real(dp),intent(in),optional :: el_temp,hyb_mixing,xc_tb09_c 156 !arrays 157 real(dp),intent(in) :: rho(npts,nspden) 158 real(dp),intent(in),optional :: grho2(npts,ngr2),lrho(npts,nspden*mgga),tau(npts,nspden*mgga) 159 real(dp),intent(out) :: exc(npts),vxcrho(npts,nspden) 160 real(dp),intent(out),optional :: dvxc(npts,ndvxc),d2vxc(npts,nd2vxc),fxcT(:),vxcgrho(npts,nvxcgrho) 161 real(dp),intent(out),optional :: vxclrho(npts,nspden*mgga),vxctau(npts,nspden*mgga) 162 type(libxc_functional_type),intent(inout),optional :: xc_funcs(2) 163 164 !Local variables------------------------------- 165 !scalars 166 real(dp) :: hyb_mixing_,xc_tb09_c_ 167 logical :: is_gga 168 169 ! ************************************************************************* 170 171 !Checks input parameters 172 if (mgga==1) then 173 if (.not.present(lrho)) then 174 MSG_BUG('lrho arg must be present in case of mGGA!') 175 end if 176 if (.not.present(tau)) then 177 MSG_BUG('tau arg must be present in case of mGGA!') 178 end if 179 if (.not.present(vxclrho)) then 180 MSG_BUG('vxclrho arg must be present in case of mGGA!') 181 end if 182 if (.not.present(vxctau)) then 183 MSG_BUG('vxctau arg must be present in case of mGGA!') 184 end if 185 end if 186 if (present(fxcT)) then 187 if (.not.present(el_temp)) then 188 MSG_BUG('el_temp arg must be present together with fxcT!') 189 end if 190 end if 191 192 xc_tb09_c_=99.99_dp;if (present(xc_tb09_c)) xc_tb09_c_=xc_tb09_c 193 194 if(ixc==41)hyb_mixing_=quarter 195 if(ixc==42)hyb_mixing_=third 196 if (present(hyb_mixing)) hyb_mixing_=hyb_mixing 197 198 !>>>>> All libXC functionals 199 200 if (ixc<0) then 201 202 if (present(xc_funcs))then 203 is_gga=libxc_functionals_isgga(xc_functionals=xc_funcs) 204 else 205 is_gga=libxc_functionals_isgga() 206 end if 207 208 if (mgga==1) then 209 if (abs(xc_tb09_c_-99.99_dp)>tol12) then 210 if (present(xc_funcs)) then 211 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 212 & lrho_updn=lrho,vxclrho=vxclrho,tau_updn=tau,vxctau=vxctau, & 213 & xc_funcs=xc_funcs,xc_tb09_c=xc_tb09_c_) 214 else 215 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 216 & grho2_updn=grho2,vxcgrho=vxcgrho, & 217 & lrho_updn=lrho,vxclrho=vxclrho,tau_updn=tau,vxctau=vxctau, & 218 & xc_tb09_c=xc_tb09_c_) 219 end if 220 else 221 if (present(xc_funcs)) then 222 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 223 & grho2_updn=grho2,vxcgrho=vxcgrho, & 224 & lrho_updn=lrho,vxclrho=vxclrho,tau_updn=tau,vxctau=vxctau, & 225 & xc_funcs=xc_funcs) 226 else 227 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 228 & grho2_updn=grho2,vxcgrho=vxcgrho, & 229 & lrho_updn=lrho,vxclrho=vxclrho,tau_updn=tau,vxctau=vxctau) 230 end if 231 end if 232 else if (is_gga) then 233 if (order**2<=1) then 234 if (present(xc_funcs)) then 235 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 236 & grho2_updn=grho2,vxcgrho=vxcgrho,xc_funcs=xc_funcs) 237 else 238 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 239 & grho2_updn=grho2,vxcgrho=vxcgrho) 240 end if 241 else 242 if (present(xc_funcs)) then 243 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 244 & grho2_updn=grho2,vxcgrho=vxcgrho,dvxc=dvxc,xc_funcs=xc_funcs) 245 else 246 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 247 & grho2_updn=grho2,vxcgrho=vxcgrho,dvxc=dvxc) 248 end if 249 end if 250 else 251 if (order**2<=1) then 252 if (present(xc_funcs)) then 253 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 254 xc_funcs=xc_funcs) 255 else 256 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho) 257 end if 258 else if (order**2<=4) then 259 if (present(xc_funcs)) then 260 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 261 & dvxc=dvxc, xc_funcs=xc_funcs) 262 else 263 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 264 & dvxc=dvxc) 265 end if 266 else 267 if (present(xc_funcs)) then 268 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 269 & dvxc=dvxc,d2vxc=d2vxc, xc_funcs=xc_funcs) 270 else 271 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 272 & dvxc=dvxc,d2vxc=d2vxc) 273 end if 274 end if 275 end if 276 277 else 278 279 ! Cases with gradient 280 if (xclevel==2)then 281 if (order**2<=1.or.ixc==16.or.ixc==17.or.ixc==26.or.ixc==27) then 282 if (ixc/=13) then 283 if (present(exexch)) then 284 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 285 & hyb_mixing=hyb_mixing_,grho2_updn=grho2,vxcgrho=vxcgrho, & 286 & exexch=exexch) 287 else 288 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 289 & hyb_mixing=hyb_mixing_,grho2_updn=grho2,vxcgrho=vxcgrho) 290 end if 291 else 292 if (present(exexch)) then 293 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 294 & grho2_updn=grho2, & 295 & exexch=exexch) 296 else 297 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 298 & grho2_updn=grho2) 299 end if 300 end if 301 else if (order/=3) then 302 if (ixc/=13) then 303 if (present(exexch)) then 304 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 305 & hyb_mixing=hyb_mixing_,dvxc=dvxc,grho2_updn=grho2,vxcgrho=vxcgrho, & 306 & exexch=exexch) 307 else 308 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 309 & hyb_mixing=hyb_mixing_,dvxc=dvxc,grho2_updn=grho2,vxcgrho=vxcgrho) 310 end if 311 else 312 if (present(exexch)) then 313 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 314 & dvxc=dvxc,grho2_updn=grho2, & 315 & exexch=exexch) 316 else 317 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 318 & dvxc=dvxc,grho2_updn=grho2) 319 end if 320 end if 321 else if (order==3) then 322 if (ixc/=13) then 323 if (present(exexch)) then 324 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 325 & dvxc=dvxc,d2vxc=d2vxc,grho2_updn=grho2,vxcgrho=vxcgrho, & 326 & hyb_mixing=hyb_mixing_,exexch=exexch) 327 else 328 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 329 & hyb_mixing=hyb_mixing_,dvxc=dvxc,d2vxc=d2vxc,grho2_updn=grho2,vxcgrho=vxcgrho) 330 end if 331 else 332 if (present(exexch)) then 333 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 334 & dvxc=dvxc,d2vxc=d2vxc,grho2_updn=grho2, & 335 & exexch=exexch) 336 else 337 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 338 & dvxc=dvxc,d2vxc=d2vxc,grho2_updn=grho2) 339 end if 340 end if 341 end if 342 343 344 ! Cases without gradient 345 else 346 if (order**2<=1) then 347 if (ixc>=31.and.ixc<=34) then !fake mgga functionals for testing purpose only (based on LDA functional) 348 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 349 & grho2_updn=grho2,vxcgrho=vxcgrho, & 350 & lrho_updn=lrho,vxclrho=vxclrho,tau_updn=tau,vxctau=vxctau) 351 else 352 if (present(fxcT)) then 353 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho,fxcT=fxcT,el_temp=el_temp) 354 else 355 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho) 356 end if 357 end if 358 else if (order==3.and.(ixc==3.or.ixc>=7.and.ixc<=10)) then 359 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 360 & dvxc=dvxc,d2vxc=d2vxc) 361 else 362 if (present(fxcT)) then 363 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 364 & dvxc=dvxc,fxcT=fxcT,el_temp=el_temp) 365 else 366 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 367 & dvxc=dvxc) 368 end if 369 end if 370 end if 371 372 end if 373 374 end subroutine drivexc_main