TABLE OF CONTENTS
- ABINIT/m_drivexc
- m_drivexc/check_kxc
- m_drivexc/drivexc
- m_drivexc/echo_xc_name
- m_drivexc/has_k3xc
- m_drivexc/has_kxc
- m_drivexc/mkdenpos
- m_drivexc/size_dvxc
- m_drivexc/xcmult
ABINIT/m_drivexc [ Modules ]
NAME
m_drivexc
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-2024 ABINIT group (MT, MJV, CE, TD, 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 .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 module m_drivexc 24 25 use defs_basis 26 use m_abicore 27 use m_errors 28 use libxc_functionals 29 use m_numeric_tools, only : invcb 30 use m_xciit, only : xciit 31 use m_xcpbe, only : xcpbe 32 use m_xchcth, only : xchcth 33 use m_xclda, only : xcpzca, xcspol, xctetr, xcwign, xchelu, xcxalp, xclb 34 35 implicit none 36 37 private
m_drivexc/check_kxc [ Functions ]
[ Top ] [ m_drivexc ] [ Functions ]
NAME
check_kxc
FUNCTION
Given a XC functional (defined by ixc), check if Kxc and/or K3xc is avalaible.
INPUTS
ixc = internal code for xc functional optdriver=type of calculation (ground-state, response function, GW, ...) [check_k3xc]= optional ; check also k3xc availability
OUTPUT
SOURCE
319 subroutine check_kxc(ixc,optdriver,check_k3xc) 320 321 !Arguments ------------------------------- 322 integer, intent(in) :: ixc,optdriver 323 logical,intent(in),optional :: check_k3xc 324 325 !Local variables ------------------------- 326 logical :: check_k3xc_,kxc_available,k3xc_available 327 character(len=500) :: msg 328 329 ! ********************************************************************* 330 331 check_k3xc_=.false. ; if (present(check_k3xc)) check_k3xc_=check_k3xc 332 333 kxc_available=has_kxc(ixc) 334 k3xc_available=has_k3xc(ixc) 335 336 if (ixc>=0) then 337 if (.not.kxc_available) then 338 write(msg,'(a,i0,3a)') & 339 & 'The selected XC functional (ixc=',ixc,')',ch10,& 340 & 'does not provide Kxc (dVxc/drho) !' 341 end if 342 if (check_k3xc_.and.(.not.k3xc_available)) then 343 write(msg,'(a,i0,3a)') & 344 & 'The selected XC functional (ixc=',ixc,')',ch10,& 345 & 'does not provide K3xc (d^2Vxc/drho^2) !' 346 end if 347 else ! ixc<0 348 if (.not.kxc_available) then 349 write(msg,'(a,i0,7a)') & 350 & 'The selected XC functional (ixc=',ixc,'):',ch10,& 351 & ' <<',trim(libxc_functionals_fullname()),'>>',ch10,& 352 & 'does not provide Kxc (dVxc/drho) !' 353 end if 354 if (check_k3xc_.and.(.not.k3xc_available)) then 355 write(msg,'(a,i0,7a)') & 356 & 'The selected XC functional (ixc=',ixc,'):',ch10,& 357 & ' <<',trim(libxc_functionals_fullname()),'>>',ch10,& 358 & 'does not provide K3xc (d^2Vxc/d^2rho) !' 359 end if 360 end if 361 362 if (.not.kxc_available) then 363 write(msg,'(7a)') trim(msg),ch10,& 364 & 'However, with the current input options, ABINIT needs Kxc.',ch10,& 365 & '>Possible action:',ch10,& 366 & 'Change the XC functional in psp file or input file.' 367 if (optdriver==0) then 368 write(msg,'(13a)') trim(msg),ch10,& 369 & '>Possible action (2):',ch10,& 370 & 'If you are using density mixing for the SCF cycle',ch10,& 371 & '(iscf>=10, which is the default for PAW),',ch10,& 372 & 'change to potential mixing (iscf=7, for instance).',ch10,& 373 & '>Possible action (3):',ch10,& 374 & 'Switch to another value of densfor_pred (=5, for instance).' 375 end if 376 ABI_ERROR(msg) 377 else if (check_k3xc_.and.(.not.k3xc_available)) then 378 write(msg,'(13a)') trim(msg),ch10,& 379 & 'However, with the current input options, ABINIT needs K3xc.',ch10,& 380 & '>Possible actions:',ch10,& 381 & '- Recompile libXC using --enable-kxc.',ch10,& 382 & ' or',ch10,& 383 & '- Change the XC functional in psp file or input file:',ch10,& 384 & ' use one of the internal LDA (ixc=3, 7 to 15, 23, 24).' 385 ABI_ERROR(msg) 386 end if 387 388 end subroutine check_kxc
m_drivexc/drivexc [ Functions ]
[ Top ] [ m_drivexc ] [ Functions ]
NAME
drivexc
FUNCTION
Driver of XC functionals. Treat spin-polarized as well as non-spin-polarized. Treat local approximations, GGAs, meta-GGAs or hybrid functionals. Optionally, deliver the XC kernel, or even the derivative of the XC kernel (the third derivative of the XC energy)
INPUTS
ixc=index of the XC functional xclevel=XC functional level (lda, gga, etc...) usegradient=[flag] 1 if the XC functional depends on density gradient (grho2_updn) uselaplacian=[flag] 1 if the XC functional depends on density laplacian (lrho_updn) usekden=[flag] 1 if the XC functional depends on kinetic energy density (tau_updn) 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) npts=number of real space points on which the density is provided nspden=number of spin-density components (1 or 2) nvxcgrho=number of components of 1st-derivative of Exc wrt density gradient (nvxcgrho) nvxclrho=number of components of 1st-derivative of Exc wrt density laplacian (nvxclrho) nvxctau=number of components of 1st-derivative of Exc wrt kinetic energy density (nvxctau) ndvxc=number of components of 1st-derivative of Vxc (dvxc) nd2vxc=number of components of 2nd-derivative of Vxc (d2vxc) rho_updn(npts,nspden)=spin-up and spin-down densities In the calling routine, spin-down density must be equal to spin-up density. If nspden=1, only spin-up density must be given (half the total density). If nspden=2, spin-up and spin-down densities must be given. === Optional input arguments === [grho2_updn(npts,(2*nspden-1)*usegradient)]=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. [lrho_updn(npts,nspden*uselaplacian)]=the Laplacian of spin-up and spin-down densities. If nspden=1, only the spin-up Laplacian density must be given and must be equal to the spin-up Laplacian density. If nspden=2, the Laplacian of spin-up and spin-down densities must be given. [tau_updn(npts,nspden*usekden)]=the spin-up and spin-down kinetic energy densities. If nspden=1, only the spin-up kinetic energy density must be given and must be equal to the half the total kinetic energy density. If nspden=2, the spin-up and spin-down kinetic energy densities must be given. [exexch]=choice of <<<local>>> exact exchange. Active if exexch=3 (only for GGA, and NOT for libxc) [el_temp]= electronic temperature (to be used for finite temperature XC functionals) [hyb_mixing]= mixing parameter for the native PBEx functionals (ixc=41 and 42) [xc_funcs(2)]= <type(libxc_functional_type)>: libxc XC functionals.
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 === [vxcgrho(npts,nvxcgrho)]=1st-derivative of the xc energy wrt density gradient> = 1/$|grad \rho_up|$ (d($\rho$*exc)/d($|grad \rho_up|$)) 1/$|grad \rho_dn|$ (d($\rho$*exc)/d($|grad \rho_dn|$)) 1/$|grad \rho|$ (d($\rho$*exc)/d($|grad \rho|$)) [vxclrho(npts,nvxclrho)]=1st-derivative of the xc energy wrt density laplacian. = d($\rho$*exc)/d($\lrho_up$) d($\rho$*exc)/d($\lrho_down$) [vxctau(npts,nvxctau)]=1st-derivative of the xc energy wrt kinetic energy density. = d($\rho$*exc)/d($\tau_up$) d($\rho$*exc)/d($\tau_down$) [dvxc(npts,ndvxc)]=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)= 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))) Note about mGGA: 2nd derivatives involving Tau or Laplacian are not output [d2vxc(npts,nd2vxc)]=second derivative of the XC potential=3rd order derivative of XC energy === Only if abs(order)>1 === === At present only available for LDA === 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)
SOURCE
899 subroutine drivexc(ixc,order,npts,nspden,usegradient,uselaplacian,usekden,& 900 & rho_updn,exc,vxcrho,nvxcgrho,nvxclrho,nvxctau,ndvxc,nd2vxc, & ! mandatory arguments 901 & grho2_updn,vxcgrho,lrho_updn,vxclrho,tau_updn,vxctau,dvxc,d2vxc, & ! optional arguments 902 & exexch,el_temp,fxcT,hyb_mixing,xc_funcs) ! optional parameters 903 904 !Arguments ------------------------------------ 905 !scalars 906 integer,intent(in) :: ixc,npts,nspden 907 integer,intent(in) :: ndvxc,nd2vxc,nvxcgrho,nvxclrho,nvxctau,order 908 integer,intent(in) :: usegradient,uselaplacian,usekden 909 integer,intent(in),optional :: exexch 910 real(dp),intent(in),optional :: el_temp,hyb_mixing 911 !arrays 912 real(dp),intent(in) :: rho_updn(npts,nspden) 913 real(dp),intent(in),optional :: grho2_updn(npts,(2*nspden-1)*usegradient) 914 real(dp),intent(in),optional :: lrho_updn(npts,nspden*uselaplacian),tau_updn(npts,nspden*usekden) 915 real(dp),intent(out) :: exc(npts),vxcrho(npts,nspden) 916 real(dp),intent(out),optional :: dvxc(npts,ndvxc),d2vxc(npts,nd2vxc),fxcT(:) 917 real(dp),intent(out),optional :: vxcgrho(npts,nvxcgrho),vxclrho(npts,nvxclrho),vxctau(npts,nvxctau) 918 type(libxc_functional_type),intent(inout),optional :: xc_funcs(2) 919 920 !Local variables------------------------------- 921 !scalars 922 integer :: ispden,ixc_from_lib,ixc1,ixc2,ndvxc_x 923 integer :: my_exexch,need_ndvxc,need_nd2vxc,need_nvxcgrho,need_nvxclrho,need_nvxctau 924 integer :: need_gradient,need_laplacian,need_kden,optpbe 925 logical :: has_gradient,has_laplacian,has_kden,libxc_test 926 real(dp) :: alpha,beta,my_hyb_mixing 927 real(dp),parameter :: rsfac=0.6203504908994000e0_dp 928 character(len=500) :: message 929 !arrays 930 real(dp),allocatable :: exci_rpa(:),rhotot(:),rspts(:),vxci_rpa(:,:),zeta(:) 931 real(dp),allocatable :: exc_c(:),exc_x(:),vxcrho_c(:,:),vxcrho_x(:,:) 932 real(dp),allocatable :: d2vxc_c(:,:),d2vxc_x(:,:),dvxc_c(:,:),dvxc_x(:,:) 933 real(dp),allocatable :: vxcgrho_x(:,:) 934 type(libxc_functional_type) :: xc_funcs_vwn3(2),xc_funcs_lyp(2) 935 936 ! ************************************************************************* 937 938 !optional arguments 939 my_exexch=0;if(present(exexch)) my_exexch=exexch 940 my_hyb_mixing=0 941 if (ixc==41) my_hyb_mixing=quarter 942 if (ixc==42) my_hyb_mixing=third 943 if (present(hyb_mixing)) my_hyb_mixing=hyb_mixing 944 945 ! ================================================= 946 ! == Compatibility tests == 947 ! ================================================= 948 949 !Check libXC initialization 950 if (ixc<0 .or. ixc==1402) then 951 libxc_test=libxc_functionals_check(stop_if_error=.true.) 952 end if 953 954 ! Check libXC consistency between ixc passed in input 955 ! and the one used to initialize the libXC library 956 if (ixc<0) then 957 ixc_from_lib=libxc_functionals_ixc() 958 if (present(xc_funcs)) then 959 ixc_from_lib=libxc_functionals_ixc(xc_functionals=xc_funcs) 960 else 961 ixc_from_lib=libxc_functionals_ixc() 962 end if 963 if (ixc/=ixc_from_lib) then 964 write(message, '(a,i0,2a,i0)')& 965 & 'The value of ixc specified in input, ixc = ',ixc,ch10,& 966 & 'differs from the one used to initialize the functional ',ixc_from_lib 967 ABI_BUG(message) 968 end if 969 end if 970 971 !Check value of order 972 if( (order<1.and.order/=-2).or.order>4)then 973 write(message, '(a,i0)' )& 974 & 'The only allowed values for order are 1, 2, -2 or 3, while it is found to be ',order 975 ABI_BUG(message) 976 end if 977 978 !Determine quantities available in input arguments 979 has_gradient=.false.;has_laplacian=.false.;has_kden=.false. 980 if (usegradient==1) then 981 if (.not.present(grho2_updn)) then 982 ABI_BUG('missing grho2_updn argument!') 983 end if 984 if (nvxcgrho>0) then 985 if (.not.present(vxcgrho)) then 986 ABI_BUG('missing vxcgrho argument!') 987 end if 988 has_gradient=.true. 989 end if 990 else if (nvxcgrho>0) then 991 ABI_BUG('nvxcgrho>0 and usegradient=0!') 992 end if 993 if (uselaplacian==1) then 994 if (.not.present(lrho_updn)) then 995 ABI_BUG('missing lrho_updn argument!') 996 end if 997 if (nvxclrho>0) then 998 if (.not.present(vxclrho)) then 999 ABI_BUG('missing vxclrho argument!') 1000 end if 1001 has_laplacian=.true. 1002 end if 1003 else if (nvxclrho>0) then 1004 ABI_BUG('nvxclrho>0 and uselaplacian=0!') 1005 end if 1006 if (usekden==1) then 1007 if (.not.present(tau_updn)) then 1008 ABI_BUG('missing tau_updn argument!') 1009 end if 1010 if (nvxctau>0) then 1011 if (.not.present(vxctau)) then 1012 ABI_BUG('missing vxctau argument!') 1013 end if 1014 has_kden=.true. 1015 end if 1016 else if (nvxctau>0) then 1017 ABI_BUG('nvxctau>0 and usekden=0!') 1018 end if 1019 if (abs(order)>=2) then 1020 if (.not.present(dvxc)) then 1021 message='order>=2 needs argument dvxc!' 1022 ABI_BUG(message) 1023 else if (ndvxc==0) then 1024 message='order>=2 needs ndvxc>0!' 1025 ABI_BUG(message) 1026 end if 1027 end if 1028 if (abs(order)>=3) then 1029 if (.not.present(d2vxc)) then 1030 message='order>=3 needs argument d2vxc!' 1031 ABI_BUG(message) 1032 else if (nd2vxc==0) then 1033 message='order>=3 needs nd2vxc>0!' 1034 ABI_BUG(message) 1035 end if 1036 end if 1037 1038 !Determine quantities needed by XC functional 1039 if (present(xc_funcs)) then 1040 call size_dvxc(ixc,order,nspden,usegradient=need_gradient,& 1041 & uselaplacian=need_laplacian,usekden=need_kden,& 1042 & nvxcgrho=need_nvxcgrho,nvxclrho=need_nvxclrho,& 1043 & nvxctau=need_nvxctau,ndvxc=need_ndvxc,nd2vxc=need_nd2vxc,& 1044 & xc_funcs=xc_funcs) 1045 else 1046 call size_dvxc(ixc,order,nspden,usegradient=need_gradient,& 1047 & uselaplacian=need_laplacian,usekden=need_kden,& 1048 & nvxcgrho=need_nvxcgrho,nvxclrho=need_nvxclrho,& 1049 & nvxctau=need_nvxctau,ndvxc=need_ndvxc,nd2vxc=need_nd2vxc) 1050 end if 1051 if ((has_gradient.and.need_gradient>usegradient).or.& 1052 & (has_laplacian.and.need_laplacian>uselaplacian).or.& 1053 & (has_kden.and.need_kden>usekden)) then 1054 write(message, '(3a)' )& 1055 & 'one of the arguments usegradient/uselaplacian/usesekden',ch10,& 1056 & 'doesnt match the requirements of the XC functional!' 1057 ABI_BUG(message) 1058 end if 1059 if ((has_gradient.and.need_nvxcgrho>nvxcgrho).or.& 1060 & (has_laplacian.and.need_nvxclrho>nvxclrho).or.& 1061 & (has_kden.and.need_nvxctau>nvxctau).or.& 1062 & need_ndvxc>ndvxc.or.need_nd2vxc>nd2vxc) then 1063 write(message, '(3a)' )& 1064 & 'one of the arguments nvxcgrho/nvxclrho/nvxctau/ndvxc/nd2vxc',ch10,& 1065 & 'doesnt match the requirements of the XC functional!' 1066 ABI_BUG(message) 1067 end if 1068 !Deactivate this test because, in case of mGGA, we can output derivatives involving 1069 ! the density and its gradient. Derivatives involving tau or Laplacian will not be output. 1070 !if (abs(order)>1.and.ixc<0.and.(need_laplacian==1.or.need_kden==1)) then 1071 ! message='Derivatives of XC potential are not available in mGGA!' 1072 ! ABI_BUG(message) 1073 !end if 1074 1075 !Check other optional arguments 1076 if (my_exexch/=0.and.usegradient==0) then 1077 message='exexch argument only valid for GGA!' 1078 ABI_BUG(message) 1079 end if 1080 if(ixc==50) then 1081 if(.not.(present(el_temp)).or.(.not.present(fxcT)))then 1082 message = 'el_temp or fxcT is not present but are needed for IIT XC functional.' 1083 ABI_BUG(message) 1084 end if 1085 if (size(fxcT)/=npts) then 1086 ABI_BUG('fxcT size must be npts!') 1087 end if 1088 end if 1089 1090 ! ================================================= 1091 ! == Intermediate quantities computation == 1092 ! ================================================= 1093 1094 !If needed, compute rhotot and rs 1095 if (ixc==1.or.ixc==2.or.ixc==3.or.ixc==4.or.ixc==5.or.& 1096 & ixc==6.or.ixc==21.or.ixc==22.or.ixc==50) then 1097 ABI_MALLOC(rhotot,(npts)) 1098 ABI_MALLOC(rspts,(npts)) 1099 if(nspden==1)then 1100 rhotot(:)=two*rho_updn(:,1) 1101 else 1102 rhotot(:)=rho_updn(:,1)+rho_updn(:,2) 1103 end if 1104 call invcb(rhotot,rspts,npts) 1105 rspts(:)=rsfac*rspts(:) 1106 end if 1107 1108 !If needed, compute zeta 1109 if (ixc==1.or.ixc==21.or.ixc==22) then 1110 ABI_MALLOC(zeta,(npts)) 1111 if(nspden==1)then 1112 zeta(:)=zero 1113 else 1114 zeta(:)=two*rho_updn(:,1)/rhotot(:)-one 1115 end if 1116 end if 1117 1118 ! ================================================= 1119 ! == XC energy, potentiel, ... computation == 1120 ! ================================================= 1121 1122 !>>>>> No exchange-correlation 1123 if (ixc==0.or.ixc==40) then 1124 exc=zero ; vxcrho=zero 1125 if (present(dvxc).and.ndvxc>0) dvxc(:,:)=zero 1126 if (present(d2vxc).and.nd2vxc>0) d2vxc(:,:)=zero 1127 if (present(vxcgrho).and.nvxcgrho>0) vxcgrho(:,:)=zero 1128 if (present(vxclrho).and.nvxclrho>0) vxclrho(:,:)=zero 1129 if (present(vxctau).and.nvxctau>0) vxctau(:,:)=zero 1130 1131 !>>>>> New Teter fit (4/93) to Ceperley-Alder data, with spin-pol option 1132 else if (ixc==1 .or. ixc==21 .or. ixc==22) then 1133 ! new Teter fit (4/93) to Ceperley-Alder data, with spin-pol option 1134 if (order**2 <= 1) then 1135 call xcspol(exc,npts,nspden,order,rspts,vxcrho,zeta,ndvxc) 1136 else 1137 call xcspol(exc,npts,nspden,order,rspts,vxcrho,zeta,ndvxc,dvxc) 1138 end if 1139 1140 !>>>>> Perdew-Zunger fit to Ceperly-Alder data (no spin-pol) 1141 else if (ixc==2) then 1142 if (order**2 <= 1) then 1143 call xcpzca(exc,npts,order,rhotot,rspts,vxcrho(:,1)) 1144 else 1145 call xcpzca(exc,npts,order,rhotot,rspts,vxcrho(:,1),dvxc) 1146 end if 1147 1148 !>>>>> Teter fit (4/91) to Ceperley-Alder values (no spin-pol) 1149 else if (ixc==3) then 1150 if (order**2 <= 1) then 1151 call xctetr(exc,npts,order,rhotot,rspts,vxcrho(:,1)) 1152 else if (order == 2) then 1153 call xctetr(exc,npts,order,rhotot,rspts,vxcrho(:,1),dvxc=dvxc) 1154 else if (order == 3) then 1155 call xctetr(exc,npts,order,rhotot,rspts,vxcrho(:,1),d2vxc=d2vxc,dvxc=dvxc) 1156 end if 1157 1158 !>>>>> Wigner xc (no spin-pol) 1159 else if (ixc==4) then 1160 if (order**2 <= 1) then 1161 call xcwign(exc,npts,order,rspts,vxcrho(:,1)) 1162 else 1163 call xcwign(exc,npts,order,rspts,vxcrho(:,1),dvxc) 1164 end if 1165 1166 !>>>>> Hedin-Lundqvist xc (no spin-pol) 1167 else if (ixc==5) then 1168 if (order**2 <= 1) then 1169 call xchelu(exc,npts,order,rspts,vxcrho(:,1)) 1170 else 1171 call xchelu(exc,npts,order,rspts,vxcrho(:,1),dvxc) 1172 end if 1173 1174 !>>>>> X-alpha (no spin-pol) 1175 else if (ixc==6) then 1176 if (order**2 <= 1) then 1177 call xcxalp(exc,npts,order,rspts,vxcrho(:,1)) 1178 else 1179 call xcxalp(exc,npts,order,rspts,vxcrho(:,1),dvxc) 1180 end if 1181 1182 !>>>>> PBE and alternatives 1183 else if (((ixc>=7.and.ixc<=15).or.(ixc>=23.and.ixc<=24)).and.ixc/=10.and.ixc/=13) then 1184 ! Perdew-Wang LSD is coded in Perdew-Burke-Ernzerhof GGA, with optpbe=1 1185 if(ixc==7)optpbe=1 1186 ! x-only part of Perdew-Wang 1187 if(ixc==8)optpbe=-1 1188 ! Exchange + RPA correlation from Perdew-Wang 1189 if(ixc==9)optpbe=3 1190 ! Perdew-Burke-Ernzerhof GGA 1191 if(ixc==11)optpbe=2 1192 ! x-only part of PBE 1193 if(ixc==12)optpbe=-2 1194 ! C09x exchange of V. R. Cooper 1195 if(ixc==24)optpbe=-4 1196 ! revPBE of Zhang and Yang 1197 if(ixc==14)optpbe=5 1198 ! RPBE of Hammer, Hansen and Norskov 1199 if(ixc==15)optpbe=6 1200 ! Wu and Cohen 1201 if(ixc==23)optpbe=7 1202 if (ixc >=7.and.ixc<=9) then 1203 if (order**2 <= 1) then 1204 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc) 1205 else if (order /=3) then 1206 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,dvxci=dvxc) 1207 else if (order ==3) then 1208 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,d2vxci=d2vxc,dvxci=dvxc) 1209 end if 1210 else if ((ixc >= 11 .and. ixc <= 15) .or. (ixc>=23 .and. ixc<=24)) then 1211 if (order**2 <= 1) then 1212 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,& 1213 & dvxcdgr=vxcgrho,exexch=my_exexch,grho2_updn=grho2_updn) 1214 else if (order /=3) then 1215 if(ixc==12 .or. ixc==24)then 1216 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,& 1217 & dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 1218 else if(ixc/=12 .or. ixc/=24) then 1219 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,& 1220 & dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 1221 end if 1222 else if (order ==3) then 1223 if(ixc==12 .or. ixc==24)then 1224 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,& 1225 & d2vxci=d2vxc,dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 1226 else if(ixc/=12 .or. ixc/=24) then 1227 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,& 1228 & d2vxci=d2vxc,dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 1229 end if 1230 end if 1231 end if 1232 1233 !>>>>> RPA correlation from Perdew-Wang 1234 else if (ixc==10) then 1235 if (order**2 <= 1) then 1236 ABI_MALLOC(exci_rpa,(npts)) 1237 ABI_MALLOC(vxci_rpa,(npts,2)) 1238 optpbe=3 1239 call xcpbe(exci_rpa,npts,nspden,optpbe,order,rho_updn,vxci_rpa,ndvxc,nd2vxc) 1240 optpbe=1 1241 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc) 1242 exc(:)=exc(:)-exci_rpa(:) 1243 ! PMA: second index of vxcrho is nspden while that of rpa is 2 they can mismatch 1244 vxcrho(:,1:min(nspden,2))=vxcrho(:,1:min(nspden,2))-vxci_rpa(:,1:min(nspden,2)) 1245 ABI_FREE(exci_rpa) 1246 ABI_FREE(vxci_rpa) 1247 else if (order /=3) then 1248 ABI_MALLOC(exci_rpa,(npts)) 1249 ABI_MALLOC(vxci_rpa,(npts,2)) 1250 optpbe=3 1251 call xcpbe(exci_rpa,npts,nspden,optpbe,order,rho_updn,vxci_rpa,ndvxc,nd2vxc,dvxci=dvxc) 1252 optpbe=1 1253 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,dvxci=dvxc) 1254 exc(:)=exc(:)-exci_rpa(:) 1255 vxcrho(:,:)=vxcrho(:,:)-vxci_rpa(:,:) 1256 ABI_FREE(exci_rpa) 1257 ABI_FREE(vxci_rpa) 1258 else if (order ==3) then 1259 ABI_MALLOC(exci_rpa,(npts)) 1260 ABI_MALLOC(vxci_rpa,(npts,2)) 1261 optpbe=3 1262 call xcpbe(exci_rpa,npts,nspden,optpbe,order,rho_updn,vxci_rpa,ndvxc,nd2vxc,& 1263 & d2vxci=d2vxc,dvxci=dvxc) 1264 optpbe=1 1265 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,& 1266 & d2vxci=d2vxc,dvxci=dvxc) 1267 exc(:)=exc(:)-exci_rpa(:) 1268 vxcrho(:,:)=vxcrho(:,:)-vxci_rpa(:,:) 1269 ABI_FREE(exci_rpa) 1270 ABI_FREE(vxci_rpa) 1271 end if 1272 1273 !>>>>> LDA xc energy like ixc==7, and Leeuwen-Baerends GGA xc potential 1274 else if(ixc==13) then 1275 if (order**2 <= 1) then 1276 optpbe=1 1277 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc) 1278 call xclb(grho2_updn,npts,nspden,rho_updn,vxcrho) 1279 else if (order /=3) then 1280 optpbe=1 1281 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,dvxci=dvxc) 1282 call xclb(grho2_updn,npts,nspden,rho_updn,vxcrho) 1283 else if (order ==3) then 1284 optpbe=1 1285 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,d2vxci=d2vxc,dvxci=dvxc) 1286 call xclb(grho2_updn,npts,nspden,rho_updn,vxcrho) 1287 end if 1288 1289 !>>>>> HTCH93, HTCH120, HTCH107, HTCH147 1290 else if(ixc==16 .or. ixc==17 .or. ixc==26 .or. ixc==27) then 1291 call xchcth(vxcgrho,exc,grho2_updn,ixc,npts,nspden,order,rho_updn,vxcrho) 1292 1293 !>>>>> Only for test purpose (test various part of MGGA implementation) 1294 else if(ixc==31 .or. ixc==32 .or. ixc==33 .or. ixc==34 .or. ixc==35) then 1295 exc(:)=zero ; vxcrho(:,:)=zero 1296 if (present(vxcgrho).and.nvxcgrho>0) vxcgrho(:,:)=zero 1297 if (present(vxclrho).and.nvxclrho>0) vxclrho(:,:)=zero 1298 if (present(vxctau).and.nvxctau>0) vxctau(:,:)=zero 1299 if (present(dvxc).and.ndvxc>0) dvxc(:,:)=zero 1300 if (present(d2vxc).and.nd2vxc>0) d2vxc(:,:)=zero 1301 1302 !>>>>> Perdew-Wang LSD is coded in Perdew-Burke-Ernzerhof GGA, with optpbe=1 1303 optpbe=1 1304 select case(ixc) 1305 case (31) 1306 alpha=1.00d0-(1.00d0/1.01d0) 1307 ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc) 1308 call xcpbe(exc,npts,nspden,optpbe,1,rho_updn,vxcrho,0,0) 1309 if (nspden==1) then 1310 exc(:)=exc(:)+alpha*tau_updn(:,1)/rho_updn(:,1) 1311 else 1312 ! It should be : exc_tot= exc_spin up + exc_spin down = 2*exc_spin up 1313 ! but this applies to tau and rho (so it cancels) 1314 do ispden=1,nspden 1315 exc(:)=exc(:)+alpha*tau_updn(:,ispden)/(rho_updn(:,1)+rho_updn(:,2)) 1316 end do 1317 end if 1318 vxctau(:,:)=alpha 1319 case (32) 1320 alpha=0.01d0 1321 ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc) 1322 call xcpbe(exc,npts,nspden,optpbe,1,rho_updn,vxcrho,0,0) 1323 if (nspden==1) then 1324 exc(:)=exc(:)+2.0d0*alpha*lrho_updn(:,1) 1325 vxcrho(:,1) =vxcrho(:,1)+2.0d0*alpha*lrho_updn(:,1) 1326 vxclrho(:,1)=alpha*2.0d0*rho_updn(:,1) 1327 else 1328 do ispden=1,nspden 1329 exc(:)=exc(:)+alpha*lrho_updn(:,ispden) 1330 vxcrho(:,ispden) =vxcrho(:,ispden)+alpha*(lrho_updn(:,1)+lrho_updn(:,2)) 1331 vxclrho(:,ispden)=alpha*(rho_updn(:,1)+rho_updn(:,2)) 1332 end do 1333 end if 1334 case (33) 1335 alpha=-0.010d0 1336 ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc) 1337 call xcpbe(exc,npts,nspden,optpbe,1,rho_updn,vxcrho,0,0) 1338 if (nspden==1) then 1339 ! it should be : exc_tot= exc_spin up + exc_spin down = 2*exc_spin up but this applies to grho2 and rho 1340 ! (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) 1341 exc(:)=exc(:)+alpha*2.0d0*grho2_updn(:,1)/rho_updn(:,1) 1342 if(nvxcgrho==2)vxcgrho(:,1:2)=2.0d0*alpha 1343 if(nvxcgrho==3)vxcgrho(:,3)=2.0d0*alpha 1344 else 1345 exc(:)=exc(:)+alpha*grho2_updn(:,3)/(rho_updn(:,1)+rho_updn(:,2)) 1346 if(nvxcgrho==2)vxcgrho(:,1:2)=2.0d0*alpha 1347 if(nvxcgrho==3)vxcgrho(:,3)=2.0d0*alpha 1348 end if 1349 case (34) 1350 alpha=-0.010d0 1351 ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc) 1352 call xcpbe(exc,npts,nspden,optpbe,1,rho_updn,vxcrho,0,0) 1353 if (nspden==1) then 1354 exc(:)=exc(:)+16.0d0*alpha*tau_updn(:,1) 1355 vxcrho(:,1)=vxcrho(:,1)+16.0d0*alpha*tau_updn(:,1) 1356 vxctau(:,1)=16.0d0*alpha*rho_updn(:,1) 1357 else 1358 do ispden=1,nspden 1359 exc(:)=exc(:)+8.0d0*alpha*tau_updn(:,ispden) 1360 vxcrho(:,ispden)=vxcrho(:,ispden)+8.0d0*alpha*(tau_updn(:,1)+tau_updn(:,2)) 1361 vxctau(:,ispden)=8.0d0*alpha*(rho_updn(:,1)+rho_updn(:,2)) 1362 end do 1363 end if 1364 case (35) 1365 alpha=0.01d0 ; beta=1.00d0-(1.00d0/1.01d0) 1366 ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc) 1367 call xcpbe(exc,npts,nspden,optpbe,1,rho_updn,vxcrho,0,0) 1368 if (nspden==1) then 1369 exc(:)=exc(:)+2.0d0*alpha*lrho_updn(:,1)+beta*tau_updn(:,1)/rho_updn(:,1) 1370 vxcrho(:,1) =vxcrho(:,1)+2.0d0*alpha*lrho_updn(:,1) 1371 vxclrho(:,1)=alpha*2.0d0*rho_updn(:,1) 1372 else 1373 do ispden=1,nspden 1374 exc(:)=exc(:)+alpha*lrho_updn(:,ispden) & 1375 & +beta*tau_updn(:,ispden)/(rho_updn(:,1)+rho_updn(:,2)) 1376 vxcrho(:,ispden) =vxcrho(:,ispden)+alpha*(lrho_updn(:,1)+lrho_updn(:,2)) 1377 vxclrho(:,ispden)=alpha*(rho_updn(:,1)+rho_updn(:,2)) 1378 end do 1379 end if 1380 vxctau(:,:)=beta 1381 end select 1382 1383 !>>>>> Hybrid PBE0 (1/4 and 1/3) 1384 else if(ixc>=41.and.ixc<=42) then 1385 ! Requires to evaluate exchange-correlation with PBE (optpbe=2) 1386 ! minus hyb_mixing*exchange with PBE (optpbe=-2) 1387 ndvxc_x=8 1388 ABI_MALLOC(exc_x,(npts)) 1389 ABI_MALLOC(vxcrho_x,(npts,nspden)) 1390 ABI_MALLOC(vxcgrho_x,(npts,nvxcgrho)) 1391 exc_x=zero;vxcrho_x=zero;vxcgrho_x=zero 1392 if (order**2 <= 1) then 1393 optpbe=2 !PBE exchange correlation 1394 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,& 1395 & dvxcdgr=vxcgrho,exexch=my_exexch,grho2_updn=grho2_updn) 1396 optpbe=-2 !PBE exchange-only 1397 call xcpbe(exc_x,npts,nspden,optpbe,order,rho_updn,vxcrho_x,ndvxc,nd2vxc,& 1398 & dvxcdgr=vxcgrho_x,exexch=my_exexch,grho2_updn=grho2_updn) 1399 exc=exc-exc_x*my_hyb_mixing 1400 vxcrho=vxcrho-vxcrho_x*my_hyb_mixing 1401 vxcgrho=vxcgrho-vxcgrho_x*my_hyb_mixing 1402 else if (order /=3) then 1403 ABI_MALLOC(dvxc_x,(npts,ndvxc_x)) 1404 optpbe=2 !PBE exchange correlation 1405 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,& 1406 dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 1407 optpbe=-2 !PBE exchange-only 1408 call xcpbe(exc_x,npts,nspden,optpbe,order,rho_updn,vxcrho_x,ndvxc_x,nd2vxc,& 1409 & dvxcdgr=vxcgrho_x,dvxci=dvxc_x,grho2_updn=grho2_updn) 1410 exc=exc-exc_x*my_hyb_mixing 1411 vxcrho=vxcrho-vxcrho_x*my_hyb_mixing 1412 vxcgrho=vxcgrho-vxcgrho_x*my_hyb_mixing 1413 dvxc(:,1:ndvxc_x)=dvxc(:,1:ndvxc_x)-dvxc_x(:,1:ndvxc_x)*my_hyb_mixing 1414 ABI_FREE(dvxc_x) 1415 else if (order ==3) then 1416 ! The size of exchange-correlation with PBE (optpbe=2) 1417 ! is the one which defines the size for ndvxc. 1418 ABI_MALLOC(dvxc_x,(npts,ndvxc_x)) 1419 ABI_MALLOC(d2vxc_x,(npts,nd2vxc)) 1420 optpbe=2 !PBE exchange correlation 1421 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,nd2vxc,& 1422 & d2vxci=d2vxc,dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 1423 optpbe=-2 !PBE exchange-only 1424 call xcpbe(exc_x,npts,nspden,optpbe,order,rho_updn,vxcrho_x,ndvxc_x,nd2vxc,& 1425 & d2vxci=d2vxc_x,dvxcdgr=vxcgrho_x,dvxci=dvxc_x,grho2_updn=grho2_updn) 1426 exc=exc-exc_x*my_hyb_mixing 1427 vxcrho=vxcrho-vxcrho_x*my_hyb_mixing 1428 vxcgrho=vxcgrho-vxcgrho_x*my_hyb_mixing 1429 d2vxc=d2vxc-d2vxc_x*my_hyb_mixing 1430 dvxc(:,1:ndvxc_x)=dvxc(:,1:ndvxc_x)-dvxc_x(:,1:ndvxc_x)*my_hyb_mixing 1431 ABI_FREE(dvxc_x) 1432 ABI_FREE(d2vxc_x) 1433 end if 1434 ABI_FREE(exc_x) 1435 ABI_FREE(vxcrho_x) 1436 ABI_FREE(vxcgrho_x) 1437 1438 !>>>>> Ichimaru,Iyetomi,Tanaka, XC at finite temp (e- gaz) 1439 else if (ixc==50) then 1440 if (order**2 <= 1) then 1441 call xciit(exc,fxcT,npts,order,rspts,el_temp,vxcrho(:,1)) 1442 else 1443 call xciit(exc,fxcT,npts,order,rspts,el_temp,vxcrho(:,1),dvxc) 1444 end if 1445 1446 !>>>>> GGA counterpart of the B3LYP functional 1447 else if(ixc==1402000) then 1448 ! Requires to evaluate exchange-correlation 1449 ! with 5/4 B3LYP - 1/4 B3LYPc, where 1450 ! B3LYPc = (0.19 Ec VWN3 + 0.81 Ec LYP) 1451 1452 ! First evaluate B3LYP. 1453 if(present(xc_funcs))then 1454 if (abs(order)==1) then 1455 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1456 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,xc_functionals=xc_funcs) 1457 else if (abs(order)==2) then 1458 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1459 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,xc_functionals=xc_funcs) 1460 else 1461 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1462 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,d2vxc=d2vxc,xc_functionals=xc_funcs) 1463 end if 1464 else 1465 if (abs(order)==1) then 1466 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1467 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho) 1468 else if (abs(order)==2) then 1469 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1470 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc) 1471 else 1472 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1473 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,d2vxc=d2vxc) 1474 end if 1475 end if 1476 1477 ! Then renormalize B3LYP and subtract VWN3 contribution 1478 ABI_MALLOC(exc_c,(npts)) 1479 ABI_MALLOC(vxcrho_c,(npts,nspden)) 1480 if(order**2>1)then 1481 ABI_MALLOC(dvxc_c,(npts,ndvxc)) 1482 end if 1483 if(order**2>4)then 1484 ABI_MALLOC(d2vxc_c,(npts,nd2vxc)) 1485 end if 1486 exc_c=zero;vxcrho_c=zero 1487 call libxc_functionals_init(-30,nspden,xc_functionals=xc_funcs_vwn3) 1488 if (order**2 <= 1) then 1489 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 1490 & vxcrho_c,xc_functionals=xc_funcs_vwn3) 1491 elseif (order**2 <= 4) then 1492 dvxc_c=zero 1493 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 1494 & vxcrho_c,dvxc=dvxc_c,xc_functionals=xc_funcs_vwn3) 1495 else 1496 dvxc_c=zero 1497 d2vxc_c=zero 1498 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 1499 & vxcrho_c,dvxc=dvxc_c,d2vxc=d2vxc,xc_functionals=xc_funcs_vwn3) 1500 end if 1501 exc=1.25d0*exc-quarter*0.19d0*exc_c 1502 vxcrho=1.25d0*vxcrho-quarter*0.19d0*vxcrho_c 1503 if(order**2>1)dvxc=1.25d0*dvxc-quarter*0.19d0*dvxc_c 1504 if(order**2>4)d2vxc=1.25d0*d2vxc-quarter*0.19d0*d2vxc_c 1505 call libxc_functionals_end(xc_functionals=xc_funcs_vwn3) 1506 1507 ! Then subtract LYP contribution 1508 call libxc_functionals_init(-131,nspden,xc_functionals=xc_funcs_lyp) 1509 if (order**2 <= 1) then 1510 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 1511 & vxcrho_c,grho2=grho2_updn,vxcgr=vxcgrho,xc_functionals=xc_funcs_lyp) 1512 elseif (order**2 <= 4) then 1513 dvxc_c=zero 1514 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 1515 & vxcrho_c,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc_c,xc_functionals=xc_funcs_lyp) 1516 else 1517 dvxc_c=zero 1518 d2vxc_c=zero 1519 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 1520 & vxcrho_c,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc_c,d2vxc=d2vxc,xc_functionals=xc_funcs_lyp) 1521 end if 1522 exc=exc-quarter*0.81d0*exc_c 1523 vxcrho=vxcrho-quarter*0.81d0*vxcrho_c 1524 if(order**2>1)dvxc=dvxc-quarter*0.81d0*dvxc_c 1525 if(order**2>4)d2vxc=d2vxc-quarter*0.81d0*d2vxc_c 1526 call libxc_functionals_end(xc_functionals=xc_funcs_lyp) 1527 1528 ABI_FREE(exc_c) 1529 ABI_FREE(vxcrho_c) 1530 if(allocated(dvxc_c))then 1531 ABI_FREE(dvxc_c) 1532 end if 1533 if(allocated(d2vxc_c))then 1534 ABI_FREE(d2vxc_c) 1535 end if 1536 1537 !>>>>> All libXC functionals 1538 else if( ixc<0 ) then 1539 1540 ! ===== meta-GGA ===== 1541 if (need_laplacian==1.or.need_kden==1) then 1542 if (need_laplacian==1.and.need_kden==1) then 1543 if (abs(order)<=1) then 1544 if (present(xc_funcs)) then 1545 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1546 & grho2=grho2_updn,vxcgr=vxcgrho,& 1547 & lrho=lrho_updn,vxclrho=vxclrho,& 1548 & tau=tau_updn,vxctau=vxctau,& 1549 & xc_functionals=xc_funcs) 1550 else 1551 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1552 & grho2=grho2_updn,vxcgr=vxcgrho,& 1553 & lrho=lrho_updn,vxclrho=vxclrho,& 1554 & tau=tau_updn,vxctau=vxctau) 1555 end if 1556 else 1557 if (present(xc_funcs)) then 1558 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1559 & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,& 1560 & lrho=lrho_updn,vxclrho=vxclrho,& 1561 & tau=tau_updn,vxctau=vxctau,& 1562 & xc_functionals=xc_funcs) 1563 else 1564 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1565 & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,& 1566 & lrho=lrho_updn,vxclrho=vxclrho,& 1567 & tau=tau_updn,vxctau=vxctau) 1568 end if 1569 end if 1570 else if (need_laplacian==1) then 1571 if (abs(order)<=1) then 1572 if (present(xc_funcs)) then 1573 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1574 & grho2=grho2_updn,vxcgr=vxcgrho,& 1575 & lrho=lrho_updn,vxclrho=vxclrho,& 1576 & xc_functionals=xc_funcs) 1577 else 1578 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1579 & grho2=grho2_updn,vxcgr=vxcgrho,& 1580 & lrho=lrho_updn,vxclrho=vxclrho) 1581 end if 1582 else 1583 if (present(xc_funcs)) then 1584 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1585 & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,& 1586 & lrho=lrho_updn,vxclrho=vxclrho,& 1587 & xc_functionals=xc_funcs) 1588 else 1589 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1590 & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,& 1591 & lrho=lrho_updn,vxclrho=vxclrho) 1592 end if 1593 end if 1594 else if (need_kden==1) then 1595 if (abs(order)<=1) then 1596 if (present(xc_funcs)) then 1597 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1598 & grho2=grho2_updn,vxcgr=vxcgrho,& 1599 & tau=tau_updn,vxctau=vxctau,& 1600 & xc_functionals=xc_funcs) 1601 else 1602 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1603 & grho2=grho2_updn,vxcgr=vxcgrho,& 1604 & tau=tau_updn,vxctau=vxctau) 1605 end if 1606 else 1607 if (present(xc_funcs)) then 1608 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1609 & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,& 1610 & tau=tau_updn,vxctau=vxctau,& 1611 & xc_functionals=xc_funcs) 1612 else 1613 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1614 & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,& 1615 & tau=tau_updn,vxctau=vxctau) 1616 end if 1617 end if 1618 end if 1619 !Some meta-GGAs can only be used with a LDA correlation (see doc) 1620 ixc1=(-ixc)/1000;ixc2=(-ixc)-ixc1*1000 1621 if (ixc1==206 .or. ixc1==207 .or. ixc1==208 .or. ixc1==209 .or. & 1622 & ixc2==206 .or. ixc2==207 .or. ixc2==208 .or. ixc2==209 )then 1623 if (present(vxcgrho)) vxcgrho(:,:)=zero 1624 if (present(vxclrho)) vxclrho(:,:)=zero 1625 if (present(vxctau)) vxctau(:,:)=zero 1626 if (present(dvxc)) dvxc(:,:)=zero 1627 end if 1628 1629 ! ===== GGA ===== 1630 else if (need_gradient==1) then 1631 if (abs(order)<=1) then 1632 if (present(xc_funcs)) then 1633 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1634 & grho2=grho2_updn,vxcgr=vxcgrho,& 1635 & xc_functionals=xc_funcs) 1636 else 1637 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1638 & grho2=grho2_updn,vxcgr=vxcgrho) 1639 end if 1640 else 1641 if (present(xc_funcs)) then 1642 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1643 & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,& 1644 & xc_functionals=xc_funcs) 1645 else 1646 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1647 & grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc) 1648 end if 1649 end if 1650 1651 ! ===== LDA ===== 1652 else 1653 if (abs(order)<=1) then 1654 if (present(xc_funcs)) then 1655 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1656 & xc_functionals=xc_funcs) 1657 else 1658 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho) 1659 end if 1660 else if (abs(order)<=2) then 1661 if (present(xc_funcs)) then 1662 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1663 & dvxc=dvxc,xc_functionals=xc_funcs) 1664 else 1665 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1666 & dvxc=dvxc) 1667 end if 1668 else if (abs(order)<=3) then 1669 if (present(xc_funcs)) then 1670 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1671 & dvxc=dvxc,d2vxc=d2vxc,xc_functionals=xc_funcs) 1672 else 1673 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,vxcrho,& 1674 & dvxc=dvxc,d2vxc=d2vxc) 1675 end if 1676 end if 1677 1678 end if ! mGGA, GGA, LDA 1679 end if ! libXC 1680 1681 ! ================================================= 1682 ! == Finalization == 1683 ! ================================================= 1684 !Deallocate arrays 1685 if(allocated(rhotot)) then 1686 ABI_FREE(rhotot) 1687 end if 1688 if(allocated(rspts)) then 1689 ABI_FREE(rspts) 1690 end if 1691 if(allocated(zeta)) then 1692 ABI_FREE(zeta) 1693 end if 1694 1695 end subroutine drivexc
m_drivexc/echo_xc_name [ Functions ]
[ Top ] [ m_drivexc ] [ Functions ]
NAME
echo_xc_name
FUNCTION
Write to log and output the xc functional which will be used for this dataset
INPUTS
ixc = internal code for xc functional
SOURCE
66 subroutine echo_xc_name (ixc) 67 68 !Arguments ------------------------------- 69 integer, intent(in) :: ixc 70 71 !Local variables ------------------------- 72 integer :: l_citation 73 character(len=500) :: message, citation 74 75 ! ********************************************************************* 76 77 message ='' 78 citation ='' 79 80 !normal case (not libxc) 81 if (ixc >= 0) then 82 83 select case (ixc) 84 case (0) 85 message = 'No xc applied (usually for testing) - ixc=0' 86 citation = '' 87 ! LDA,LSD 88 case (1) 89 message = 'LDA: new Teter (4/93) with spin-polarized option - ixc=1' 90 citation = 'S. Goedecker, M. Teter, J. Huetter, PRB 54, 1703 (1996)' ! [[cite:Goedecker1996]] 91 case (2) 92 message = 'LDA: Perdew-Zunger-Ceperley-Alder - ixc=2' 93 citation = 'J.P.Perdew and A.Zunger, PRB 23, 5048 (1981) ' ! [[cite:Perdew1981]] 94 case (3) 95 message = 'LDA: old Teter (4/91) fit to Ceperley-Alder data - ixc=3' 96 citation = '' 97 case (4) 98 message = 'LDA: Wigner - ixc=4' 99 citation = 'E.P.Wigner, Trans. Faraday Soc. 34, 67 (1938)' ! [[cite:Wigner1938]] 100 case (5) 101 message = 'LDA: Hedin-Lundqvist - ixc=5' 102 citation = 'L.Hedin and B.I.Lundqvist, J. Phys. C4, 2064 (1971)' ! [[cite:Hedin1971]] 103 case (6) 104 message = 'LDA: "X-alpha" xc - ixc=6' 105 citation = 'Slater J. C., Phys. Rev. 81, 385 (1951)' ! [[cite:Slater1951]] 106 case (7) 107 message = 'LDA: Perdew-Wang 92 LSD fit to Ceperley-Alder data - ixc=7' 108 citation = 'J.P.Perdew and Y.Wang, PRB 45, 13244 (1992)' ! [[cite:Perdew1992a]] 109 case (8) 110 message = 'LDA: Perdew-Wang 92 LSD , exchange-only - ixc=8' 111 citation = 'J.P.Perdew and Y.Wang, PRB 45, 13244 (1992)' ! [[cite:Perdew1992a]] 112 case (9) 113 message = 'LDA: Perdew-Wang 92 Ex+Ec_RPA energy - ixc=9' 114 citation = 'J.P.Perdew and Y.Wang, PRB 45, 13244 (1992)' ! [[cite:Perdew1992]] 115 case (10) 116 message = 'LDA: RPA LSD energy (only the energy !!) - ixc=10' 117 citation = '' 118 ! GGA 119 case (11) 120 message = 'GGA: Perdew-Burke-Ernzerhof functional - ixc=11' 121 citation = 'J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996)' ! [[cite:Perdew1996]] 122 case (12) 123 message = 'GGA: x-only Perdew-Burke-Ernzerhof functional - ixc=12' 124 citation = 'J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996)' ! [[cite:Perdew1996]] 125 case (13) 126 message = 'GGA: LDA (ixc==7) energy, and the xc _potential_ is given by van Leeuwen-Baerends GGA - ixc=13' 127 citation = 'R. van Leeuwen and E. J. Baerends PRA 49, 2421 (1994)' ! [[cite:VanLeeuwen1994]] 128 case (14) 129 message = 'GGA: revPBE functional - ixc=14' 130 citation = 'Zhang and Yang, PRL 80, 890 (1998)' ! [[cite:Zhang1998]] 131 case (15) 132 message = 'GGA: RPBE functional - ixc=15' 133 citation = 'Hammer, L. B. Hansen, and J. K. Norskov, PRB 59, 7413 (1999)' ! [[cite:Hammer1999]] 134 case (16) 135 message = 'GGA: HCTH93 functional - ixc=16' 136 citation = 'F.A. Hamprecht, A.J. Cohen, D.J. Tozer, N.C. Handy, JCP 109, 6264 (1998)' ! [[cite:Hamprecht1998]] 137 case (17) 138 message = 'GGA: HCTH120 functional - ixc=17' 139 citation = 'A.D. Boese, N.L. Doltsinis, N.C. Handy, and M. Sprik, JCP 112, 1670 (2000)' ! [[cite:Boese2000]] 140 case (23) 141 message = 'GGA: Wu Cohen functional - ixc=23' 142 citation = 'Z. Wu and R. E. Cohen, PRB 73, 235116 (2006)' ! [[cite:Wu2006]] 143 case (24) 144 message = 'GGA: C09x exchange functional - ixc=24' 145 citation = 'Valentino R. Cooper, PRB 81, 161104(R) (2010)' ! [[cite:Cooper2010]] 146 case (26) 147 message = 'GGA: HCTH147 functional - ixc=26' 148 citation = 'A.D. Boese, N.L. Doltsinis, N.C. Handy, and M. Sprik, JCP 112, 1670 (2000)' ! [[cite:Boese2000]] 149 case (27) 150 message = 'GGA: HCTH407 functional - ixc=27' 151 citation = 'A.D. Boese, and N.C. Handy, JCP 114, 5497 (2001)' ! [[cite:Boese2001]] 152 ! Fermi-Amaldi 153 case (20) 154 message = 'Fermi-Amaldi correction - ixc=20' 155 citation = '' 156 case (21) 157 message = 'Fermi-Amaldi correction with LDA(ixc=1) kernel - ixc=21' 158 citation = '' 159 case (22) 160 message = 'Fermi-Amaldi correction with hybrid BPG kernel - ixc=22' 161 citation = '' 162 case (31) 163 message = 'Meta-GGA fake1 - ixc=31' 164 citation = '' 165 case (32) 166 message = 'Meta-GGA fake2 - ixc=32' 167 citation = '' 168 case (33) 169 message = 'Meta-GGA fake3 - ixc=33' 170 citation = '' 171 case (34) 172 message = 'Meta-GGA fake4 - ixc=34' 173 citation = '' 174 case (35) 175 message = 'Meta-GGA fake5 - ixc=35' 176 citation = '' 177 case (40) 178 message = 'Hartree-Fock with mixing coefficient alpha=1' 179 citation = '' 180 case (41) 181 message = 'PBE0 with alpha=0.25' 182 citation = '' 183 case (42) 184 message = 'modified PBE0 with alpha=0.33' 185 citation = '' 186 case (50) 187 message = 'LDA at finite T Ichimaru-Iyetomy-Tanaka - ixc=50' 188 citation = 'Ichimaru S., Iyetomi H., Tanaka S., Phys. Rep. 149, 91-205 (1987) ' ! [[cite:Ichimaru1987]] 189 case default 190 write(message,'(a,i0)')" echo_xc_name does not know how to handle ixc = ",ixc 191 ABI_WARNING(message) 192 end select 193 194 message = " Exchange-correlation functional for the present dataset will be:" // ch10 & 195 & // " " // trim(message) 196 197 l_citation=len_trim(citation) 198 citation = " Citation for XC functional:" // ch10 // " " // trim(citation) 199 200 call wrtout(ab_out,message,'COLL') 201 call wrtout(std_out,message,'COLL') 202 203 if(l_citation/=0)then 204 call wrtout(ab_out,citation,'COLL') 205 call wrtout(std_out,citation,'COLL') 206 end if 207 208 message =' ' 209 call wrtout(ab_out,message,'COLL') 210 call wrtout(std_out,message,'COLL') 211 212 end if ! end libxc if 213 214 end subroutine echo_xc_name
m_drivexc/has_k3xc [ Functions ]
[ Top ] [ m_drivexc ] [ Functions ]
NAME
has_k3xc
FUNCTION
Given a XC functional (defined by ixc), return TRUE if K3xc (d2Vxc/drho2) is avalaible.
INPUTS
ixc = internal code for xc functional [xc_funcs(2)]= <type(libxc_functional_type)> = optional - libXC set of functionals
OUTPUT
SOURCE
275 logical function has_k3xc(ixc,xc_funcs) 276 277 !Arguments ------------------------------- 278 integer, intent(in) :: ixc 279 type(libxc_functional_type),intent(in),optional :: xc_funcs(2) 280 281 !Local variables ------------------------- 282 283 ! ********************************************************************* 284 285 has_k3xc=.false. 286 287 if (ixc>=0) then 288 has_k3xc=(ixc==0.or.ixc==3.or.(ixc>=7.and.ixc<=15).or. & 289 & ixc==23.or.ixc==24.or.ixc==41.or.ixc==42.or.ixc==1402000) 290 else if (ixc==-406.or.ixc==-427.or.ixc==-428.or.ixc==-456)then 291 has_k3xc=.false. 292 else ! ixc<0 and not one of the allowed hybrids 293 if (present(xc_funcs)) then 294 has_k3xc=libxc_functionals_has_k3xc(xc_funcs) 295 else 296 has_k3xc=libxc_functionals_has_k3xc() 297 end if 298 end if 299 300 end function has_k3xc
m_drivexc/has_kxc [ Functions ]
[ Top ] [ m_drivexc ] [ Functions ]
NAME
has_kxc
FUNCTION
Given a XC functional (defined by ixc), return TRUE if Kxc (dVxc/drho) is avalaible.
INPUTS
ixc = internal code for xc functional [xc_funcs(2)]= <type(libxc_functional_type)> = optional - libXC set of functionals
OUTPUT
SOURCE
232 logical function has_kxc(ixc,xc_funcs) 233 234 !Arguments ------------------------------- 235 integer, intent(in) :: ixc 236 type(libxc_functional_type),intent(in),optional :: xc_funcs(2) 237 238 !Local variables ------------------------- 239 240 ! ********************************************************************* 241 242 has_kxc=.false. 243 244 if (ixc>=0) then 245 has_kxc=(ixc/=16.and.ixc/=17.and.ixc/=26.and.ixc/=27) 246 else if (ixc==-406.or.ixc==-427.or.ixc==-428.or.ixc==-456)then 247 has_kxc=.true. 248 else ! ixc<0 and not one of the allowed hybrids 249 if (present(xc_funcs)) then 250 has_kxc=libxc_functionals_has_kxc(xc_funcs) 251 else 252 has_kxc=libxc_functionals_has_kxc() 253 end if 254 end if 255 256 end function has_kxc
m_drivexc/mkdenpos [ Functions ]
[ Top ] [ m_drivexc ] [ Functions ]
NAME
mkdenpos
FUNCTION
Make a density positive everywhere: when the density (or spin-density) is smaller than xc_denpos, set it to the value of xc_denpos
INPUTS
nfft=(effective) number of FFT grid points (for this processor)
nspden=number of spin-density components (max. 2)
option=0 if density rhonow is stored as (up,dn)
1 if density rhonow is stored as (up+dn,up)
Active only when nspden=2
xc_denpos= lowest allowed density (usually for the computation of the XC functionals)
OUTPUT
(see side effects)
SIDE EFFECTS
Input/output
iwarn=At input: iwarn=0 a warning will be printed when rho is negative
iwarn>0 no warning will be printed out
At output: iwarn is increased by 1
rhonow(nfft,nspden)=electron (spin)-density in real space,
either on the unshifted grid (if ishift==0,
then equal to rhor),or on the shifted grid
NOTES
At this stage, rhonow(:,1:nspden) contains the density in real space, on the unshifted or shifted grid. Now test for negative densities Note that, ignoring model core charge, as long as boxcut>=2 the shifted density is derivable from the square of a Fourier interpolated charge density => CANNOT go < 0. However, actually can go < 0 to within machine precision; do not print useless warnings in this case, just fix it. Fourier interpolated core charge can go < 0 due to Gibbs oscillations; could avoid this by recomputing the model core charge at the new real space grid points (future work).
SOURCE
681 subroutine mkdenpos(iwarn,nfft,nspden,option,rhonow,xc_denpos) 682 683 !Arguments ------------------------------------ 684 !scalars 685 integer,intent(in) :: nfft,nspden,option 686 integer,intent(inout) :: iwarn 687 real(dp),intent(in) :: xc_denpos 688 !arrays 689 real(dp),intent(inout) :: rhonow(nfft,nspden) 690 691 !Local variables------------------------------- 692 !scalars 693 integer :: ifft,ispden,numneg 694 real(dp) :: rhotmp,worst 695 character(len=600) :: message 696 !arrays 697 real(dp) :: rho(2) 698 699 ! ************************************************************************* 700 701 numneg=0 702 worst=zero 703 704 if(nspden==1)then 705 706 ! Non spin-polarized 707 !$OMP PARALLEL DO PRIVATE(ifft,rhotmp) REDUCTION(MIN:worst) REDUCTION(+:numneg) SHARED(nfft,rhonow) 708 do ifft=1,nfft 709 rhotmp=rhonow(ifft,1) 710 if(rhotmp<xc_denpos)then 711 if(rhotmp<-xc_denpos)then 712 ! This case is probably beyond machine precision considerations 713 worst=min(worst,rhotmp) 714 numneg=numneg+1 715 end if 716 rhonow(ifft,1)=xc_denpos 717 end if 718 end do 719 else if (nspden==2) then 720 721 ! Spin-polarized 722 723 ! rhonow is stored as (up,dn) 724 if (option==0) then 725 726 !$OMP PARALLEL DO PRIVATE(ifft,ispden,rho,rhotmp) REDUCTION(MIN:worst) REDUCTION(+:numneg) & 727 !$OMP&SHARED(nfft,nspden,rhonow) 728 do ifft=1,nfft 729 ! For polarized case, rho(1) is spin-up density, rho(2) is spin-down density 730 rho(1)=rhonow(ifft,1) 731 rho(2)=rhonow(ifft,2) 732 do ispden=1,nspden 733 if (rho(ispden)<xc_denpos) then 734 if (rho(ispden)<-xc_denpos) then 735 ! This case is probably beyond machine precision considerations 736 worst=min(worst,rho(ispden)) 737 numneg=numneg+1 738 end if 739 rhonow(ifft,ispden)=xc_denpos 740 end if 741 end do 742 end do 743 744 ! rhonow is stored as (up+dn,up) 745 else if (option==1) then 746 747 !$OMP PARALLEL DO PRIVATE(ifft,ispden,rho,rhotmp) & 748 !$OMP&REDUCTION(MIN:worst) REDUCTION(+:numneg) & 749 !$OMP&SHARED(nfft,nspden,rhonow) 750 do ifft=1,nfft 751 ! For polarized case, rho(1) is spin-up density, rho(2) is spin-down density 752 rho(1)=rhonow(ifft,2) 753 rho(2)=rhonow(ifft,1)-rho(1) 754 do ispden=1,nspden 755 if (rho(ispden)<xc_denpos) then 756 if (rho(ispden)<-xc_denpos) then 757 ! This case is probably beyond machine precision considerations 758 worst=min(worst,rho(ispden)) 759 numneg=numneg+1 760 end if 761 rho(ispden)=xc_denpos 762 rhonow(ifft,1)=rho(1)+rho(2) 763 rhonow(ifft,2)=rho(1) 764 end if 765 end do 766 end do 767 768 end if ! option 769 770 else 771 ABI_BUG('nspden>2 not allowed !') 772 end if ! End choice between non-spin polarized and spin-polarized. 773 774 if (numneg>0) then 775 if (iwarn==0) then 776 write(message,'(a,i0,a,a,a,es10.2,a,e10.2,11a)')& 777 & 'Density went too small (lower than xc_denpos) at ',numneg,' points',ch10,& 778 & 'and was set to xc_denpos = ',xc_denpos,'. Lowest was ',worst,'.',ch10,& 779 & 'This might be due to (1) too low boxcut or (2) too low ecut for',ch10,& 780 & ' pseudopotential core charge, or (3) too low ecut for estimated initial density.',ch10,& 781 & ' Possible workarounds : increase ecut, or define the input variable densty,',ch10,& 782 & ' with a value larger than the guess for the decay length, or initialize your,',ch10,& 783 & ' density with a preliminary LDA or GGA-PBE if you are using a more exotic xc functional.' 784 ABI_WARNING(message) 785 end if 786 iwarn=iwarn+1 787 end if 788 789 end subroutine mkdenpos
m_drivexc/size_dvxc [ Functions ]
[ Top ] [ m_drivexc ] [ Functions ]
NAME
size_dvxc
FUNCTION
Give the sizes of the several arrays involved in exchange-correlation calculation needed to allocated them for the drivexc routine
INPUTS
ixc= choice of exchange-correlation scheme 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) nspden= number of spin components [xc_funcs(2)]= <type(libxc_functional_type)> [add_tfw]= optional flag controling the addition of Weiszacker gradient correction to Thomas-Fermi XC energy
OUTPUT
--- All optionals [usegradient]= [flag] 1 if the XC functional needs the gradient of the density (grho2_updn) [uselaplacian]= [flag] 1 if the XC functional needs the laplacian of the density (lrho_updn) [usekden]= [flag] 1 if the XC functional needs the kinetic energy density (lrho_updn) [nvxcgrho]= size of the array dvxcdgr(npts,nvxcgrho) (derivative of Exc wrt to gradient) [nvxclrho]= size of the array dvxclpl(npts,nvxclrho) (derivative of Exc wrt to laplacian) [nvxctau]= size of the array dvxctau(npts,nvxctau) (derivative of Exc wrt to kin. ener. density) [ndvxc]= size of the array dvxc(npts,ndvxc) (second derivatives of Exc wrt to density and gradient) [nd2vxc]= size of the array d2vxc(npts,nd2vxc) (third derivatives of Exc wrt density)
SOURCE
423 subroutine size_dvxc(ixc,order,nspden,& 424 & usegradient,uselaplacian,usekden,& 425 & nvxcgrho,nvxclrho,nvxctau,ndvxc,nd2vxc,& 426 & add_tfw,xc_funcs) ! Optional 427 428 !Arguments---------------------- 429 integer,intent(in) :: ixc,nspden,order 430 integer,intent(out),optional :: nvxcgrho,nvxclrho,nvxctau,ndvxc,nd2vxc 431 integer,intent(out),optional :: usegradient,uselaplacian,usekden 432 logical, intent(in),optional :: add_tfw 433 type(libxc_functional_type),intent(in),optional :: xc_funcs(2) 434 435 !Local variables---------------- 436 logical :: libxc_has_kxc,libxc_has_k3xc,libxc_isgga,libxc_ismgga,libxc_ishybrid,my_add_tfw 437 logical :: need_gradient,need_laplacian,need_kden 438 439 ! ************************************************************************* 440 441 !Several flags 442 my_add_tfw=.false.;if (present(add_tfw)) my_add_tfw=add_tfw 443 libxc_isgga=.false. ; libxc_ismgga=.false. ; libxc_ishybrid=.false. 444 if(ixc<0)then 445 if(present(xc_funcs))then 446 libxc_has_kxc=libxc_functionals_has_kxc(xc_funcs) 447 libxc_has_k3xc=libxc_functionals_has_k3xc(xc_funcs) 448 libxc_isgga=libxc_functionals_isgga(xc_functionals=xc_funcs) 449 libxc_ismgga=libxc_functionals_ismgga(xc_functionals=xc_funcs) 450 libxc_ishybrid=libxc_functionals_is_hybrid(xc_functionals=xc_funcs) 451 else 452 libxc_has_kxc=libxc_functionals_has_kxc() 453 libxc_has_k3xc=libxc_functionals_has_k3xc() 454 libxc_isgga=libxc_functionals_isgga() 455 libxc_ismgga=libxc_functionals_ismgga() 456 libxc_ishybrid=libxc_functionals_is_hybrid() 457 end if 458 end if 459 460 !Do we use the gradient? 461 need_gradient=((ixc>=11.and.ixc<=17).or.(ixc==23.or.ixc==24).or. & 462 & (ixc==26.or.ixc==27).or.(ixc>=31.and.ixc<=35).or. & 463 & (ixc==41.or.ixc==42).or.ixc==1402000) 464 if (ixc<0.and.(libxc_isgga.or.libxc_ismgga.or.libxc_ishybrid)) need_gradient=.true. 465 if (my_add_tfw) need_gradient=.true. 466 if (present(usegradient)) usegradient=merge(1,0,need_gradient) 467 468 !Do we use the laplacian? 469 need_laplacian=(ixc==32.or.ixc==35) 470 if (ixc<0) then 471 if(present(xc_funcs)) need_laplacian=libxc_functionals_needs_laplacian(xc_functionals=xc_funcs) 472 if(.not.present(xc_funcs)) need_laplacian=libxc_functionals_needs_laplacian() 473 end if 474 if (present(uselaplacian)) uselaplacian=merge(1,0,need_laplacian) 475 476 !Do we use the kinetic energy density? 477 need_kden=(ixc==31.or.ixc==34.or.ixc==35) 478 if (ixc<0) need_kden=libxc_ismgga 479 if (present(usekden)) usekden=merge(1,0,need_kden) 480 481 !First derivative(s) of XC functional wrt gradient of density 482 if (present(nvxcgrho)) then 483 nvxcgrho=0 484 if (abs(order)>=1) then 485 if (need_gradient) nvxcgrho=3 486 if (ixc==13) nvxcgrho=0 487 if (ixc==16.or.ixc==17.or.ixc==26.or.ixc==27) nvxcgrho=2 488 end if 489 end if 490 491 !First derivative(s) of XC functional wrt laplacian of density 492 if (present(nvxclrho)) then 493 nvxclrho=0 494 if (abs(order)>=1) then 495 if (need_laplacian) nvxclrho=min(nspden,2) 496 end if 497 end if 498 499 !First derivative(s) of XC functional wrt kinetic energy density 500 if (present(nvxctau)) then 501 nvxctau=0 502 if (abs(order)>=1) then 503 if (need_kden) nvxctau=min(nspden,2) 504 end if 505 end if 506 507 !Second derivative(s) of XC functional wrt density 508 if (present(ndvxc)) then 509 ndvxc=0 510 if (abs(order)>=2) then 511 if (ixc==1.or.ixc==7.or.ixc==8.or.ixc==9.or.ixc==10.or.ixc==13.or. & 512 & ixc==21.or.ixc==22) then 513 ndvxc=min(nspden,2)+1 514 else if ((ixc>=2.and.ixc<=6).or.(ixc>=31.and.ixc<=35).or.ixc==50) then 515 ndvxc=1 516 else if (ixc==12.or.ixc==24) then 517 ndvxc=8 518 else if (ixc==11.or.ixc==12.or.ixc==14.or.ixc==15.or. & 519 & ixc==23.or.ixc==41.or.ixc==42.or.ixc==1402000) then 520 ndvxc=15 521 else if (ixc<0) then 522 if (libxc_has_kxc.or.ixc==-406.or.ixc==-427.or.ixc==-428.or.ixc==-456) then 523 ndvxc=2*min(nspden,2)+1 ; if (order==-2) ndvxc=2 524 if (need_gradient) ndvxc=15 ! This is for GGA, but also for mGGA 525 ! (we dont consider derivatives wrt Tau or Laplacian) 526 end if 527 end if 528 end if 529 end if 530 531 !Third derivative(s) of XC functional wrt density 532 if (present(nd2vxc)) then 533 nd2vxc=0 534 if (abs(order)>=3) then 535 if (ixc==3.or.(ixc>=11.and.ixc<=15.and.ixc/=13).or. & 536 & ixc==23.or.ixc==24.or.ixc==41.or.ixc==42) then 537 nd2vxc=1 538 else if ((ixc>=7.and.ixc<=10).or.ixc==13.or.ixc==1402000) then 539 nd2vxc=3*min(nspden,2)-2 540 else if (ixc<0) then 541 if (libxc_has_k3xc) then 542 if (.not.need_gradient) nd2vxc=3*min(nspden,2)-2 543 end if 544 end if 545 end if 546 end if 547 548 end subroutine size_dvxc
m_drivexc/xcmult [ Functions ]
[ Top ] [ m_drivexc ] [ Functions ]
NAME
xcmult
FUNCTION
In the case of GGA, multiply the different gradient of spin-density by the derivative of the XC functional with respect to the norm of the gradient, then divide it by the norm of the gradient
INPUTS
depsxc(nfft,nspgrad)=derivative of Exc with respect to the (spin-)density,
or to the norm of the gradient of the (spin-)density,
further divided by the norm of the gradient of the (spin-)density
The different components of depsxc will be
for nspden=1, depsxc(:,1)=d(rho.exc)/d(rho)
and if ngrad=2, depsxc(:,2)=1/2*1/|grad rho_up|*d(rho.exc)/d(|grad rho_up|)
+ 1/|grad rho|*d(rho.exc)/d(|grad rho|)
(do not forget : |grad rho| /= |grad rho_up| + |grad rho_down|
for nspden=2, depsxc(:,1)=d(rho.exc)/d(rho_up)
depsxc(:,2)=d(rho.exc)/d(rho_down)
and if ngrad=2, depsxc(:,3)=1/|grad rho_up|*d(rho.exc)/d(|grad rho_up|)
depsxc(:,4)=1/|grad rho_down|*d(rho.exc)/d(|grad rho_down|)
depsxc(:,5)=1/|grad rho|*d(rho.exc)/d(|grad rho|)
nfft=(effective) number of FFT grid points (for this processor)
ngrad = must be 2
nspden=number of spin-density components
nspgrad=number of spin-density and spin-density-gradient components
OUTPUT
(see side effects)
SIDE EFFECTS
rhonow(nfft,nspden,ngrad*ngrad)=
at input :
electron (spin)-density in real space and its gradient,
either on the unshifted grid (if ishift==0,
then equal to rhor), or on the shifted grid
rhonow(:,:,1)=electron density in electrons/bohr**3
rhonow(:,:,2:4)=gradient of electron density in el./bohr**4
at output :
rhonow(:,:,2:4) has been multiplied by the proper factor,
described above.
SOURCE
596 subroutine xcmult (depsxc,nfft,ngrad,nspden,nspgrad,rhonow) 597 598 !Arguments ------------------------------------ 599 !scalars 600 integer,intent(in) :: nfft,ngrad,nspden,nspgrad 601 !arrays 602 real(dp),intent(in) :: depsxc(nfft,nspgrad) 603 real(dp),intent(inout) :: rhonow(nfft,nspden,ngrad*ngrad) 604 605 !Local variables------------------------------- 606 !scalars 607 integer :: idir,ifft 608 real(dp) :: rho_tot,rho_up 609 610 ! ************************************************************************* 611 612 do idir=1,3 613 614 if(nspden==1)then 615 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(depsxc,idir,nfft,rhonow) 616 do ifft=1,nfft 617 rhonow(ifft,1,1+idir)=rhonow(ifft,1,1+idir)*depsxc(ifft,2) 618 end do 619 620 else 621 622 ! In the spin-polarized case, there are more factors to take into account 623 !$OMP PARALLEL DO PRIVATE(ifft,rho_tot,rho_up) SHARED(depsxc,idir,nfft,rhonow) 624 do ifft=1,nfft 625 rho_tot=rhonow(ifft,1,1+idir) 626 rho_up =rhonow(ifft,2,1+idir) 627 rhonow(ifft,1,1+idir)=rho_up *depsxc(ifft,3) + rho_tot*depsxc(ifft,5) 628 rhonow(ifft,2,1+idir)=(rho_tot-rho_up)*depsxc(ifft,4)+ rho_tot*depsxc(ifft,5) 629 end do 630 631 end if ! nspden==1 632 633 end do ! End loop on directions 634 635 end subroutine xcmult