1 2! Copyright (C) 2009 T. McQueen and J. K. Dewhurst. 3! This file is distributed under the terms of the GNU General Public License. 4! See the file COPYING for license details. 5 6module libxcifc 7 8use xc_f90_lib_m 9 10! libxc version number 11integer libxcv(3) 12 13contains 14 15!BOP 16! !ROUTINE: xcifc_libxc 17! !INTERFACE: 18subroutine xcifc_libxc(xctype,n,c_tb09,tempa,rho,rhoup,rhodn,g2rho,g2up,g2dn, & 19 grho2,gup2,gdn2,gupdn,tau,tauup,taudn,ex,ec,vx,vc,vxup,vxdn,vcup,vcdn,dxdgr2, & 20 dxdgu2,dxdgd2,dxdgud,dcdgr2,dcdgu2,dcdgd2,dcdgud,dxdg2r,dxdg2u,dxdg2d,dcdg2r, & 21 dcdg2u,dcdg2d,wx,wxup,wxdn,wc,wcup,wcdn) 22! !USES: 23use, intrinsic :: iso_c_binding 24! !INPUT/OUTPUT PARAMETERS: 25! xctype : type of exchange-correlation functional (in,integer(3)) 26! n : number of density points (in,integer) 27! c_tb09 : Tran-Blaha '09 constant c (in,real,optional) 28! tempa : temperature in atomic units (in,real,optional) 29! rho : spin-unpolarised charge density (in,real(n),optional) 30! rhoup : spin-up charge density (in,real(n),optional) 31! rhodn : spin-down charge density (in,real(n),optional) 32! g2rho : grad^2 rho (in,real(n),optional) 33! g2up : grad^2 rhoup (in,real(n),optional) 34! g2dn : grad^2 rhodn (in,real(n),optional) 35! grho2 : |grad rho|^2 (in,real(n),optional) 36! gup2 : |grad rhoup|^2 (in,real(n),optional) 37! gdn2 : |grad rhodn|^2 (in,real(n),optional) 38! gupdn : (grad rhoup).(grad rhodn) (in,real(n),optional) 39! tau : kinetic energy density (in,real(n),optional) 40! tauup : spin-up kinetic energy density (in,real(n),optional) 41! taudn : spin-down kinetic energy density (in,real(n),optional) 42! ex : exchange energy density (out,real(n),optional) 43! ec : correlation energy density (out,real(n),optional) 44! vx : spin-unpolarised exchange potential (out,real(n),optional) 45! vc : spin-unpolarised correlation potential (out,real(n),optional) 46! vxup : spin-up exchange potential (out,real(n),optional) 47! vxdn : spin-down exchange potential (out,real(n),optional) 48! vcup : spin-up correlation potential (out,real(n),optional) 49! vcdn : spin-down correlation potential (out,real(n),optional) 50! dxdgr2 : de_x/d(|grad rho|^2) (out,real(n),optional) 51! dxdgu2 : de_x/d(|grad rhoup|^2) (out,real(n),optional) 52! dxdgd2 : de_x/d(|grad rhodn|^2) (out,real(n),optional) 53! dxdgud : de_x/d((grad rhoup).(grad rhodn)) (out,real(n),optional) 54! dcdgr2 : de_c/d(|grad rho|^2) (out,real(n),optional) 55! dcdgu2 : de_c/d(|grad rhoup|^2) (out,real(n),optional) 56! dcdgd2 : de_c/d(|grad rhodn|^2) (out,real(n),optional) 57! dcdgud : de_c/d((grad rhoup).(grad rhodn)) (out,real(n),optional) 58! dxdg2r : de_x/d(grad^2 rho) (out,real(n),optional) 59! dxdg2u : de_x/d(grad^2 rhoup) (out,real(n),optional) 60! dxdg2d : de_x/d(grad^2 rhodn) (out,real(n),optional) 61! dcdg2r : de_c/d(grad^2 rho) (out,real(n),optional) 62! dcdg2u : de_c/d(grad^2 rhoup) (out,real(n),optional) 63! dcdg2d : de_c/d(grad^2 rhodn) (out,real(n),optional) 64! wx : de_x/dtau (out,real(n),optional) 65! wxup : de_x/dtauup (out,real(n),optional) 66! wxdn : de_x/dtaudn (out,real(n),optional) 67! wc : de_c/dtau (out,real(n),optional) 68! wcup : de_c/dtauup (out,real(n),optional) 69! wcdn : de_c/dtaudn (out,real(n),optional) 70! !DESCRIPTION: 71! Interface to the ETSF {\tt libxc} exchange-correlation functional library: 72! \newline{\tt http://www.tddft.org/programs/octopus/wiki/index.php/Libxc}. 73! The second and third integers in {\tt xctype} define the exchange and 74! correlation functionals in {\tt libxc}, respectively. 75! 76! !REVISION HISTORY: 77! Created April 2009 (Tyrel McQueen) 78! Modified September 2009 (JKD and TMQ) 79! Updated for Libxc 1, July 2010 (JKD) 80! Updated for Libxc 4, March 2018 (JKD) 81! Updated for Libxc 5, May 2020 (JKD) 82!EOP 83!BOC 84implicit none 85! mandatory arguments 86integer, intent(in) :: xctype(3),n 87! optional arguments 88real(8), optional, intent(in) :: c_tb09,tempa 89real(8), optional, intent(in) :: rho(n),rhoup(n),rhodn(n) 90real(8), optional, intent(in) :: g2rho(n),g2up(n),g2dn(n) 91real(8), optional, intent(in) :: grho2(n),gup2(n),gdn2(n),gupdn(n) 92real(8), optional, intent(in) :: tau(n),tauup(n),taudn(n) 93real(8), optional, intent(out) :: ex(n),ec(n),vx(n),vc(n) 94real(8), optional, intent(out) :: vxup(n),vxdn(n),vcup(n),vcdn(n) 95real(8), optional, intent(out) :: dxdgr2(n),dxdgu2(n),dxdgd2(n),dxdgud(n) 96real(8), optional, intent(out) :: dxdg2r(n),dxdg2u(n),dxdg2d(n) 97real(8), optional, intent(out) :: wx(n),wxup(n),wxdn(n) 98real(8), optional, intent(out) :: dcdgr2(n),dcdgu2(n),dcdgd2(n),dcdgud(n) 99real(8), optional, intent(out) :: dcdg2r(n),dcdg2u(n),dcdg2d(n) 100real(8), optional, intent(out) :: wc(n),wcup(n),wcdn(n) 101! local variables 102integer nspin,fmly,id,k 103integer(c_size_t) np 104type(xc_f90_func_t) p 105real(8) ta 106! allocatable arrays 107real(8), allocatable :: r(:,:),sigma(:,:),vrho(:,:),vsigma(:,:) 108real(8), allocatable :: lapl(:,:),t(:,:),vlapl(:,:),vtau(:,:) 109if (present(rho)) then 110 nspin=XC_UNPOLARIZED 111else if (present(rhoup).and.present(rhodn)) then 112 nspin=XC_POLARIZED 113else 114 write(*,*) 115 write(*,'("Error(xcifc_libxc): missing arguments")') 116 write(*,*) 117 stop 118end if 119if (xctype(2).ne.0) then 120 if (xctype(2).eq.xctype(3)) then 121 write(*,*) 122 write(*,'("Error(xcifc_libxc): Libxc exchange and correlation & 123 &functionals")') 124 write(*,'(" are the same : ",2I8)') xctype(2:3) 125 write(*,*) 126 stop 127 end if 128end if 129! convert number of points to long integer 130np=n 131! loop over functional kinds (exchange or correlation) 132do k=2,3 133 id=xctype(k) 134 if (id.gt.0) then 135 fmly=xc_f90_family_from_id(id) 136! initialise functional 137 call xc_f90_func_init(p,id,nspin) 138 select case(fmly) 139 case(XC_FAMILY_LDA) 140!-------------------------! 141! LDA functionals ! 142!-------------------------! 143! set temperature for free energy functional 144 if ((id.eq.XC_LDA_XC_KSDT).or.(id.eq.XC_LDA_XC_GDSMFB)) then 145 call xc_f90_func_set_ext_params(p,[tempa]) 146 end if 147 if (k.eq.2) then 148! exchange 149 if (present(rho)) then 150 call xc_f90_lda_exc_vxc(p,np,rho(1),ex(1),vx(1)) 151 else 152 allocate(r(2,n),vrho(2,n)) 153 r(1,:)=rhoup(:); r(2,:)=rhodn(:) 154 call xc_f90_lda_exc_vxc(p,np,r(1,1),ex(1),vrho(1,1)) 155 vxup(:)=vrho(1,:); vxdn(:)=vrho(2,:) 156 deallocate(r,vrho) 157 end if 158 else 159! correlation 160 if (present(rho)) then 161 call xc_f90_lda_exc_vxc(p,np,rho(1),ec(1),vc(1)) 162 else 163 allocate(r(2,n),vrho(2,n)) 164 r(1,:)=rhoup(:); r(2,:)=rhodn(:) 165 call xc_f90_lda_exc_vxc(p,np,r(1,1),ec(1),vrho(1,1)) 166 vcup(:)=vrho(1,:); vcdn=vrho(2,:) 167 deallocate(r,vrho) 168 end if 169 end if 170 case(XC_FAMILY_GGA,XC_FAMILY_HYB_GGA) 171!-------------------------! 172! GGA functionals ! 173!-------------------------! 174 if (k.eq.2) then 175! exchange 176 if (present(rho)) then 177 call xc_f90_gga_exc_vxc(p,np,rho(1),grho2(1),ex(1),vx(1),dxdgr2(1)) 178 else 179 allocate(r(2,n),sigma(3,n),vrho(2,n),vsigma(3,n)) 180 r(1,:)=rhoup(:); r(2,:)=rhodn(:) 181 sigma(1,:)=gup2(:); sigma(2,:)=gupdn(:); sigma(3,:)=gdn2(:) 182 call xc_f90_gga_exc_vxc(p,np,r(1,1),sigma(1,1),ex(1),vrho(1,1), & 183 vsigma(1,1)) 184 vxup(:)=vrho(1,:); vxdn(:)=vrho(2,:) 185 dxdgu2(:)=vsigma(1,:); dxdgud(:)=vsigma(2,:); dxdgd2(:)=vsigma(3,:) 186 deallocate(r,sigma,vrho,vsigma) 187 end if 188 else 189! correlation 190 if (present(rho)) then 191 call xc_f90_gga_exc_vxc(p,np,rho(1),grho2(1),ec(1),vc(1),dcdgr2(1)) 192 else 193 allocate(r(2,n),sigma(3,n),vrho(2,n),vsigma(3,n)) 194 r(1,:)=rhoup(:); r(2,:)=rhodn(:) 195 sigma(1,:)=gup2(:); sigma(2,:)=gupdn(:); sigma(3,:)=gdn2(:) 196 call xc_f90_gga_exc_vxc(p,np,r(1,1),sigma(1,1),ec(1),vrho(1,1), & 197 vsigma(1,1)) 198 vcup(:)=vrho(1,:); vcdn(:)=vrho(2,:) 199 dcdgu2(:)=vsigma(1,:); dcdgud(:)=vsigma(2,:); dcdgd2(:)=vsigma(3,:) 200 deallocate(r,sigma,vrho,vsigma) 201 end if 202 end if 203 case(XC_FAMILY_MGGA) 204!------------------------------! 205! meta-GGA functionals ! 206!------------------------------! 207! set Tran-Blaha '09 constant if required 208 if (id.eq.XC_MGGA_X_TB09) then 209 if (present(c_tb09)) call xc_f90_func_set_ext_params(p,[c_tb09]) 210 end if 211 if (k.eq.2) then 212! exchange 213 if (present(rho)) then 214 if (present(ex)) then 215! spin-unpolarised energy functional 216 call xc_f90_mgga_exc_vxc(p,np,rho(1),grho2(1),g2rho(1),tau(1), & 217 ex(1),vx(1),dxdgr2(1),dxdg2r(1),wx(1)) 218 else 219! spin-unpolarised potential-only functional 220 allocate(vsigma(1,n),vlapl(1,n),vtau(1,n)) 221 call xc_f90_mgga_vxc(p,np,rho(1),grho2(1),g2rho(1),tau(1),vx(1), & 222 vsigma(1,1),vlapl(1,1),vtau(1,1)) 223 deallocate(vsigma,vlapl,vtau) 224 end if 225 else 226 allocate(r(2,n),sigma(3,n),lapl(2,n),t(2,n)) 227 allocate(vrho(2,n),vsigma(3,n),vlapl(2,n),vtau(2,n)) 228 r(1,:)=rhoup(:); r(2,:)=rhodn(:) 229 sigma(1,:)=gup2(:); sigma(2,:)=gupdn(:); sigma(3,:)=gdn2(:) 230 lapl(1,:)=g2up(:); lapl(2,:)=g2dn(:) 231 t(1,:)=tauup(:); t(2,:)=taudn(:) 232 if (present(ex)) then 233! spin-polarised energy functional 234 call xc_f90_mgga_exc_vxc(p,np,r(1,1),sigma(1,1),lapl(1,1),t(1,1), & 235 ex(1),vrho(1,1),vsigma(1,1),vlapl(1,1),vtau(1,1)) 236 dxdgu2(:)=vsigma(1,:); dxdgud(:)=vsigma(2,:); dxdgd2(:)=vsigma(3,:) 237 dxdg2u(:)=vlapl(1,:); dxdg2d(:)=vlapl(2,:) 238 wxup(:)=vtau(1,:); wxdn(:)=vtau(2,:) 239 else 240! spin-polarised potential-only functional 241 call xc_f90_mgga_vxc(p,np,r(1,1),sigma(1,1),lapl(1,1),t(1,1), & 242 vrho(1,1),vsigma(1,1),vlapl(1,1),vtau(1,1)) 243 end if 244 vxup(:)=vrho(1,:); vxdn(:)=vrho(2,:) 245 deallocate(r,sigma,lapl,t) 246 deallocate(vrho,vsigma,vlapl,vtau) 247 end if 248 else 249! correlation 250 if (present(rho)) then 251 if (present(ec)) then 252 call xc_f90_mgga_exc_vxc(p,np,rho(1),grho2(1),g2rho(1),tau(1), & 253 ec(1),vc(1),dcdgr2(1),dcdg2r(1),wc(1)) 254 else 255 allocate(vsigma(1,n),vlapl(1,n),vtau(1,n)) 256 call xc_f90_mgga_vxc(p,np,rho(1),grho2(1),g2rho(1),tau(1),vc(1), & 257 vsigma(1,1),vlapl(1,1),vtau(1,1)) 258 deallocate(vsigma,vlapl,vtau) 259 end if 260 else 261 allocate(r(2,n),sigma(3,n),lapl(2,n),t(2,n)) 262 allocate(vrho(2,n),vsigma(3,n),vlapl(2,n),vtau(2,n)) 263 r(1,:)=rhoup(:); r(2,:)=rhodn(:) 264 sigma(1,:)=gup2(:); sigma(2,:)=gupdn(:); sigma(3,:)=gdn2(:) 265 lapl(1,:)=g2up(:); lapl(2,:)=g2dn(:) 266 t(1,:)=tauup(:); t(2,:)=taudn(:) 267 if (present(ec)) then 268 call xc_f90_mgga_exc_vxc(p,np,r(1,1),sigma(1,1),lapl(1,1),t(1,1), & 269 ec(1),vrho(1,1),vsigma(1,1),vlapl(1,1),vtau(1,1)) 270 dcdgu2(:)=vsigma(1,:); dcdgud(:)=vsigma(2,:); dcdgd2(:)=vsigma(3,:) 271 dcdg2u(:)=vlapl(1,:); dcdg2d(:)=vlapl(2,:) 272 wcup(:)=vtau(1,:); wcdn(:)=vtau(2,:) 273 else 274 call xc_f90_mgga_vxc(p,np,r(1,1),sigma(1,1),lapl(1,1),t(1,1), & 275 vrho(1,1),vsigma(1,1),vlapl(1,1),vtau(1,1)) 276 end if 277 vcup(:)=vrho(1,:); vcdn(:)=vrho(2,:) 278 deallocate(r,sigma,lapl,t) 279 deallocate(vrho,vsigma,vlapl,vtau) 280 end if 281 end if 282 case default 283 write(*,*) 284 write(*,'("Error(xcifc_libxc): unsupported Libxc functional family : ",& 285 &I8)') fmly 286 write(*,*) 287 stop 288 end select 289! destroy functional 290 call xc_f90_func_end(p) 291 else 292! case when id=0 293 if (k.eq.2) then 294 if (present(ex)) ex(:)=0.d0 295 if (present(vx)) vx(:)=0.d0 296 if (present(vxup)) vxup(:)=0.d0 297 if (present(vxdn)) vxdn(:)=0.d0 298 if (present(dxdgr2)) dxdgr2(:)=0.d0 299 if (present(dxdgu2)) dxdgu2(:)=0.d0 300 if (present(dxdgd2)) dxdgd2(:)=0.d0 301 if (present(dxdgud)) dxdgud(:)=0.d0 302 else 303 if (present(ec)) ec(:)=0.d0 304 if (present(vc)) vc(:)=0.d0 305 if (present(vcup)) vcup(:)=0.d0 306 if (present(vcdn)) vcdn(:)=0.d0 307 if (present(dcdgr2)) dcdgr2(:)=0.d0 308 if (present(dcdgu2)) dcdgu2(:)=0.d0 309 if (present(dcdgd2)) dcdgd2(:)=0.d0 310 if (present(dcdgud)) dcdgud(:)=0.d0 311 end if 312 end if 313end do 314end subroutine 315 316subroutine fxcifc_libxc(fxctype,n,rho,rhoup,rhodn,fxc,fxcuu,fxcud,fxcdd) 317use, intrinsic :: iso_c_binding 318implicit none 319! mandatory arguments 320integer, intent(in) :: fxctype(3),n 321! optional arguments 322real(8), optional, intent(in) :: rho(n),rhoup(n),rhodn(n) 323real(8), optional, intent(out) :: fxc(n),fxcuu(n),fxcud(n),fxcdd(n) 324! local variables 325integer nspin,fmly,id,k 326integer(c_size_t) np 327type(xc_f90_func_t) p 328! allocatable arrays 329real(8), allocatable :: r(:,:),f(:,:) 330np=n 331if (present(rho)) then 332 nspin=XC_UNPOLARIZED 333else if (present(rhoup).and.present(rhodn)) then 334 nspin=XC_POLARIZED 335else 336 write(*,*) 337 write(*,'("Error(fxcifc_libxc): missing arguments")') 338 write(*,*) 339 stop 340end if 341! zero the kernel 342if (present(fxc)) fxc(:)=0.d0 343if (present(fxcuu)) fxcuu(:)=0.d0 344if (present(fxcud)) fxcud(:)=0.d0 345if (present(fxcdd)) fxcdd(:)=0.d0 346! loop over functional kinds (exchange or correlation) 347do k=2,3 348 id=fxctype(k) 349 if (id.le.0) cycle 350 fmly=xc_f90_family_from_id(id) 351! initialise functional 352 call xc_f90_func_init(p,id,nspin) 353 select case(fmly) 354 case(XC_FAMILY_LDA) 355!-------------------------! 356! LDA functionals ! 357!-------------------------! 358 if (present(rho)) then 359 allocate(f(1,n)) 360 call xc_f90_lda_fxc(p,np,rho(1),f(1,1)) 361 fxc(:)=fxc(:)+f(1,:) 362 deallocate(f) 363 else 364 allocate(r(2,n),f(3,n)) 365 r(1,:)=rhoup(:); r(2,:)=rhodn(:) 366 call xc_f90_lda_fxc(p,np,r(1,1),f(1,1)) 367 fxcuu(:)=fxcuu(:)+f(1,:) 368 fxcud(:)=fxcud(:)+f(2,:) 369 fxcdd(:)=fxcdd(:)+f(3,:) 370 deallocate(r,f) 371 end if 372 case default 373 write(*,*) 374 write(*,'("Error(fxcifc_libxc): unsupported Libxc functional family : ",& 375 &I8)') fmly 376 write(*,*) 377 stop 378 end select 379! destroy functional 380 call xc_f90_func_end(p) 381end do 382end subroutine 383 384subroutine xcdata_libxc(xctype,xcdescr,xcspin,xcgrad,hybrid,hybridc) 385implicit none 386! arguments 387integer, intent(in) :: xctype(3) 388character(512), intent(out) :: xcdescr 389integer, intent(out) :: xcspin 390integer, intent(out) :: xcgrad 391logical, intent(inout) :: hybrid 392real(8), intent(inout) :: hybridc 393! local variables 394integer j,k,id 395integer fmly,flg 396character(128) name 397type(xc_f90_func_t) p 398type(xc_f90_func_info_t) info 399! check version is compatible 400call xc_f90_version(libxcv(1),libxcv(2),libxcv(3)) 401if (libxcv(1).ne.5) then 402 write(*,*) 403 write(*,'("Error(xcdata_libxc): incompatible Libxc version : ",I2.2,".",& 404 &I3.3,".",I3.3)') libxcv(:) 405 write(*,*) 406 stop 407end if 408! unknown spin polarisation 409xcspin=-1 410! no gradients by default 411xcgrad=0 412! not hybrid by default 413hybrid=.false. 414do k=2,3 415 id=xctype(k) 416 if (id.gt.0) then 417 call xc_f90_func_init(p,id,XC_UNPOLARIZED) 418 name=xc_f90_functional_get_name(id) 419 fmly=xc_f90_family_from_id(id) 420 info=xc_f90_func_get_info(p) 421 if (fmly.eq.XC_FAMILY_HYB_GGA) then 422 hybridc=xc_f90_hyb_exx_coef(p) 423 hybrid=.true. 424 end if 425! hybrids should have only correlation part to avoid double-scaling with hybridc 426 if ((fmly.eq.XC_FAMILY_HYB_GGA).and.(k.eq.2)) then 427 write(*,*) 428 write(*,'("Error(xcdata_libxc): set only correlation part of xctype for & 429 &Libxc hybrids")') 430 write(*,*) 431 stop 432 end if 433! functional family 434 if ((fmly.ne.XC_FAMILY_LDA).and.(fmly.ne.XC_FAMILY_GGA).and. & 435 (fmly.ne.XC_FAMILY_HYB_GGA).and.(fmly.ne.XC_FAMILY_MGGA)) then 436 write(*,*) 437 write(*,'("Error(xcdata_libxc): unsupported Libxc family : ",I8)') fmly 438 write(*,*) 439 stop 440 end if 441! post-processed gradients required for GGA functionals 442 if (fmly.eq.XC_FAMILY_GGA.or.fmly.eq.XC_FAMILY_HYB_GGA) xcgrad=2 443! kinetic energy density required for meta-GGA functionals 444 if (fmly.eq.XC_FAMILY_MGGA) then 445 flg=xc_f90_func_info_get_flags(info) 446 if ((iand(flg,XC_FLAGS_HAVE_VXC).ne.0).and. & 447 (iand(flg,XC_FLAGS_HAVE_EXC).eq.0)) then 448! potential-only functional 449 xcgrad=3 450 else if (iand(flg,XC_FLAGS_HAVE_EXC).ne.0) then 451! energy functional 452 xcgrad=4 453 else 454 write(*,*) 455 write(*,'("Error(xcdata_libxc): unsupported Libxc meta-GGA type")') 456 write(*,*) 457 stop 458 end if 459! check if the density laplacian is required 460 if ((xcgrad.eq.4).and.(iand(flg,XC_FLAGS_NEEDS_LAPLACIAN).ne.0)) then 461 write(*,*) 462 write(*,'("Error(xcdata_libxc): energy meta-GGAs requiring the density & 463 &laplacian are not supported")') 464 write(*,*) 465 stop 466 end if 467 end if 468! destroy functional 469 call xc_f90_func_end(p) 470 else 471 name='none' 472 end if 473 if (k.eq.2) then 474 xcdescr='exchange: '//trim(name) 475 else 476 xcdescr=trim(xcdescr)//'; correlation: '//trim(name) 477 end if 478end do 479end subroutine 480!EOC 481 482end module 483 484