1!! NAME 2!! libxc_mod 3!! 4!! FUNCTION 5!! This module contains routines using features from libXC library 6!! (http://www.tddft.org/programs/octopus/wiki/index.php/Libxc) 7!! 8!! COPYRIGHT 9!! Copyright (C) 2015-2020 ABINIT group (MT) - Adapted for ATOMPAW 10!! This file is distributed under the terms of the 11!! GNU General Public License. 12!! See http://www.gnu.org/copyleft/gpl.txt. 13 14#if defined HAVE_CONFIG_H 15#include "config.h" 16#endif 17 18module libxc_mod 19 20 use io_tools 21 use Tools 22 use globalmath 23 24!ISO C bindings are mandatory 25#ifdef HAVE_FC_ISO_C_BINDING 26 use iso_c_binding 27#endif 28 29 implicit none 30 31!!================================================================= 32!! CONSTANTS 33!!================================================================= 34 35#if defined HAVE_LIBXC 36 logical,parameter,public :: have_libxc=.true. 37#else 38 logical,parameter,public :: have_libxc=.false. 39#endif 40 41!Public constants (use libxc_functionals_constants_load to init them) 42 integer,public,save :: XC_FAMILY_UNKNOWN = -1 43 integer,public,save :: XC_FAMILY_LDA = 1 44 integer,public,save :: XC_FAMILY_GGA = 2 45 integer,public,save :: XC_FAMILY_MGGA = 4 46 integer,public,save :: XC_FAMILY_LCA = 8 47 integer,public,save :: XC_FAMILY_OEP = 16 48 integer,public,save :: XC_FAMILY_HYB_GGA = 32 49 integer,public,save :: XC_FAMILY_HYB_MGGA = 64 50 integer,public,save :: XC_FLAGS_HAVE_EXC = 1 51 integer,public,save :: XC_FLAGS_HAVE_VXC = 2 52 integer,public,save :: XC_FLAGS_HAVE_FXC = 4 53 integer,public,save :: XC_FLAGS_HAVE_KXC = 8 54 integer,public,save :: XC_FLAGS_HAVE_LXC = 16 55 integer,public,save :: XC_FLAGS_NEEDS_LAPLACIAN = 32768 56 integer,public,save :: XC_EXCHANGE = 0 57 integer,public,save :: XC_CORRELATION = 1 58 integer,public,save :: XC_EXCHANGE_CORRELATION = 2 59 integer,public,save :: XC_KINETIC = 3 60 integer,public,save :: XC_HYB_NONE = 0 61 integer,public,save :: XC_HYB_FOCK = 1 62 integer,public,save :: XC_HYB_PT2 = 2 63 integer,public,save :: XC_HYB_ERF_SR = 4 64 integer,public,save :: XC_HYB_YUKAWA_SR = 8 65 integer,public,save :: XC_HYB_GAUSSIAN_SR = 16 66 integer,public,save :: XC_HYB_SEMILOCAL = 0 67 integer,public,save :: XC_HYB_HYBRID = 1 68 integer,public,save :: XC_HYB_CAM = 2 69 integer,public,save :: XC_HYB_CAMY = 3 70 integer,public,save :: XC_HYB_CAMG = 4 71 integer,public,save :: XC_HYB_DOUBLE_HYBRID = 5 72 integer,public,save :: XC_HYB_MIXTURE = 32768 73 integer,public,save :: XC_SINGLE_PRECISION = 0 74 logical,private,save :: libxc_constants_initialized=.false. 75 76!!================================================================= 77!! FUNCTIONS 78!!================================================================= 79 80!Public functions 81 public :: libxc_init_func ! Initialize libXC functional(s) 82 public :: libxc_end_func ! Destroy libXC functional(s) 83 public :: libxc_print_func ! Print libXC functionnal(s) details 84 public :: libxc_getid_fromInput ! From XC input string, gives the libXC id(s) 85 public :: libxc_getid_fromName ! From a char. string, gives the libXC id 86 public :: libxc_getshortname ! Gives a short version of the libXC name (w/o XC_) 87 public :: libxc_getid ! From LibXC datastructure, gives the libXC id(s) 88 public :: libxc_family_from_id ! Retrieve family of a XC functional from its libXC id 89 public :: libxc_islda ! Return TRUE if the XC functional is LDA 90 public :: libxc_isgga ! Return TRUE if the XC functional is GGA 91 public :: libxc_ismgga ! Return TRUE if the XC functional is MGGA 92 public :: libxc_needs_laplacian ! Return TRUE if functional uses LAPLACIAN 93 public :: libxc_ishybrid ! Return TRUE if the XC functional is Hybrid 94 public :: libxc_is_hybrid_from_id ! Return TRUE if a XC functional is hybrid, from its id 95 public :: libxc_has_kxc ! Return TRUE if Kxc (3rd der) is available 96 public :: libxc_has_k3xc ! Return TRUE if K3xc (4th der) is available 97 public :: libxc_nspin ! The number of spin components 98 public :: libxc_getvxc ! Return XC potential and energy, from input density 99 100!Private functions 101 private :: libxc_check ! Check if the code has been compiled with libXC 102 private :: libxc_constants_load ! Load libXC constants from C headers 103#if defined HAVE_FC_ISO_C_BINDING 104 private :: xc_char_to_c ! Convert a string from Fortran to C 105 private :: xc_char_to_f ! Convert a string from C to Fortran 106#endif 107 108!!================================================================= 109!! GLOBAL VARIABLE FOR XC FUNCTIONAL 110!!================================================================= 111 112!XC functional public type 113 type,public :: libxc_functional_t 114 integer :: id ! identifier 115 integer :: family ! LDA, GGA, etc. 116 integer :: xckind ! EXCHANGE, CORRELATION, etc. 117 integer :: nspin ! # of spin components 118 integer :: abi_ixc ! Abinit IXC id for this functional 119 logical :: has_exc ! TRUE is exc is available for the functional 120 logical :: has_vxc ! TRUE is vxc is available for the functional 121 logical :: has_fxc ! TRUE is fxc is available for the functional 122 logical :: has_kxc ! TRUE is kxc is available for the functional 123 logical :: needs_laplacian ! TRUE is functional needs laplacian of density 124 logical :: is_hybrid ! TRUE is functional is a hybrid functional 125 real(8) :: hyb_mixing ! Hybrid functional: mixing factor of Fock contribution (default=0) 126 real(8) :: hyb_mixing_sr ! Hybrid functional: mixing factor of SR Fock contribution (default=0) 127 real(8) :: hyb_range ! Range (for separation) for a hybrid functional (default=0) 128#ifdef HAVE_FC_ISO_C_BINDING 129 type(C_PTR),pointer :: conf => null() ! C pointer to the functional itself 130#endif 131 end type libxc_functional_t 132 133!Private global XC functional 134 type(libxc_functional_t),target,save,private :: libxc_funcs(2) 135 136!!================================================================= 137!! INTERFACES for C BINDINGS 138!!================================================================= 139#ifdef HAVE_FC_ISO_C_BINDING 140 141 interface 142 integer(C_INT) function xc_func_init(xc_func,functional,nspin) bind(C) 143 use iso_c_binding, only : C_INT,C_PTR 144 type(C_PTR) :: xc_func 145 integer(C_INT),value :: functional,nspin 146 end function xc_func_init 147 end interface 148 149 interface 150 subroutine xc_func_end(xc_func) bind(C) 151 use iso_c_binding, only : C_PTR 152 type(C_PTR) :: xc_func 153 end subroutine xc_func_end 154 end interface 155 156 interface 157 integer(C_INT) function xc_functional_get_number(name) bind(C) 158 use iso_c_binding, only : C_INT,C_PTR 159 type(C_PTR),value :: name 160 end function xc_functional_get_number 161 end interface 162 163 interface 164 type(C_PTR) function xc_functional_get_name(number) bind(C) 165 use iso_c_binding, only : C_INT,C_PTR 166 integer(C_INT),value :: number 167 end function xc_functional_get_name 168 end interface 169 170 interface 171 integer(C_INT) function xc_family_from_id(id,family,number) bind(C) 172 use iso_c_binding, only : C_INT,C_PTR 173 integer(C_INT),value :: id 174 type(C_PTR),value :: family,number 175 end function xc_family_from_id 176 end interface 177 178 interface 179 subroutine xc_hyb_cam_coef(xc_func,omega,alpha,beta) bind(C) 180 use iso_c_binding, only : C_DOUBLE,C_PTR 181 type(C_PTR) :: xc_func 182 real(C_DOUBLE) :: omega,alpha,beta 183 end subroutine xc_hyb_cam_coef 184 end interface 185 186 interface 187 subroutine xc_get_lda(xc_func,np,rho,zk,vrho,v2rho2,v3rho3) bind(C) 188 use iso_c_binding, only : C_INT,C_PTR 189 type(C_PTR) :: xc_func 190 integer(C_INT),value :: np 191 type(C_PTR),value :: rho,zk,vrho,v2rho2,v3rho3 192 end subroutine xc_get_lda 193 end interface 194 195 interface 196 subroutine xc_get_gga(xc_func,np,rho,sigma,zk,vrho,vsigma,v2rho2,v2rhosigma,v2sigma2, & 197& v3rho3,v3rho2sigma,v3rhosigma2,v3sigma3) bind(C) 198 use iso_c_binding, only : C_INT,C_PTR 199 type(C_PTR) :: xc_func 200 integer(C_INT),value :: np 201 type(C_PTR),value :: rho,sigma,zk,vrho,vsigma,v2rho2,v2rhosigma,v2sigma2, & 202& v3rho3,v3rho2sigma,v3rhosigma2,v3sigma3 203 end subroutine xc_get_gga 204 end interface 205 206 interface 207 subroutine xc_get_mgga(xc_func,np,rho,sigma,lapl,tau,zk,vrho,vsigma,vlapl,vtau, & 208& v2rho2,v2rhosigma,v2rholapl,v2rhotau,v2sigma2,v2sigmalapl, & 209& v2sigmatau,v2lapl2,v2lapltau,v2tau2) bind(C) 210 use iso_c_binding, only : C_INT,C_PTR 211 integer(C_INT),value :: np 212 type(C_PTR),value :: rho,sigma,lapl,tau,zk,vrho,vsigma,vlapl,vtau, & 213& v2rho2,v2sigma2,v2lapl2,v2tau2,v2rhosigma,v2rholapl,v2rhotau, & 214& v2sigmalapl,v2sigmatau,v2lapltau 215 type(C_PTR) :: xc_func 216 end subroutine xc_get_mgga 217 end interface 218 219 interface 220 subroutine xc_func_set_params(xc_func,params,n_params) bind(C) 221 use iso_c_binding, only : C_INT,C_DOUBLE,C_PTR 222 type(C_PTR) :: xc_func 223 integer(C_INT),value :: n_params 224 real(C_DOUBLE) :: params(*) 225 end subroutine xc_func_set_params 226 end interface 227 228 interface 229 subroutine xc_func_set_density_threshold(xc_func,dens_threshold) bind(C) 230 use iso_c_binding, only : C_DOUBLE,C_PTR 231 type(C_PTR) :: xc_func 232 real(C_DOUBLE) :: dens_threshold 233 end subroutine xc_func_set_density_threshold 234 end interface 235 236 interface 237 integer(C_INT) function xc_func_is_hybrid_from_id(func_id) bind(C) 238 use iso_c_binding, only : C_INT 239 integer(C_INT),value :: func_id 240 end function xc_func_is_hybrid_from_id 241 end interface 242 243 interface 244 subroutine xc_get_singleprecision_constant(xc_cst_singleprecision) bind(C) 245 use iso_c_binding, only : C_INT 246 integer(C_INT) :: xc_cst_singleprecision 247 end subroutine xc_get_singleprecision_constant 248 end interface 249 250 interface 251 subroutine xc_get_family_constants(xc_cst_unknown,xc_cst_lda,xc_cst_gga,xc_cst_mgga, & 252& xc_cst_lca,xc_cst_oep,xc_cst_hyb_gga,xc_cst_hyb_mgga) & 253& bind(C) 254 use iso_c_binding, only : C_INT 255 integer(C_INT) :: xc_cst_unknown,xc_cst_lda,xc_cst_gga,xc_cst_mgga, & 256& xc_cst_lca,xc_cst_oep,xc_cst_hyb_gga,xc_cst_hyb_mgga 257 end subroutine xc_get_family_constants 258 end interface 259 260 interface 261 subroutine xc_get_flags_constants(xc_cst_flags_have_exc,xc_cst_flags_have_vxc, & 262& xc_cst_flags_have_fxc,xc_cst_flags_have_kxc, & 263& xc_cst_flags_have_lxc,xc_cst_flags_needs_laplacian) bind(C) 264 use iso_c_binding, only : C_INT 265 integer(C_INT) :: xc_cst_flags_have_exc,xc_cst_flags_have_vxc,xc_cst_flags_have_fxc, & 266& xc_cst_flags_have_kxc,xc_cst_flags_have_lxc,xc_cst_flags_needs_laplacian 267 end subroutine xc_get_flags_constants 268 end interface 269 270 interface 271 subroutine xc_get_kind_constants(xc_cst_exchange,xc_cst_correlation, & 272& xc_cst_exchange_correlation,xc_cst_kinetic) bind(C) 273 use iso_c_binding, only : C_INT 274 integer(C_INT) :: xc_cst_exchange,xc_cst_correlation, & 275& xc_cst_exchange_correlation,xc_cst_kinetic 276 end subroutine xc_get_kind_constants 277 end interface 278 279 interface 280 subroutine xc_get_hybrid_constants(xc_cst_hyb_none, & 281 xc_cst_hyb_fock,xc_cst_hyb_pt2,xc_cst_hyb_erf_sr,xc_cst_hyb_yukawa_sr, & 282 xc_cst_hyb_gaussian_sr,xc_cst_hyb_semilocal, xc_cst_hyb_hybrid,xc_cst_hyb_cam, & 283 xc_cst_hyb_camy,xc_cst_hyb_camg,xc_cst_hyb_double_hybrid, & 284 xc_cst_hyb_mixture) bind(C) 285 use iso_c_binding, only : C_INT 286 integer(C_INT) :: xc_cst_hyb_none, xc_cst_hyb_fock,xc_cst_hyb_pt2, xc_cst_hyb_erf_sr, & 287 xc_cst_hyb_yukawa_sr,xc_cst_hyb_gaussian_sr,xc_cst_hyb_semilocal, & 288 xc_cst_hyb_hybrid,xc_cst_hyb_cam,xc_cst_hyb_camy,xc_cst_hyb_camg, & 289 xc_cst_hyb_double_hybrid,xc_cst_hyb_mixture 290 end subroutine xc_get_hybrid_constants 291 end interface 292 293 interface 294 type(C_PTR) function xc_func_type_malloc() bind(C) 295 use iso_c_binding, only : C_PTR 296 end function xc_func_type_malloc 297 end interface 298 299 interface 300 subroutine xc_func_type_free(xc_func) bind(C) 301 use iso_c_binding, only : C_PTR 302 type(C_PTR) :: xc_func 303 end subroutine xc_func_type_free 304 end interface 305 306 interface 307 type(C_PTR) function xc_get_info_name(xc_func) bind(C) 308 use iso_c_binding, only : C_PTR 309 type(C_PTR) :: xc_func 310 end function xc_get_info_name 311 end interface 312 313 interface 314 type(C_PTR) function xc_get_info_refs(xc_func,iref) bind(C) 315 use iso_c_binding, only : C_INT,C_PTR 316 type(C_PTR) :: xc_func 317 integer(C_INT) :: iref 318 end function xc_get_info_refs 319 end interface 320 321 interface 322 integer(C_INT) function xc_get_info_flags(xc_func) bind(C) 323 use iso_c_binding, only : C_INT,C_PTR 324 type(C_PTR) :: xc_func 325 end function xc_get_info_flags 326 end interface 327 328 interface 329 integer(C_INT) function xc_get_info_kind(xc_func) bind(C) 330 use iso_c_binding, only : C_INT,C_PTR 331 type(C_PTR) :: xc_func 332 end function xc_get_info_kind 333 end interface 334 335#endif 336 337 private 338 339 CONTAINS 340 341 342!!================================================================= 343!! NAME 344!! libxc_getid_fromInput 345!! 346!! FUNCTION 347!! From a character string (as given in ATOMPAW input file), 348!! gives the libXC id(s) 349!! 350!! INPUTS 351!! xcname= string containing the name of a XC functional 352!! 353!! OUTPUT 354!! id(2)= id(s) of the libXC functional 355!! [xcname_short]= short version of the libXC name (optional) 356!! 357!!================================================================= 358 subroutine libxc_getid_fromInput(xcname,id,xcname_short) 359 360 implicit none 361 integer,intent(inout) :: id(2) 362 character*(*),intent(in) :: xcname 363 character*(*),intent(out),optional :: xcname_short 364 365#if defined HAVE_LIBXC 366 367!---- Local variables 368 integer :: ii,i_plus 369 character*50 :: xcstrg(2) 370 371!------------------------------------------------------------------ 372!---- Executable code 373 374 i_plus=index(xcname,'+') 375 if (i_plus<=0) then 376 xcstrg(1)=trim(xcname) 377 xcstrg(2)="" 378 else 379 xcstrg(1)=trim(xcname(1:i_plus-1)) 380 xcstrg(2)=trim(xcname(i_plus+1:)) 381 end if 382 383 do ii=1,2 384 id(ii)=-1 385 call uppercase(xcstrg(ii)) 386 387 id(ii)=libxc_getid_fromName(xcstrg(ii)) 388 389 if (id(ii)==-1.and.ii==2) exit 390 391 if (id(ii)==-1.and.xcstrg(ii)(1:6)=="LIBXC_") then 392 read(unit=xcstrg(ii)(7:),fmt=*,err=333,end=333) id(ii) 393333 continue 394 end if 395 396 if (id(ii)==-1) then 397 write(std_out,'(/,2x,a)') "Error in get_id_from_name:" 398 write(std_out,'(2x,3a)') " Unknown X, C or XC functionnal (", & 399& trim(xcstrg(ii)),") !" 400 stop 401 end if 402 end do 403 404 if (present(xcname_short)) then 405 xcname_short="" 406 if (id(1)>0.and.xcstrg(1)(1:3)=="XC_") xcname_short=trim(xcname_short) //trim(xcstrg(1)(4:)) 407 if (id(1)>0.and.xcstrg(1)(1:6)=="LIBXC_") xcname_short=trim(xcname_short) //trim(xcstrg(1)) 408 if (id(2)>0.and.xcstrg(2)(1:3)=="XC_") xcname_short=trim(xcname_short)//"+"//trim(xcstrg(2)(4:)) 409 if (id(2)>0.and.xcstrg(2)(1:6)=="LIBXC_") xcname_short=trim(xcname_short)//"+"//trim(xcstrg(2)) 410 if (trim(xcname_short)=="") xcname_short=xcname 411 end if 412 413#else 414 id(1:2)=-2 415 if (present(xcname_short)) xcname_short=xcname 416#endif 417 418 end subroutine libxc_getid_fromInput 419 420 421!!================================================================= 422!! NAME 423!! libxc_getid_fromName 424!! 425!! FUNCTION 426!! Return identifer of a XC functional from its name 427!! Return -1 if undefined 428!! 429!! INPUTS 430!! xcname= string containing the name of a XC functional 431!! 432!!================================================================= 433 function libxc_getid_fromName(xcname) 434 435 implicit none 436 integer :: libxc_getid_fromName 437 character(len=*),intent(in) :: xcname 438 439#if defined HAVE_LIBXC 440 441!---- Local variables 442#if defined HAVE_FC_ISO_C_BINDING 443 character(len=256) :: str 444 character(kind=C_CHAR,len=1),target :: name_c(len_trim(xcname)+1) 445 character(kind=C_CHAR,len=1),target :: name_c_xc(len_trim(xcname)-2) 446 type(C_PTR) :: name_c_ptr 447#endif 448 449!------------------------------------------------------------------ 450!---- Executable code 451 452#if defined HAVE_FC_ISO_C_BINDING 453 str=trim(xcname) 454 if (xcname(1:3)=="XC_".or.xcname(1:3)=="xc_") then 455 str=xcname(4:);name_c_xc=xc_char_to_c(str) 456 name_c_ptr=c_loc(name_c_xc) 457 else 458 name_c=xc_char_to_c(str) 459 name_c_ptr=c_loc(name_c) 460 end if 461 libxc_getid_fromName=int(xc_functional_get_number(name_c_ptr)) 462#endif 463 464#else 465 libxc_getid_fromName=-1 466#endif 467 468end function libxc_getid_fromName 469 470 471!!================================================================= 472!! NAME 473!! libxc_getid 474!! 475!! FUNCTION 476!! From LibXC datastructure, gives the libXC id(s) 477!! 478!! OUTPUT 479!! id(2)= id(s) of the XC functional 480!! 481!!================================================================= 482 subroutine libxc_getid(id) 483 484 implicit none 485 integer :: id(2) 486 487#if defined HAVE_LIBXC 488 489!---- Local variables 490 integer :: ii 491 492!------------------------------------------------------------------ 493!---- Executable code 494 495 do ii=1,2 496 id(ii)=libxc_funcs(ii)%id 497 end do 498 499#else 500 id(1:2)=-2 501#endif 502 503 end subroutine libxc_getid 504 505 506!!================================================================= 507!! NAME 508!! libxc_getshortname 509!! 510!! FUNCTION 511!! From a character string (given in input file), gives a short 512!! version of the libXC name (without XC_) 513!! 514!! INPUTS 515!! xcname= string containing the name of a XC functional 516!! 517!! OUTPUT 518!! xcname_short= short version of the libXC name 519!! 520!!================================================================= 521 subroutine libxc_getshortname(xcname,xcname_short) 522 523 implicit none 524 character*(*),intent(in) :: xcname 525 character*(*),intent(out) :: xcname_short 526 527#if defined HAVE_LIBXC 528!---- Local variables 529 integer :: i_plus 530 character*50 :: xcstrg(2) 531 532!------------------------------------------------------------------ 533!---- Executable code 534 535 i_plus=index(xcname,'+') 536 if (i_plus<=0) then 537 xcstrg(1)=trim(xcname) 538 xcstrg(2)="" 539 else 540 xcstrg(1)=trim(xcname(1:i_plus-1)) 541 xcstrg(2)=trim(xcname(i_plus+1:)) 542 end if 543 call uppercase(xcstrg(1)) 544 call uppercase(xcstrg(2)) 545 546 xcname_short="" 547 if (xcstrg(1)(1:3)=="XC_") xcname_short=trim(xcname_short) //trim(xcstrg(1)(4:)) 548 if (xcstrg(1)(1:6)=="LIBXC_") xcname_short=trim(xcname_short) //trim(xcstrg(1)) 549 if (xcstrg(2)(1:3)=="XC_") xcname_short=trim(xcname_short)//"+"//trim(xcstrg(2)(4:)) 550 if (xcstrg(2)(1:6)=="LIBXC_") xcname_short=trim(xcname_short)//"+"//trim(xcstrg(2)) 551 if (trim(xcname_short)=="") xcname_short=xcname 552 553#else 554 xcname_short=xcname 555#endif 556 557 end subroutine libxc_getshortname 558 559 560!!================================================================= 561!! NAME 562!! libxc_init_func 563!! 564!! FUNCTION 565!! Initialize libXC functional(s) 566!! 567!! INPUTS 568!! id(2)= libXC id(s) a XC functional 569!! nsp=number of psin component 570!! 571!!================================================================= 572 subroutine libxc_init_func(id,nsp) 573 574 implicit none 575 integer,intent(in) :: id(2) 576 integer,intent(in) :: nsp 577 578!---- Local variables 579 integer :: ii 580 type(libxc_functional_t),pointer :: xc_func 581#if defined HAVE_LIBXC && defined HAVE_FC_ISO_C_BINDING 582 integer :: flags 583 integer(C_INT) :: func_id_c,npar_c,nspin_c,success_c 584 real(C_DOUBLE) :: alpha_c,beta_c,omega_c 585 real(C_DOUBLE) :: param_c(3) 586 character(kind=C_CHAR,len=1),pointer :: strg_c 587 type(C_PTR) :: func_ptr_c 588#endif 589 590!------------------------------------------------------------------ 591!---- Executable code 592 593!Check libXC 594 if (.not.libxc_check(stop_if_error=.true.)) return 595 if (.not.libxc_constants_initialized) call libxc_constants_load() 596 597 libxc_funcs(1:2)%id=id(1:2) 598 599 do ii=1,2 600 601 xc_func => libxc_funcs(ii) 602 603 xc_func%family=XC_FAMILY_UNKNOWN 604 xc_func%xckind=-1 605 xc_func%nspin=nsp 606 xc_func%has_exc=.false. 607 xc_func%has_vxc=.false. 608 xc_func%has_fxc=.false. 609 xc_func%has_kxc=.false. 610 xc_func%needs_laplacian=.false. 611 xc_func%is_hybrid=.false. 612 xc_func%hyb_mixing=0.d0 613 xc_func%hyb_mixing_sr=0.d0 614 xc_func%hyb_range=0.d0 615 616 if (xc_func%id<=0) cycle 617 618! Get XC functional family 619 libxc_funcs%family=libxc_family_from_id(xc_func%id) 620! Live dangerously 621! if (xc_func%family/=XC_FAMILY_LDA .and. & 622!& xc_func%family/=XC_FAMILY_GGA .and. & 623!& xc_func%family/=XC_FAMILY_MGGA.and. & 624!& xc_func%family/=XC_FAMILY_HYB_GGA) then 625! write(std_out,'(a,i4,a)') 'The LibXC functional family ',xc_func%family, & 626!& ' is currently unsupported by ATOMPAW!' 627! write(std_out,'(a)') '(-1 means the family is unknown to the LibXC itself)' 628! write(std_out,'(a)') 'Please consult the LibXC documentation.' 629! stop 630! end if 631 632#if defined HAVE_LIBXC && defined HAVE_FC_ISO_C_BINDING 633 634! Allocate functional 635 func_ptr_c=xc_func_type_malloc() 636 call c_f_pointer(func_ptr_c,xc_func%conf) 637 638! Initialize functional 639 func_id_c=int(xc_func%id,kind=C_INT) 640 nspin_c=int(nsp,kind=C_INT) 641 success_c=xc_func_init(xc_func%conf,func_id_c,nspin_c) 642 if (success_c/=0) then 643 write(std_out,'(a)') 'Error in libXC functional initialization!' 644 stop 645 end if 646 647! Special treatment for LDA_C_XALPHA functional 648 if (xc_func%id==libxc_getid_fromName('XC_LDA_C_XALPHA')) then 649 param_c(1)=real(0.d0,kind=C_DOUBLE);npar_c=int(1,kind=C_INT) 650 call xc_func_set_params(xc_func%conf,param_c,npar_c) 651 end if 652 653! Get functional kind 654 xc_func%xckind=int(xc_get_info_kind(xc_func%conf)) 655 656! Get functional flags 657 flags=int(xc_get_info_flags(xc_func%conf)) 658 xc_func%has_exc=(iand(flags,XC_FLAGS_HAVE_EXC)>0) 659 xc_func%has_vxc=(iand(flags,XC_FLAGS_HAVE_VXC)>0) 660 xc_func%has_fxc=(iand(flags,XC_FLAGS_HAVE_FXC)>0) 661 xc_func%has_kxc=(iand(flags,XC_FLAGS_HAVE_KXC)>0) 662 xc_func%needs_laplacian=(iand(flags,XC_FLAGS_NEEDS_LAPLACIAN)>0) 663 664! Retrieve parameters for hybrid functionals 665 xc_func%is_hybrid=(xc_func_is_hybrid_from_id(xc_func%id)==1) 666 if (xc_func%is_hybrid) then 667 call xc_hyb_cam_coef(xc_func%conf,omega_c,alpha_c,beta_c) 668 xc_func%hyb_mixing=real(alpha_c,kind=8) 669 xc_func%hyb_mixing_sr=real(beta_c,kind=8) 670 xc_func%hyb_range=real(omega_c,kind=8) 671 endif 672 673#endif 674 675 end do ! loop on functionals 676 677!Print functional(s) information 678 call libxc_print_func(6) 679 680 end subroutine libxc_init_func 681 682 683!!================================================================= 684!! NAME 685!! libxc_end_func 686!! 687!! FUNCTION 688!! Destroy libXC functional(s) 689!! 690!!================================================================= 691 subroutine libxc_end_func() 692 693 implicit none 694 695#if defined HAVE_LIBXC 696 697!---- Local variables 698 integer :: ii 699 type(libxc_functional_t),pointer :: xc_func 700 701!------------------------------------------------------------------ 702!---- Executable code 703 704 do ii = 1,2 705 706 xc_func => libxc_funcs(ii) 707 708 if (xc_func%id <= 0) cycle 709 xc_func%id=-1 710 xc_func%family=-1 711 xc_func%xckind=-1 712 xc_func%nspin=1 713 xc_func%abi_ixc=huge(0) 714 xc_func%has_exc=.false. 715 xc_func%has_vxc=.false. 716 xc_func%has_fxc=.false. 717 xc_func%has_kxc=.false. 718 xc_func%needs_laplacian=.false. 719 xc_func%is_hybrid=.false. 720 xc_func%hyb_mixing=0.d0 721 xc_func%hyb_mixing_sr=0.d0 722 xc_func%hyb_range=0.d0 723#if defined HAVE_FC_ISO_C_BINDING 724 if (associated(xc_func%conf)) then 725 call xc_func_end(xc_func%conf) 726 call xc_func_type_free(c_loc(xc_func%conf)) 727 end if 728#endif 729 730 end do 731 732#endif 733 734 end subroutine libxc_end_func 735 736 737!!================================================================= 738!! NAME 739!! libxc_print_func 740!! 741!! FUNCTION 742!! Print libXC functionnal(s) details 743!! 744!! INPUTS 745!! unt= logical unit to print 746!! 747!!================================================================= 748 subroutine libxc_print_func(unt) 749 750 implicit none 751 integer :: unt 752 753#if defined HAVE_LIBXC 754 755!---- Local variables 756 integer :: ii 757 character(len=500) :: msg 758 type(libxc_functional_t),pointer :: xc_func 759#if defined HAVE_FC_ISO_C_BINDING 760 integer(C_INT) :: iref_c 761 character(kind=C_CHAR,len=1),pointer :: strg_c 762#endif 763 764!------------------------------------------------------------------ 765!---- Executable code 766 767 do ii=1,2 768 769 xc_func => libxc_funcs(ii) 770 if (xc_func%id<=0) cycle 771 772 if (xc_func%xckind==XC_EXCHANGE) then 773 write(unt,'(a)') 'Exchange functional (LibXC):' 774 else if (xc_func%xckind==XC_CORRELATION) then 775 write(unt,'(a)') 'Correlation functional (LibXC):' 776 else if (xc_func%xckind==XC_EXCHANGE_CORRELATION) then 777 write(unt,'(a)') 'Exchange-Correlation functional (LibXC):' 778 end if 779 780#if defined HAVE_FC_ISO_C_BINDING 781 call c_f_pointer(xc_get_info_name(xc_func%conf),strg_c) 782 call xc_char_to_f(strg_c,msg) 783 784 iref_c=0 785 do while (iref_c>=0) 786 call c_f_pointer(xc_get_info_refs(xc_func%conf,iref_c),strg_c) 787 if (associated(strg_c)) then 788 call xc_char_to_f(strg_c,msg) 789 write(unt,'(2x,a)') trim(msg) 790 iref_c=iref_c+1 791 else 792 iref_c=-1 793 end if 794 end do 795#endif 796 797 end do 798 799#endif 800 801 end subroutine libxc_print_func 802 803 804!!================================================================= 805!! libxc_family_from_id 806!! 807!! FUNCTION 808!! Return family of a XC functional from its id 809!! 810!! INPUTS 811!! xcid= id of a LibXC functional 812!! 813!!================================================================= 814 function libxc_family_from_id(xcid) 815 816 implicit none 817 integer :: libxc_family_from_id 818 integer,intent(in) :: xcid 819 820#if defined HAVE_LIBXC 821 822!---- Local variables 823#if defined HAVE_FC_ISO_C_BINDING 824 integer(C_INT) :: xcid_c 825#endif 826 827!------------------------------------------------------------------ 828!---- Executable code 829 830#if defined HAVE_FC_ISO_C_BINDING 831 xcid_c=int(xcid,kind=C_INT) 832 libxc_family_from_id=int(xc_family_from_id(xcid_c,C_NULL_PTR,C_NULL_PTR)) 833#endif 834 835#else 836 libxc_family_from_id=-1 837#endif 838 839end function libxc_family_from_id 840 841 842!!================================================================= 843!! NAME 844!! libxc_islda 845!! 846!! FUNCTION 847!! Return TRUE is LibXC functional is LDA 848!! 849!!================================================================= 850 function libxc_islda() 851 852 implicit none 853 logical :: libxc_islda 854 855!------------------------------------------------------------------ 856!---- Executable code 857 858 libxc_islda = .false. 859 if (.not.libxc_constants_initialized) call libxc_constants_load() 860 861 libxc_islda = (any(libxc_funcs(:)%family == XC_FAMILY_LDA)) 862 863 end function libxc_islda 864 865 866!!================================================================= 867!! NAME 868!! libxc_isgga 869!! 870!! FUNCTION 871!! Return TRUE is LibXC functional is GGA 872!! 873!!================================================================= 874 function libxc_isgga() 875 876 implicit none 877 logical :: libxc_isgga 878 879!------------------------------------------------------------------ 880!---- Executable code 881 882 libxc_isgga = .false. 883 if (.not.libxc_constants_initialized) call libxc_constants_load() 884 885 libxc_isgga = (any(libxc_funcs(:)%family == XC_FAMILY_GGA)) 886 887 end function libxc_isgga 888 889 890!!================================================================= 891!! libxc_ismgga 892!! 893!! FUNCTION 894!! Return TRUE is LibXC functional is meta-GGA 895!! 896!!================================================================= 897 function libxc_ismgga() 898 899 implicit none 900 logical :: libxc_ismgga 901 902!------------------------------------------------------------------ 903!---- Executable code 904 905 libxc_ismgga = .false. 906 if (.not.libxc_constants_initialized) call libxc_constants_load() 907 908 libxc_ismgga = (any(libxc_funcs(:)%family == XC_FAMILY_MGGA)) 909 910 end function libxc_ismgga 911 912 913!!================================================================= 914 function libxc_needs_laplacian() 915 916 implicit none 917 logical :: libxc_needs_laplacian 918 919!------------------------------------------------------------------ 920!---- Executable code 921 922 libxc_needs_laplacian = .false. 923 924 libxc_needs_laplacian = (any(libxc_funcs(:)%needs_laplacian)) 925 926 end function libxc_needs_laplacian 927 928 929!!================================================================= 930!! NAME 931!! libxc_ishybrid 932!! 933!! FUNCTION 934!! Return TRUE is LibXC functional is Hybrid 935!! 936!!================================================================= 937 function libxc_ishybrid() 938 939 implicit none 940 logical :: libxc_ishybrid 941 942!------------------------------------------------------------------ 943!---- Executable code 944 945 libxc_ishybrid = .false. 946 947 libxc_ishybrid=(any(libxc_funcs(:)%is_hybrid)) 948 949 end function libxc_ishybrid 950 951 952!!================================================================= 953!! NAME 954!! libxc_is_hybrid_from_id 955!! 956!! FUNCTION 957!! Test function to identify whether a functional is hybrid or not, from its id 958!! 959!! INPUTS 960!! xcid= id of a LibXC functional 961!! 962!!================================================================= 963 function libxc_is_hybrid_from_id(xcid) 964 965 implicit none 966 logical :: libxc_is_hybrid_from_id 967 integer,intent(in) :: xcid 968 969!---- Local variables 970#if defined HAVE_FC_ISO_C_BINDING 971 integer(C_INT) :: xcid_c 972#endif 973 974!------------------------------------------------------------------ 975!---- Executable code 976 977#if defined HAVE_FC_ISO_C_BINDING 978 xcid_c=int(xcid,kind=C_INT) 979 libxc_is_hybrid_from_id =(xc_func_is_hybrid_from_id(xcid_c)==1) 980#else 981 libxc_is_hybrid_from_id = .false. if (.false.) write(std_out,*) xcid 982#endif 983 984end function libxc_is_hybrid_from_id 985 986 987!!================================================================= 988!! NAME 989!! libxc_functionals_has_kxc 990!! 991!! FUNCTION 992!! Test function to identify whether the presently used functional 993!! provides Kxc or not (fxc in the libXC convention) 994!! 995!!================================================================= 996function libxc_has_kxc() 997 998 implicit none 999 logical :: libxc_has_kxc 1000 1001!---- Local variables 1002 integer :: ii 1003 1004!------------------------------------------------------------------ 1005!---- Executable code 1006 1007 libxc_has_kxc=.true. 1008 1009 do ii=1,2 1010 if (.not.libxc_funcs(ii)%has_fxc) libxc_has_kxc=.false. 1011 end do 1012 1013end function libxc_has_kxc 1014 1015 1016!!================================================================= 1017!! NAME 1018!! libxc_functionals_has_k3xc 1019!! 1020!! FUNCTION 1021!! Test function to identify whether the presently used functional 1022!! provides K3xc or not (kxc in the libXC convention) 1023!! 1024!!================================================================= 1025function libxc_has_k3xc() 1026 1027 implicit none 1028 logical :: libxc_has_k3xc 1029 1030!---- Local variables 1031 integer :: ii 1032 1033!------------------------------------------------------------------ 1034!---- Executable code 1035 1036 libxc_has_k3xc=.true. 1037 1038 do ii=1,2 1039 if (.not.libxc_funcs(ii)%has_kxc) libxc_has_k3xc=.false. 1040 end do 1041 1042end function libxc_has_k3xc 1043 1044 1045!!================================================================= 1046!! NAME 1047!! libxc_functionals_nspin 1048!! 1049!! FUNCTION 1050!! Returns the number of spin components for the XC functionals 1051!! 1052!!================================================================= 1053function libxc_nspin() 1054 1055 implicit none 1056 integer :: libxc_nspin 1057 1058!------------------------------------------------------------------ 1059!---- Executable code 1060 1061 libxc_nspin = 1 1062 if (any(libxc_funcs(:)%nspin==2)) libxc_nspin=2 1063 1064end function libxc_nspin 1065 1066 1067!!================================================================= 1068!! NAME 1069!! libxc_getvxc 1070!! 1071!! FUNCTION 1072!! Return XC potential and energy, from input density (gradient etc...) 1073!! 1074!! INPUTS 1075!! npts= number of of points for the density 1076!! nsp= number of spin-density components 1077!! rho(npts,nsp)= electronic density 1078!! [grho(npts,2*nsp-1)]= gradient of the density (optional) 1079!! [lrho(npts,2*nsp-1)]= laplacian of the density (optional) 1080!! [tau(npts,nsp)]= sum of squared gradient of occ wf's (optional) 1081!! 1082!! OUTPUT 1083!! exc(npts)=XC energy density 1084!! vxc(npts,nsp)=derivative of the energy density wrt to the density 1085!! (df_xc/dn) 1086!! [vxclrho(npts,nsp)]=derivative of the energy density wrt to the density laplacian 1087!! (=df_xc/dLapl(n)) 1088!! [vxctau(npts,nsp)]=derivative of the energy density wrt to the kinetic energy density 1089!! (=df_xc/dtau) 1090!! [vxcgr(npts,2*nsp-1)]=derivative of the energy density wrt to the gradient of the density 1091!! (=1/g * dfxc/dg = 2*df_xc/dsigma) 1092!! 1093!! NOTES 1094!! nsp=1 : rho is total density (not half) 1095!! grho is abs(grad(rho)) 1096!! nsp=2 : rho is [rho^up,rho^dn] 1097!! grho is [abs(grad(rho^up)),abs(grad(rho^dn)),abs(grad(rho^tot))] 1098!! 1099!! All energies (input/output) should be given in Rydberg units 1100 1101!! 1102!!================================================================= 1103 subroutine libxc_getvxc(npts,exc,vxc,nsp,rho, & 1104& grho,lrho,tau,vxcgr,vxclrho,vxctau) ! Optional arguments 1105 1106 implicit none 1107 integer, intent(in) :: npts,nsp 1108 real(8),intent(inout) :: exc(npts),vxc(npts,nsp) 1109 real(8),intent(in) :: rho(npts,nsp) 1110 real(8),intent(in),optional :: grho(npts,2*nsp-1) 1111 real(8),intent(in),optional :: lrho(npts,nsp) 1112 real(8),intent(in),optional :: tau(npts,nsp) 1113 real(8),intent(inout),optional :: vxcgr(npts,2*nsp-1) 1114 real(8),intent(inout), optional :: vxclrho(npts,nsp) 1115 real(8),intent(inout), optional :: vxctau(npts,nsp) 1116 1117#if defined HAVE_LIBXC 1118 1119!---- Local variables 1120 real(8),parameter :: tol=1.d-14 1121 1122 integer :: ii,ipts,izero 1123 logical :: is_lda,is_gga,is_mgga,needs_laplacian 1124 real(8),target :: exctmp 1125 real(8),target :: rhotmp(nsp),vxctmp(nsp),sigma(3),vsigma(3) 1126 real(8),target :: v2rho2(3),v2rhosigma(6),v2sigma2(6) 1127 real(8),target :: v3rho3(4),v3rho2sigma(9),v3rhosigma2(12),v3sigma3(10) 1128 real(8),target :: lrhotmp(nsp),tautmp(nsp),vlrho(nsp),vtau(nsp) 1129 1130#if defined HAVE_LIBXC && defined HAVE_FC_ISO_C_BINDING 1131 type(C_PTR) :: rho_c,sigma_c,lrho_c,tau_c 1132 type(C_PTR) :: exc_c(2),vxc_c(2),vsigma_c(2),vlrho_c(2),vtau_c(2) 1133 type(C_PTR) :: v2rho2_c(2),v2rhosigma_c(2),v2sigma2_c(2) 1134 type(C_PTR) :: v3rho3_c(2),v3rho2sigma_c(2),v3rhosigma2_c(2),v3sigma3_c(2) 1135#endif 1136 1137!------------------------------------------------------------------ 1138!---- Executable code 1139 1140 if (.not.libxc_constants_initialized) call libxc_constants_load() 1141 1142 is_lda=libxc_islda() 1143 is_gga=libxc_isgga() 1144 is_mgga=libxc_ismgga() 1145 needs_laplacian=libxc_needs_laplacian() 1146 1147!Why? 1148 if (nsp==2.and.is_mgga) stop 'spin not available for mGGA!' 1149 1150 if (is_gga.and.((.not.present(vxcgr).or.(.not.present(grho))))) then 1151 write(std_out,'(/,2x,a)') "Bug in libxc_getvxc:" 1152 write(std_out,'(2x,a)') " GGA called without grho or vxcgr!" 1153 stop 1154 end if 1155 if(needs_laplacian.and.(.not.present(lrho))) then 1156 write(std_out,'(/,2x,a)') "Error in libxc_getvxc:" 1157 write(std_out,'(2x,a)') " getvxc need Laplacian of density!" 1158 stop 1159 endif 1160!Need to add more consistency tests 1161 1162!Initialize all output arrays to zero 1163 exc=0.d0 ; vxc=0.d0 1164 if (is_gga.or.is_mgga.and.present(vxcgr)) vxcgr=0.d0 1165 if (is_mgga.and.present(vxclrho)) vxclrho=0.d0 1166 if (is_mgga.and.present(vxctau)) vxctau=0.d0 1167 1168 1169!Filter density/gradient when density goes to zero 1170!This is useless ; libxc has its own filtering process 1171 izero=npts 1172!do ipts=npts,2,-1 1173! if (all(rho(ipts,:)<tol)) izero=ipts-1 1174!end do 1175 1176!Adjust zero-density threshold 1177 do ii = 1,2 1178 if (libxc_funcs(ii)%id>0) then 1179 call xc_func_set_density_threshold(libxc_funcs(ii)%conf,tol) 1180 end if 1181 end do 1182 1183!Define C pointers to libXC routine arguments 1184#if defined HAVE_FC_ISO_C_BINDING 1185 do ii = 1,2 1186 exc_c(ii)=C_NULL_PTR 1187 vxc_c(ii)=C_NULL_PTR 1188 vsigma_c(ii)=C_NULL_PTR 1189 vlrho_c(ii)=C_NULL_PTR 1190 vtau_c(ii)=C_NULL_PTR 1191 v2rho2_c(ii)=C_NULL_PTR 1192 v2sigma2_c(ii)=C_NULL_PTR 1193 v2rhosigma_c(ii)=C_NULL_PTR 1194 v3rho3_c(ii)=C_NULL_PTR 1195 v3sigma3_c(ii)=C_NULL_PTR 1196 v3rho2sigma_c(ii)=C_NULL_PTR 1197 v3rhosigma2_c(ii)=C_NULL_PTR 1198 if (libxc_funcs(ii)%has_exc) then 1199 exc_c(ii)=c_loc(exctmp) 1200 end if 1201 if (libxc_funcs(ii)%has_vxc) then 1202 vxc_c(ii)=c_loc(vxctmp) 1203 vsigma_c(ii)=c_loc(vsigma) 1204 vlrho_c(ii)=c_loc(vlrho) 1205 vtau_c(ii)=c_loc(vtau) 1206 endif 1207 enddo 1208#endif 1209 1210!Initialize temporary arrays 1211#if defined HAVE_LIBXC && defined HAVE_FC_ISO_C_BINDING 1212 rhotmp=0.d0 ; rho_c=c_loc(rhotmp) 1213 if (is_gga.or.is_mgga) then 1214 sigma=0.d0 ; sigma_c=c_loc(sigma) 1215 endif 1216 if (is_mgga) then 1217 tautmp=0.d0; tau_c=c_loc(tautmp) 1218 lrhotmp=0.d0; lrho_c=c_loc(lrhotmp) 1219 end if 1220#endif 1221 1222!Loop over points 1223 do ipts=1,npts 1224 1225 !Load density (and gradient) for this point 1226 vxctmp=0.d0;exctmp=0.d0 1227 if (ipts<=izero) then 1228 rhotmp(1:nsp)=rho(ipts,1:nsp) 1229 else 1230 rhotmp=tol 1231 end if 1232 if (is_gga.or.is_mgga) then 1233 if (ipts<=izero) then 1234 if (nsp==1) then 1235 !AtomPAW passes |grho| while LibXC needs |grho|^2 1236 sigma(1)=grho(ipts,1)**2 1237 else 1238 !AtomPAW passes |grho_up|, |grho_dn|, and |grho_tot| 1239 !while Libxc needs |grho_up|^2, grho_up.grho_dn, and |grho_dn|^2 1240 sigma(1)= grho(ipts,1)**2 1241 sigma(3)= grho(ipts,2)**2 1242 sigma(2)=(grho(ipts,3)**2-sigma(1)-sigma(3))*0.5d0 1243 end if 1244 else 1245 sigma=0.d0 1246 end if 1247 end if 1248 if (is_mgga) then 1249 if (ipts<=izero) then 1250 !AtomPAW passes tau (Ry) while LibXC needs tau (Ha) 1251 if (nsp==1) then 1252 tautmp(1)=tau(ipts,1)/2.d0 ! From Ry to Ha units 1253 lrhotmp(1)=0.d0 1254 if (needs_laplacian) lrhotmp(1)=lrho(ipts,1) 1255 else 1256 tautmp(1)= tau(ipts,1)/2.d0 ! From Ry to Ha units 1257 tautmp(2)= tau(ipts,2)/2.d0 ! From Ry to Ha units 1258 if (needs_laplacian) then 1259 lrhotmp(1)= lrho(ipts,1) 1260 lrhotmp(2)= lrho(ipts,2) 1261 endif 1262 end if 1263 else 1264 tautmp=0.d0 1265 if (needs_laplacian) lrhotmp=0.d0 1266 end if 1267 end if 1268 1269! Loop over functionals 1270 do ii=1,2 1271 if (libxc_funcs(ii)%id<=0) cycle 1272 1273! Get the energy and the potential (and possibly the other derivatives) 1274#if defined HAVE_FC_ISO_C_BINDING 1275! ===== LDA ===== 1276 if (libxc_funcs(ii)%family==XC_FAMILY_LDA) then 1277 exctmp=0.d0 ; vxctmp=0.d0 1278 call xc_get_lda(libxc_funcs(ii)%conf,1,rho_c, & 1279& exc_c(ii),vxc_c(ii),v2rho2_c(ii),v3rho3_c(ii)) 1280! ===== GGA ===== 1281 else if (libxc_funcs(ii)%family==XC_FAMILY_GGA.or. & 1282& libxc_funcs(ii)%family==XC_FAMILY_HYB_GGA) then 1283 exctmp=0.d0 ; vxctmp=0.d0 ; vsigma=0.d0 1284 call xc_get_gga(libxc_funcs(ii)%conf,1,rho_c,sigma_c, & 1285& exc_c(ii),vxc_c(ii),vsigma_c(ii), & 1286& v2rho2_c(ii),v2rhosigma_c(ii),v2sigma2_c(ii), & 1287& v3rho3_c(ii),v3rho2sigma_c(ii),v3rhosigma2_c(ii),v3sigma3_c(ii)) 1288! ===== mGGA ===== 1289 else if (libxc_funcs(ii)%family==XC_FAMILY_MGGA) then 1290 exctmp=0.d0 ; vxctmp=0.d0 ; vsigma=0.d0 ; vlrho=0.d0 ; vtau=0.d0 1291 call xc_get_mgga(libxc_funcs(ii)%conf,1,rho_c,sigma_c,lrho_c,tau_c, & 1292& exc_c(ii),vxc_c(ii),vsigma_c(ii),vlrho_c(ii),vtau_c(ii), & 1293& C_NULL_PTR,C_NULL_PTR,C_NULL_PTR,C_NULL_PTR,C_NULL_PTR, & 1294& C_NULL_PTR,C_NULL_PTR,C_NULL_PTR,C_NULL_PTR,C_NULL_PTR) 1295 end if 1296#endif 1297 1298 exc(ipts)=exc(ipts)+2.d0*exctmp ! From Ha to Ry 1299 vxc(ipts,1:nsp)=vxc(ipts,1:nsp)+2.d0*vxctmp(1:nsp) ! From Ha to Ry 1300 1301! Additional output in case of GGA or meta-GGA 1302 if (is_gga.or.is_mgga.and.present(vxcgr)) then 1303 if (nsp==1) then 1304 !Note: for nsp=1, factor 2 comes from sigma=grad**2 1305 vxcgr(ipts,1)=vxcgr(ipts,1)+2.d0*vsigma(1) *2.d0 ! From Ha to Ry 1306 else 1307 ! vxcgrho(npts,3)= 1/|grad_rho_up| dfx/d(|grad_rho_up|) 1308 ! 1/|grad_rho_dn| dfx/d(|grad_rho_dn|) 1309 ! and 1/|grad_rho| dfc/d(|grad_rho|) 1310 ! These formulas have been checked! 1311 if (libxc_funcs(ii)%xckind==XC_EXCHANGE) then 1312 vxcgr(ipts,1) = vxcgr(ipts,1) + 2.d0*vsigma(1) *2.d0 ! From Ha to Ry 1313 vxcgr(ipts,2) = vxcgr(ipts,2) + 2.d0*vsigma(3) *2.d0 ! From Ha to Ry 1314 else 1315 vxcgr(ipts,1) = vxcgr(ipts,1) + (2.d0*vsigma(1) - vsigma(2)) *2.d0 ! From Ha to Ry 1316 vxcgr(ipts,2) = vxcgr(ipts,2) + (2.d0*vsigma(3) - vsigma(2)) *2.d0 ! From Ha to Ry 1317 vxcgr(ipts,3) = vxcgr(ipts,3) + vsigma(2) *2.d0 ! From Ha to Ry 1318 end if 1319 end if 1320 end if 1321 1322! Additional output in case of meta-GGA 1323 if (is_mgga.and.present(vxctau)) then 1324 vxctau(ipts,1:nsp)=vxctau(ipts,1:nsp)+2.d0*vtau(1:nsp) ! From Ha to Ry 1325 end if 1326 if (is_mgga.and.present(vxclrho)) then 1327 vxclrho(ipts,1:nsp)=vxclrho(ipts,1:nsp)+2.d0*vlrho(1:nsp) ! From Ha to Ry 1328 end if 1329 1330 end do ! loop over functional(s) 1331 end do ! loop over points 1332 1333#endif 1334 end subroutine libxc_getvxc 1335 1336 1337!!================================================================= 1338!! NAME 1339!! libxc_constants_load 1340!! 1341!! FUNCTION 1342!! Load libXC constants from C headers 1343!! 1344!!================================================================= 1345 subroutine libxc_constants_load() 1346 1347 implicit none 1348 1349#if defined HAVE_LIBXC 1350 1351!---- Local variables 1352#if defined HAVE_FC_ISO_C_BINDING 1353 integer(C_INT) :: i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13 1354#endif 1355 1356!------------------------------------------------------------------ 1357!---- Executable code 1358 1359#if defined HAVE_FC_ISO_C_BINDING 1360 call xc_get_singleprecision_constant(i1) 1361 XC_SINGLE_PRECISION = int(i1) 1362 call xc_get_family_constants(i1,i2,i3,i4,i5,i6,i7,i8) 1363 XC_FAMILY_UNKNOWN = int(i1) 1364 XC_FAMILY_LDA = int(i2) 1365 XC_FAMILY_GGA = int(i3) 1366 XC_FAMILY_MGGA = int(i4) 1367 XC_FAMILY_LCA = int(i5) 1368 XC_FAMILY_OEP = int(i6) 1369 XC_FAMILY_HYB_GGA = int(i7) 1370 XC_FAMILY_HYB_MGGA = int(i8) 1371 call xc_get_flags_constants(i1,i2,i3,i4,i5,i6) 1372 XC_FLAGS_HAVE_EXC = int(i1) 1373 XC_FLAGS_HAVE_VXC = int(i2) 1374 XC_FLAGS_HAVE_FXC = int(i3) 1375 XC_FLAGS_HAVE_KXC = int(i4) 1376 XC_FLAGS_HAVE_LXC = int(i5) 1377 XC_FLAGS_NEEDS_LAPLACIAN = int(i6) 1378 call xc_get_kind_constants(i1,i2,i3,i4) 1379 XC_EXCHANGE = int(i1) 1380 XC_CORRELATION = int(i2) 1381 XC_EXCHANGE_CORRELATION = int(i3) 1382 XC_KINETIC = int(i4) 1383 call xc_get_hybrid_constants(i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12,i13) 1384 XC_HYB_NONE = int(i1) 1385 XC_HYB_FOCK = int(i2) 1386 XC_HYB_PT2 = int(i3) 1387 XC_HYB_ERF_SR = int(i4) 1388 XC_HYB_YUKAWA_SR = int(i5) 1389 XC_HYB_GAUSSIAN_SR = int(i6) 1390 XC_HYB_SEMILOCAL = int(i7) 1391 XC_HYB_HYBRID = int(i8) 1392 XC_HYB_CAM = int(i9) 1393 XC_HYB_CAMY = int(i10) 1394 XC_HYB_CAMG = int(i11) 1395 XC_HYB_DOUBLE_HYBRID = int(i12) 1396 XC_HYB_MIXTURE = int(i13) 1397 libxc_constants_initialized=.true. 1398#endif 1399 1400#else 1401 libxc_constants_initialized=.false. 1402#endif 1403 1404 end subroutine libxc_constants_load 1405 1406 1407!!================================================================= 1408!! NAME 1409!! libxc_check 1410!! 1411!! FUNCTION 1412!! Check if the code has been compiled with libXC 1413!! 1414!! INPUTS 1415!! [stop_if_error]=optional flag; if TRUE the code stops 1416!! if libXC is not correctly used 1417!! 1418!!================================================================= 1419 function libxc_check(stop_if_error) 1420 1421 implicit none 1422 logical :: libxc_check 1423 logical,intent(in),optional :: stop_if_error 1424 1425!---- Local variables 1426 character(len=100) :: msg 1427 1428!------------------------------------------------------------------ 1429!---- Executable code 1430 1431 libxc_check=.true. ; msg="" 1432 1433#if defined HAVE_LIBXC 1434#if defined FC_G95 1435 libxc_check=.false. 1436 msg='LibXC cannot be used with G95 Fortran compiler!' 1437#endif 1438#if defined HAVE_FC_ISO_C_BINDING 1439 if (.not.libxc_constants_initialized) call libxc_constants_load() 1440 if (XC_SINGLE_PRECISION==1) then 1441 libxc_check=.false. 1442 msg='LibXC should be compiled with double precision!' 1443 end if 1444#else 1445 libxc_check=.false. 1446 msg='LibXC cannot be used without ISO_C_BINDING support by the Fortran compiler!' 1447#endif 1448#else 1449 libxc_check=.false. 1450 msg='ATOMPAW was not compiled with LibXC support.' 1451#endif 1452 1453 if (present(stop_if_error)) then 1454 if (stop_if_error.and.trim(msg)/="") then 1455 write(std_out,'(a)') trim(msg) ; stop 1456 end if 1457 end if 1458 1459 end function libxc_check 1460 1461 1462!!================================================================= 1463!! NAME 1464!! xc_char_to_c 1465!! 1466!! FUNCTION 1467!! Helper function to convert a Fortran string to a C string 1468!! Based on a routine by Joseph M. Krahn 1469!! 1470!! INPUTS 1471!! f_string=Fortran string 1472!! 1473!! OUTPUT 1474!! c_string=C string 1475!! 1476!!================================================================= 1477#if defined HAVE_FC_ISO_C_BINDING 1478function xc_char_to_c(f_string) result(c_string) 1479 1480 character(len=*),intent(in) :: f_string 1481 character(kind=C_CHAR,len=1) :: c_string(len_trim(f_string)+1) 1482 1483!---- Local variables 1484 integer :: ii,strlen 1485 1486!------------------------------------------------------------------ 1487!---- Executable code 1488 1489 strlen=len_trim(f_string) 1490 forall(ii=1:strlen) 1491 c_string(ii)=f_string(ii:ii) 1492 end forall 1493 c_string(strlen+1)=C_NULL_CHAR 1494 1495 end function xc_char_to_c 1496#endif 1497 1498 1499!!================================================================= 1500!! NAME 1501!! xc_char_to_f 1502!! 1503!! FUNCTION 1504!! Helper function to convert a C string to a Fortran string 1505!! Based on a routine by Joseph M. Krahn 1506!! 1507!! INPUTS 1508!! c_string=C string 1509!! 1510!! OUTPUT 1511!! f_string=Fortran string 1512!! 1513!!================================================================= 1514#if defined HAVE_FC_ISO_C_BINDING 1515subroutine xc_char_to_f(c_string,f_string) 1516 1517 character(kind=C_CHAR,len=1),intent(in) :: c_string(*) 1518 character(len=*),intent(out) :: f_string 1519 1520!---- Local variables 1521 integer :: ii 1522 1523!------------------------------------------------------------------ 1524!---- Executable code 1525 1526 ii=1 1527 do while(c_string(ii)/=C_NULL_CHAR.and.ii<=len(f_string)) 1528 if (iachar(c_string(ii)) <= 127) then 1529 f_string(ii:ii)=c_string(ii) 1530 else 1531 f_string(ii:ii)="?" 1532 end if 1533 ii=ii+1 1534 end do 1535 if (ii<len(f_string)) f_string(ii:)=' ' 1536 1537 end subroutine xc_char_to_f 1538#endif 1539 1540 1541!!================================================================= 1542 1543end Module libxc_mod 1544