TABLE OF CONTENTS
- ABINIT/check_kxc
- ABINIT/drivexc
- ABINIT/drivexc_main
- ABINIT/echo_xc_name
- ABINIT/m_drivexc
- ABINIT/mkdenpos
- ABINIT/size_dvxc
- ABINIT/xcmult
ABINIT/check_kxc [ Functions ]
NAME
check_kxc
FUNCTION
Given a XC functional (defined by ixc), check if Kxc (dVxc/drho) is avalaible.
INPUTS
ixc = internal code for xc functional optdriver=type of calculation (ground-state, response function, GW, ...)
OUTPUT
PARENTS
respfn,scfcv
CHILDREN
SOURCE
610 subroutine check_kxc(ixc,optdriver) 611 612 613 !This section has been created automatically by the script Abilint (TD). 614 !Do not modify the following lines by hand. 615 #undef ABI_FUNC 616 #define ABI_FUNC 'check_kxc' 617 !End of the abilint section 618 619 implicit none 620 621 !Arguments ------------------------------- 622 integer, intent(in) :: ixc,optdriver 623 624 !Local variables ------------------------- 625 logical :: kxc_available 626 character(len=500) :: msg 627 628 ! ********************************************************************* 629 630 kxc_available=.false. 631 632 if (ixc>=0) then 633 kxc_available=(ixc/=16.and.ixc/=17.and.ixc/=26.and.ixc/=27) 634 if (.not.kxc_available) then 635 write(msg,'(a,i0,3a)') & 636 & 'The selected XC functional (ixc=',ixc,')',ch10,& 637 & 'does not provide Kxc (dVxc/drho) !' 638 end if 639 else if (ixc==-406.or.ixc==-427.or.ixc==-428.or.ixc==-456)then 640 kxc_available=.true. 641 else ! ixc<0 and not one of the allowed hybrids 642 kxc_available=libxc_functionals_has_kxc() 643 if (.not.kxc_available) then 644 write(msg,'(a,i0,7a)') & 645 & 'The selected XC functional (ixc=',ixc,'):',ch10,& 646 & ' <<',trim(libxc_functionals_fullname()),'>>',ch10,& 647 & 'does not provide Kxc (dVxc/drho) !' 648 end if 649 end if 650 651 if (.not.kxc_available) then 652 write(msg,'(7a)') trim(msg),ch10,& 653 & 'However, with the current input options, ABINIT needs Kxc.',ch10,& 654 & '>Possible action:',ch10,& 655 & 'Change the XC functional in psp file or input file.' 656 if (optdriver==0) then 657 write(msg,'(13a)') trim(msg),ch10,& 658 & '>Possible action (2):',ch10,& 659 & 'If you are using density mixing for the SCF cycle',ch10,& 660 & '(iscf>=10, which is the default for PAW),',ch10,& 661 & 'change to potential mixing (iscf=7, for instance).',ch10,& 662 & '>Possible action (3):',ch10,& 663 & 'Switch to another value of densfor_pred (=5, for instance).' 664 end if 665 MSG_ERROR(msg) 666 end if 667 668 end subroutine check_kxc
ABINIT/drivexc [ Functions ]
NAME
drivexc
FUNCTION
Driver of XC functionals. Treat spin-polarized as well as non-spin-polarized. Treat local approximations or GGAs. Optionally, deliver the XC kernel, or even the derivative of the XC kernel (the third derivative of the XC energy)
INPUTS
ixc=number of the XC functional ndvxc= size of dvxc(npts,ndvxc) ngr2= size of grho2_updn(npts,ngr2) nvxcgrho= size of vxcgrho(npts,nvxcgrho) npts=number of real space points on which the density (and its gradients, if needed) is provided nspden=number of spin-density components (1 or 2) order=gives the maximal derivative of Exc computed. 1=usual value (return exc and vxc) 2=also computes the kernel (return exc,vxc,kxc) -2=like 2, except (to be described) 3=also computes the derivative of the kernel (return exc,vxc,kxc,k3xc) rho_updn(npts,nspden)=the spin-up and spin-down densities If nspden=1, only the spin-up density must be given. In the calling routine, the spin-down density must be equal to the spin-up density, and both are half the total density. If nspden=2, the spin-up and spin-down density must be given Optional inputs: [el_temp]= electronic temperature (to be used for finite temperature XC functionals) [exexch]= choice of <<<local>>> exact exchange. Active if exexch=3 (only for GGA, and NOT for libxc) [grho2_updn(npts,ngr2)]=the square of the gradients of spin-up, spin-down, and total density If nspden=1, only the square of the gradient of the spin-up density must be given. In the calling routine, the square of the gradient of the spin-down density must be equal to the square of the gradient of the spin-up density, and both must be equal to one-quarter of the square of the gradient of the total density. If nspden=2, the square of the gradients of spin-up, spin-down, and total density must be given. Note that the square of the gradient of the total density is usually NOT related to the square of the gradient of the spin-up and spin-down densities, because the gradients are not usually aligned. This is not the case when nspden=1. [hyb_mixing]= mixing parameter for the native PBEx functionals (ixc=41 and 42) [lrho_updn(npts,nspden)]=the Laplacian of spin-up and spin-down densities If nspden=1, only the spin-up Laplacian density must be given. In the calling routine, the spin-down Laplacian density must be equal to the spin-up Laplacian density, and both are half the total Laplacian density. If nspden=2, the Laplacian of spin-up and spin-down densities must be given [tau_updn(npts,nspden)]=the spin-up and spin-down kinetic energy densities If nspden=1, only the spin-up kinetic energy density must be given. In the calling routine, the spin-down kinetic energy density must be equal to the spin-up kinetic energy density, and both are half the total kinetic energy density. If nspden=2, the spin-up and spin-down kinetic energy densities must be given [xc_funcs(2)]= <type(libxc_functional_type)>, optional : libxc XC functionals. If not specified, the underlying xc_global(2) is used by libxc. [xc_tb09_c]= c parameter for the Tran-Blaha mGGA functional, within libxc
OUTPUT
exc(npts)=exchange-correlation energy density (hartree) vxcrho(npts,nspden)= (d($\rho$*exc)/d($\rho_up$)) (hartree) and (d($\rho$*exc)/d($\rho_down$)) (hartree) vxcgrho(npts,3)= 1/$|grad \rho_up|$ (d($\rho$*exc)/d($|grad \rho_up|$)) (hartree) 1/$|grad \rho_dn|$ (d($\rho$*exc)/d($|grad \rho_dn|$)) (hartree) and 1/$|grad \rho|$ (d($\rho$*exc)/d($|grad \rho|$)) (hartree) (will be zero if a LDA functional is used) vxclrho(npts,nspden)=(only for meta-GGA, i.e. optional output)= (d($\rho$*exc)/d($\lrho_up$)) (hartree) and (d($\rho$*exc)/d($\lrho_down$)) (hartree) vxctau(npts,nspden)=(only for meta-GGA, i.e. optional output)= derivative of XC energy density with respect to kinetic energy density (depsxcdtau). (d($\rho$*exc)/d($\tau_up$)) (hartree) and (d($\rho$*exc)/d($\tau_down$)) (hartree) Optional output: if(abs(order)>1) [dvxc(npts,ndvxc)]=partial second derivatives of the xc energy (This is a mess, to be rationalized !!) In case of local energy functional (option=1,-1 or 3): dvxc(npts,1+nspden)= (Hartree*bohr^3) if(nspden=1 .and. order==2): dvxci(:,1)=dvxc/d$\rho$ , dvxc(:,2) empty if(nspden=1 .and. order==-2): also compute dvxci(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$ if(nspden=2): dvxc(:,1)=dvxc($\uparrow$)/d$\rho(\downarrow)$, dvxc(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$, dvxc(:,3)=dvxc($\downarrow$)/d$\rho(\downarrow)$ In case of gradient corrected functional (option=2,-2, 4, -4, 5, 6, 7): dvxc(npts,15)= dvxc(:,1)= d2Ex/drho_up drho_up dvxc(:,2)= d2Ex/drho_dn drho_dn dvxc(:,3)= dEx/d(abs(grad(rho_up))) / abs(grad(rho_up)) dvxc(:,4)= dEx/d(abs(grad(rho_dn))) / abs(grad(rho_dn)) dvxc(:,5)= d2Ex/d(abs(grad(rho_up))) drho_up / abs(grad(rho_up)) dvxc(:,6)= d2Ex/d(abs(grad(rho_dn))) drho_dn / abs(grad(rho_dn)) dvxc(:,7)= 1/abs(grad(rho_up)) * d/d(abs(grad(rho_up)) (dEx/d(abs(grad(rho_up))) /abs(grad(rho_up))) dvxc(:,8)= 1/abs(grad(rho_dn)) * d/d(abs(grad(rho_dn)) (dEx/d(abs(grad(rho_dn))) /abs(grad(rho_dn))) dvxc(:,9)= d2Ec/drho_up drho_up dvxc(:,10)=d2Ec/drho_up drho_dn dvxc(:,11)=d2Ec/drho_dn drho_dn dvxc(:,12)=dEc/d(abs(grad(rho))) / abs(grad(rho)) dvxc(:,13)=d2Ec/d(abs(grad(rho))) drho_up / abs(grad(rho)) dvxc(:,14)=d2Ec/d(abs(grad(rho))) drho_dn / abs(grad(rho)) dvxc(:,15)=1/abs(grad(rho)) * d/d(abs(grad(rho)) (dEc/d(abs(grad(rho))) /abs(grad(rho))) if(abs(order)>2) (only available for LDA and nspden=1) [d2vxc(npts,nd2vxc)]=partial third derivatives of the xc energy if nspden=1 d2vxc(npts,1)=second derivative of the XC potential=3rd order derivative of energy if nspden=2 d2vxc(npts,1), d2vxc(npts,2), d2vxc(npts,3), d2vxc(npts,4) (3rd derivative of energy) [fxcT(npts)]=XC free energy of the electron gaz at finite temperature (to be used for plasma systems)
PARENTS
drivexc_main
CHILDREN
invcb,libxc_functionals_end,libxc_functionals_getvxc libxc_functionals_init,xchcth,xchelu,xciit,xclb,xcpbe,xcpzca,xcspol xctetr,xcwign,xcxalp
SOURCE
1208 subroutine drivexc(exc,ixc,npts,nspden,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 1209 & dvxc,d2vxc,grho2_updn,vxcgrho,el_temp,exexch,fxcT,& !Optional arguments 1210 & hyb_mixing,lrho_updn,vxclrho,tau_updn,vxctau,xc_funcs,xc_tb09_c) !Optional arguments 1211 1212 1213 !This section has been created automatically by the script Abilint (TD). 1214 !Do not modify the following lines by hand. 1215 #undef ABI_FUNC 1216 #define ABI_FUNC 'drivexc' 1217 !End of the abilint section 1218 1219 implicit none 1220 1221 !Arguments ------------------------------------ 1222 !scalars 1223 integer,intent(in) :: ixc,ndvxc,ngr2,nd2vxc,npts,nspden,nvxcgrho,order 1224 integer,intent(in),optional :: exexch 1225 real(dp),intent(in),optional :: el_temp,hyb_mixing,xc_tb09_c 1226 !arrays 1227 real(dp),intent(in) :: rho_updn(npts,nspden) 1228 real(dp),intent(in),optional :: grho2_updn(npts,ngr2) 1229 real(dp),intent(in),optional :: lrho_updn(npts,nspden), tau_updn(npts,nspden) 1230 real(dp),intent(out) :: exc(npts),vxcrho(npts,nspden) 1231 real(dp),intent(out),optional :: d2vxc(npts,nd2vxc),dvxc(npts,ndvxc),fxcT(:) 1232 real(dp),intent(out),optional :: vxcgrho(npts,nvxcgrho) 1233 real(dp),intent(out),optional :: vxclrho(npts,nspden),vxctau(npts,nspden) 1234 type(libxc_functional_type),intent(inout),optional :: xc_funcs(2) 1235 1236 !Local variables------------------------------- 1237 !scalars 1238 integer :: exexch_,ixc_from_lib,ixc1,ixc2,ndvxc_x,optpbe,ispden 1239 logical :: libxc_test,xc_err_ndvxc,xc_err_nvxcgrho1,xc_err_nvxcgrho2 1240 logical :: is_gga,is_mgga 1241 real(dp) :: alpha 1242 real(dp),parameter :: rsfac=0.6203504908994000e0_dp 1243 character(len=500) :: message 1244 !arrays 1245 real(dp),allocatable :: exci_rpa(:) 1246 real(dp),allocatable :: rhotot(:),rspts(:),vxci_rpa(:,:),zeta(:) 1247 real(dp),allocatable :: exc_c(:),exc_x(:),vxcrho_c(:,:),vxcrho_x(:,:) 1248 real(dp),allocatable :: d2vxc_c(:,:),d2vxc_x(:,:),dvxc_c(:,:),dvxc_x(:,:) 1249 real(dp),allocatable :: vxcgrho_x(:,:) 1250 type(libxc_functional_type) :: xc_funcs_vwn3(2),xc_funcs_lyp(2) 1251 1252 ! ************************************************************************* 1253 1254 ! ================================================= 1255 ! == Compatibility tests == 1256 ! ================================================= 1257 1258 !Checks the values of order and other size parameters 1259 if( (order<1 .and. order/=-2) .or. order>4)then 1260 write(message, '(a,i0)' )& 1261 & 'The only allowed values for order are 1,2,-2 or 3, while it is found to be ',order 1262 MSG_BUG(message) 1263 end if 1264 xc_err_ndvxc=.false. 1265 xc_err_nvxcgrho1=.false. 1266 xc_err_nvxcgrho2=.false. 1267 if (ixc==16.or.ixc==17.or.ixc==26.or.ixc==27) xc_err_nvxcgrho1=(nvxcgrho/=2) 1268 if (order**2>1) then 1269 if ((ixc>=2.and.ixc<=6).or.(ixc==50)) xc_err_ndvxc=(ndvxc/=1) 1270 if (ixc==1.or.ixc==21.or.ixc==22) xc_err_ndvxc=(ndvxc/=nspden+1) 1271 end if 1272 if (order==2) then 1273 if ((ixc>=7.and.ixc<=10).or.(ixc==13)) xc_err_nvxcgrho1=(ndvxc/=1+nspden.or.nvxcgrho/=0) 1274 if (ixc==12.or.ixc==24) xc_err_nvxcgrho1=(ndvxc/=8.or.nvxcgrho/=3) 1275 if (ixc==11.or.ixc==14.or.ixc==15.or.ixc==23.or.ixc==41.or.ixc==42) xc_err_nvxcgrho1=(ndvxc/=15.or.nvxcgrho/=3) 1276 end if 1277 if (order==3) then 1278 if (ixc==3) xc_err_ndvxc=(ndvxc/=1) 1279 if (ixc==10.or.ixc==13) xc_err_nvxcgrho1=(ndvxc/=1+nspden.or.nvxcgrho/=0) 1280 if (ixc==12.or.ixc==24) xc_err_nvxcgrho1=(ndvxc/=8.or.nvxcgrho/=3) 1281 if (ixc==11.or.ixc==14.or.ixc==15.or.ixc==23.or.ixc==41.or.ixc==42) xc_err_nvxcgrho1=(ndvxc/=15.or.nvxcgrho/=3) 1282 if (ixc>=7.and.ixc<=9) xc_err_nvxcgrho2=(ndvxc/=1+nspden.or.nvxcgrho/=0.or.nd2vxc/=(3*nspden-2)) 1283 if (ixc==50) xc_err_ndvxc=(nd2vxc/=0) 1284 end if 1285 if (xc_err_ndvxc) then 1286 write(message, '(7a,i0,a,i0,a)' )& 1287 & 'Wrong value of ndvxc:',ch10,& 1288 & 'the value of the spin polarization',ch10,& 1289 & 'is not compatible with the value of ixc:',ch10,& 1290 & 'ixc=',ixc,' (nspden=',nspden,')' 1291 MSG_BUG(message) 1292 end if 1293 if (xc_err_nvxcgrho1) then 1294 write(message,'(3a,i0,a,i0,a,i0)' )& 1295 & 'Wrong value of ndvxc or nvxcgrho:',ch10,& 1296 & 'ixc=',ixc,'ndvxc=',ndvxc,'nvxcgrho=',nvxcgrho 1297 MSG_BUG(message) 1298 end if 1299 if (xc_err_nvxcgrho2) then 1300 write(message, '(3a,i0,a,i0,a,i0,a,i0)' )& 1301 & 'Wrong value of ndvxc or nvxcgrho or nd2vxc:',ch10,& 1302 & 'ixc=',ixc,'ndvxc=',ndvxc,'nvxcgrho=',nvxcgrho,'nd2vxc=',nd2vxc 1303 MSG_BUG(message) 1304 end if 1305 1306 !Check libXC 1307 if (ixc<0 .or. ixc==1402) then 1308 libxc_test=libxc_functionals_check(stop_if_error=.true.) 1309 end if 1310 1311 if (ixc<0) then 1312 ! Prepare the tests 1313 if(present(xc_funcs))then 1314 is_gga=libxc_functionals_isgga(xc_functionals=xc_funcs) 1315 is_mgga=libxc_functionals_ismgga(xc_functionals=xc_funcs) 1316 ixc_from_lib=libxc_functionals_ixc(xc_functionals=xc_funcs) 1317 else 1318 is_gga=libxc_functionals_isgga() 1319 is_mgga=libxc_functionals_ismgga() 1320 ixc_from_lib=libxc_functionals_ixc() 1321 end if 1322 ! Check consistency between ixc passed in input and the one used to initialize the library. 1323 if (ixc /= ixc_from_lib) then 1324 write(message, '(a,i0,2a,i0,2a)')& 1325 & 'The value of ixc specified in input, ixc = ',ixc,ch10,& 1326 & 'differs from the one used to initialize the functional ',ixc_from_lib,ch10,& 1327 & 'Action: reinitialize the global structure funcs, see NOTES in m_libxc_functionals' 1328 MSG_BUG(message) 1329 end if 1330 else if (ixc==1420)then 1331 is_gga=.true. 1332 is_mgga=.false. 1333 end if 1334 1335 if (ixc<0 .or. ixc==1402) then 1336 ! Check whether all the necessary arrays are present and have the correct dimensions 1337 if (is_gga .or. is_mgga) then 1338 if ( (.not. present(grho2_updn)) .or. (.not. present(vxcgrho))) then 1339 write(message, '(5a,2L2,2a,2L2,2a,i10,a,i5,a,i5)' )& 1340 & 'At least one of the functionals is a GGA or a MGGA,',ch10,& 1341 & 'but not all the necessary arrays are present.',ch10,& 1342 & 'is_gga, is_mgga=',is_gga,is_mgga,ch10,& 1343 & 'present(grho2_updn),present(vxcgrho)=',present(grho2_updn),present(vxcgrho),ch10,& 1344 & 'ixc=',ixc,' nvxcgrho=',nvxcgrho,' ngr2=',ngr2 1345 MSG_BUG(message) 1346 end if 1347 if (ngr2==0.or.nvxcgrho/=3) then 1348 write(message, '(5a,i7,a,i6,a,i6)' )& 1349 & 'The value of the number of the XC functional ixc',ch10,& 1350 & 'is not compatible with the value of nvxcgrho or ngr2',ch10,& 1351 & 'ixc=',ixc,' nvxcgrho=',nvxcgrho,' ngr2=',ngr2 1352 MSG_BUG(message) 1353 end if 1354 end if 1355 if (is_mgga) then 1356 if ( (.not. present(lrho_updn)) .or. (.not. present(vxclrho)) .or. & 1357 (.not. present(tau_updn)) .or. (.not. present(vxctau)) ) then 1358 write(message, '(5a,i7)' )& 1359 & 'At least one of the functionals is a MGGA,',ch10,& 1360 & 'but not all the necessary arrays are present.',ch10,& 1361 & 'ixc=',ixc 1362 MSG_BUG(message) 1363 end if 1364 end if 1365 end if 1366 1367 !Checks the compatibility between the inputs and the presence of the optional arguments 1368 if(present(dvxc))then 1369 if(order**2<=1.or.ixc==16.or.ixc==17.or.ixc==26.or.ixc==27)then 1370 write(message, '(5a,i0,a,i0)' )& 1371 & 'The value of the number of the XC functional ixc',ch10,& 1372 & 'or the value of order is not compatible with the presence of the array dvxc',ch10,& 1373 & 'ixc=',ixc,'order=',order 1374 MSG_BUG(message) 1375 end if 1376 end if 1377 if(present(d2vxc))then 1378 if(order/=3.or.(ixc/=3.and.(((ixc>15).and.(ixc/=23)) .or.& 1379 & (ixc>=0.and.ixc<7).or.(ixc>=40 .and. ixc/=1402))) )then 1380 write(message, '(5a,i6,a,i0)' )& 1381 & 'The value of the number of the XC functional ixc',ch10,& 1382 & 'or the value of order is not compatible with the presence of the array d2vxc',ch10,& 1383 & 'ixc=',ixc,'order=',order 1384 MSG_BUG(message) 1385 end if 1386 end if 1387 if(present(vxcgrho))then 1388 if(nvxcgrho==0.or. & 1389 & ((((ixc>17.and.ixc/=23.and.ixc/=24.and.ixc/=26.and.ixc/=27) .or.& 1390 & (ixc>= 0.and.ixc<7).or.(ixc>=40 .and. ixc/=1402)) .and. nvxcgrho /=3 ))) then 1391 if(ixc<31.and.ixc>34)then !! additional if to include ixc 31 to 34 in the list (these ixc are used for mgga test see below) 1392 write(message, '(5a,i0,a,i0)' )& 1393 & 'The value of the number of the XC functional ixc',ch10,& 1394 & 'or the value of nvxcgrho is not compatible with the presence of the array vxcgrho',ch10,& 1395 & 'ixc=',ixc,' nvxcgrho=',nvxcgrho 1396 MSG_BUG(message) 1397 end if 1398 end if 1399 end if 1400 if(present(grho2_updn))then 1401 if (ngr2/=2*nspden-1 ) then 1402 write(message, '(4a)' ) ch10,& 1403 & 'drivexc : BUG -',ch10,& 1404 & 'ngr2 must be 2*nspden-1 !' 1405 ! MSG_BUG(message) 1406 end if 1407 if ((ixc>0.and.ixc<11).or.(ixc>17.and.ixc<23).or.ixc==25.or.(ixc>27.and.ixc<31).or. & 1408 & (ixc>34.and.ixc<41).or.ixc==50) then 1409 write(message, '(5a,i0)' )& 1410 & 'The value of the number of the XC functional ixc',ch10,& 1411 & 'is not compatible with the presence of the array grho2_updn',ch10,& 1412 & 'ixc=',ixc 1413 MSG_BUG(message) 1414 end if 1415 end if 1416 if (present(exexch)) then 1417 if(exexch/=0.and.(.not.present(grho2_updn))) then 1418 message='exexch argument only valid for GGA!' 1419 MSG_BUG(message) 1420 end if 1421 end if 1422 if(ixc==31.or.ixc==32.or.ixc==33.or.ixc==34) then 1423 if((.not.(present(vxcgrho))) .or. (.not.(present(vxclrho))) .or. (.not.(present(vxctau))))then 1424 message = 'vxcgrho or vxclrho or vxctau is not present but they are all needed for MGGA XC tests.' 1425 MSG_BUG(message) 1426 end if 1427 end if 1428 if(ixc==50) then 1429 if(.not.(present(el_temp)).or.(.not.present(fxcT)))then 1430 message = 'el_temp or fxcT is not present but are needed for IIT XC functional.' 1431 MSG_BUG(message) 1432 end if 1433 if (size(fxcT)/=npts) then 1434 MSG_BUG('fxcT size must be npts!') 1435 end if 1436 end if 1437 1438 ! ================================================= 1439 ! == Intermediate quantities computation == 1440 ! ================================================= 1441 1442 !If needed, compute rhotot and rs 1443 if (ixc==1.or.ixc==2.or.ixc==3.or.ixc==4.or.ixc==5.or.ixc==6.or.ixc==21.or.ixc==22.or.ixc==50) then 1444 ABI_ALLOCATE(rhotot,(npts)) 1445 ABI_ALLOCATE(rspts,(npts)) 1446 if(nspden==1)then 1447 rhotot(:)=two*rho_updn(:,1) 1448 else 1449 rhotot(:)=rho_updn(:,1)+rho_updn(:,2) 1450 end if 1451 call invcb(rhotot,rspts,npts) 1452 rspts(:)=rsfac*rspts(:) 1453 end if 1454 1455 !If needed, compute zeta 1456 if (ixc==1.or.ixc==21.or.ixc==22) then 1457 ABI_ALLOCATE(zeta,(npts)) 1458 if(nspden==1)then 1459 zeta(:)=zero 1460 else 1461 zeta(:)=two*rho_updn(:,1)/rhotot(:)-one 1462 end if 1463 end if 1464 1465 ! ================================================= 1466 ! == XC energy, potentiel, ... computation == 1467 ! ================================================= 1468 1469 exexch_=0;if(present(exexch)) exexch_=exexch 1470 1471 !>>>>> No exchange-correlation 1472 if (ixc==0.or.ixc==40) then 1473 exc=zero ; vxcrho=zero 1474 if(present(d2vxc)) d2vxc(:,:)=zero 1475 if(present(dvxc)) dvxc(:,:)=zero 1476 if(present(vxcgrho)) vxcgrho(:,:)=zero 1477 if(present(vxclrho)) vxclrho(:,:)=zero 1478 if(present(vxctau)) vxctau(:,:)=zero 1479 1480 !>>>>> New Teter fit (4/93) to Ceperley-Alder data, with spin-pol option 1481 else if (ixc==1 .or. ixc==21 .or. ixc==22) then 1482 ! new Teter fit (4/93) to Ceperley-Alder data, with spin-pol option 1483 if (order**2 <= 1) then 1484 call xcspol(exc,npts,nspden,order,rspts,vxcrho,zeta,ndvxc) 1485 else 1486 call xcspol(exc,npts,nspden,order,rspts,vxcrho,zeta,ndvxc,dvxc) 1487 end if 1488 1489 !>>>>> Perdew-Zunger fit to Ceperly-Alder data (no spin-pol) 1490 else if (ixc==2) then 1491 if (order**2 <= 1) then 1492 call xcpzca(exc,npts,order,rhotot,rspts,vxcrho(:,1)) 1493 else 1494 call xcpzca(exc,npts,order,rhotot,rspts,vxcrho(:,1),dvxc) 1495 end if 1496 1497 !>>>>> Teter fit (4/91) to Ceperley-Alder values (no spin-pol) 1498 else if (ixc==3) then 1499 if (order**2 <= 1) then 1500 call xctetr(exc,npts,order,rhotot,rspts,vxcrho(:,1)) 1501 else if (order == 2) then 1502 call xctetr(exc,npts,order,rhotot,rspts,vxcrho(:,1),dvxc=dvxc) 1503 else if (order == 3) then 1504 call xctetr(exc,npts,order,rhotot,rspts,vxcrho(:,1),d2vxc=d2vxc,dvxc=dvxc) 1505 end if 1506 1507 !>>>>> Wigner xc (no spin-pol) 1508 else if (ixc==4) then 1509 if (order**2 <= 1) then 1510 call xcwign(exc,npts,order,rspts,vxcrho(:,1)) 1511 else 1512 call xcwign(exc,npts,order,rspts,vxcrho(:,1),dvxc) 1513 end if 1514 1515 !>>>>> Hedin-Lundqvist xc (no spin-pol) 1516 else if (ixc==5) then 1517 if (order**2 <= 1) then 1518 call xchelu(exc,npts,order,rspts,vxcrho(:,1)) 1519 else 1520 call xchelu(exc,npts,order,rspts,vxcrho(:,1),dvxc) 1521 end if 1522 1523 !>>>>> X-alpha (no spin-pol) 1524 else if (ixc==6) then 1525 if (order**2 <= 1) then 1526 call xcxalp(exc,npts,order,rspts,vxcrho(:,1)) 1527 else 1528 call xcxalp(exc,npts,order,rspts,vxcrho(:,1),dvxc) 1529 end if 1530 1531 !>>>>> PBE and alternatives 1532 else if (((ixc>=7.and.ixc<=15).or.(ixc>=23.and.ixc<=24)).and.ixc/=10.and.ixc/=13) then 1533 ! Perdew-Wang LSD is coded in Perdew-Burke-Ernzerhof GGA, with optpbe=1 1534 if(ixc==7)optpbe=1 1535 ! x-only part of Perdew-Wang 1536 if(ixc==8)optpbe=-1 1537 ! Exchange + RPA correlation from Perdew-Wang 1538 if(ixc==9)optpbe=3 1539 ! Perdew-Burke-Ernzerhof GGA 1540 if(ixc==11)optpbe=2 1541 ! x-only part of PBE 1542 if(ixc==12)optpbe=-2 1543 ! C09x exchange of V. R. Cooper 1544 if(ixc==24)optpbe=-4 1545 ! revPBE of Zhang and Yang 1546 if(ixc==14)optpbe=5 1547 ! RPBE of Hammer, Hansen and Norskov 1548 if(ixc==15)optpbe=6 1549 ! Wu and Cohen 1550 if(ixc==23)optpbe=7 1551 if (ixc >=7.and.ixc<=9) then 1552 if (order**2 <= 1) then 1553 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 1554 else if (order /=3) then 1555 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,dvxci=dvxc) 1556 else if (order ==3) then 1557 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,d2vxci=d2vxc,dvxci=dvxc) 1558 end if 1559 else if ((ixc >= 11 .and. ixc <= 15) .or. (ixc>=23 .and. ixc<=24)) then 1560 if (order**2 <= 1) then 1561 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 1562 & dvxcdgr=vxcgrho,exexch=exexch_,grho2_updn=grho2_updn) 1563 else if (order /=3) then 1564 if(ixc==12 .or. ixc==24)then 1565 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 1566 & dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 1567 else if(ixc/=12 .or. ixc/=24) then 1568 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 1569 & dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 1570 end if 1571 else if (order ==3) then 1572 if(ixc==12 .or. ixc==24)then 1573 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 1574 & d2vxci=d2vxc,dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 1575 else if(ixc/=12 .or. ixc/=24) then 1576 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 1577 & d2vxci=d2vxc,dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 1578 end if 1579 end if 1580 end if 1581 1582 !>>>>> RPA correlation from Perdew-Wang 1583 else if (ixc==10) then 1584 if (order**2 <= 1) then 1585 ABI_ALLOCATE(exci_rpa,(npts)) 1586 ABI_ALLOCATE(vxci_rpa,(npts,2)) 1587 optpbe=3 1588 call xcpbe(exci_rpa,npts,nspden,optpbe,order,rho_updn,vxci_rpa,ndvxc,ngr2,nd2vxc) 1589 optpbe=1 1590 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 1591 exc(:)=exc(:)-exci_rpa(:) 1592 ! PMA: second index of vxcrho is nspden while that of rpa is 2 they can mismatch 1593 vxcrho(:,1:min(nspden,2))=vxcrho(:,1:min(nspden,2))-vxci_rpa(:,1:min(nspden,2)) 1594 ABI_DEALLOCATE(exci_rpa) 1595 ABI_DEALLOCATE(vxci_rpa) 1596 else if (order /=3) then 1597 ABI_ALLOCATE(exci_rpa,(npts)) 1598 ABI_ALLOCATE(vxci_rpa,(npts,2)) 1599 optpbe=3 1600 call xcpbe(exci_rpa,npts,nspden,optpbe,order,rho_updn,vxci_rpa,ndvxc,ngr2,nd2vxc,dvxci=dvxc) 1601 optpbe=1 1602 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,dvxci=dvxc) 1603 exc(:)=exc(:)-exci_rpa(:) 1604 vxcrho(:,:)=vxcrho(:,:)-vxci_rpa(:,:) 1605 ABI_DEALLOCATE(exci_rpa) 1606 ABI_DEALLOCATE(vxci_rpa) 1607 else if (order ==3) then 1608 ABI_ALLOCATE(exci_rpa,(npts)) 1609 ABI_ALLOCATE(vxci_rpa,(npts,2)) 1610 optpbe=3 1611 call xcpbe(exci_rpa,npts,nspden,optpbe,order,rho_updn,vxci_rpa,ndvxc,ngr2,nd2vxc,& 1612 & d2vxci=d2vxc,dvxci=dvxc) 1613 optpbe=1 1614 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 1615 & d2vxci=d2vxc,dvxci=dvxc) 1616 exc(:)=exc(:)-exci_rpa(:) 1617 vxcrho(:,:)=vxcrho(:,:)-vxci_rpa(:,:) 1618 ABI_DEALLOCATE(exci_rpa) 1619 ABI_DEALLOCATE(vxci_rpa) 1620 end if 1621 1622 !>>>>> LDA xc energy like ixc==7, and Leeuwen-Baerends GGA xc potential 1623 else if(ixc==13) then 1624 if (order**2 <= 1) then 1625 optpbe=1 1626 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 1627 call xclb(grho2_updn,npts,nspden,rho_updn,vxcrho) 1628 else if (order /=3) then 1629 optpbe=1 1630 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,dvxci=dvxc) 1631 call xclb(grho2_updn,npts,nspden,rho_updn,vxcrho) 1632 else if (order ==3) then 1633 optpbe=1 1634 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,d2vxci=d2vxc,dvxci=dvxc) 1635 call xclb(grho2_updn,npts,nspden,rho_updn,vxcrho) 1636 end if 1637 1638 !>>>>> HTCH93, HTCH120, HTCH107, HTCH147 1639 else if(ixc==16 .or. ixc==17 .or. ixc==26 .or. ixc==27) then 1640 call xchcth(vxcgrho,exc,grho2_updn,ixc,npts,nspden,order,rho_updn,vxcrho) 1641 1642 !>>>>> Only for test purpose (test various part of MGGA implementation) 1643 else if(ixc==31 .or. ixc==32 .or. ixc==33 .or. ixc==34) then 1644 exc(:)=zero 1645 vxcrho(:,:)=zero 1646 vxcgrho(:,:)=zero 1647 vxclrho(:,:)=zero 1648 vxctau(:,:)=zero 1649 1650 !>>>>> Perdew-Wang LSD is coded in Perdew-Burke-Ernzerhof GGA, with optpbe=1 1651 optpbe=1 1652 select case(ixc) 1653 case (31) 1654 alpha=1.00d0-(1.00d0/1.01d0) 1655 ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc) 1656 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 1657 if (nspden==1) then 1658 ! it should be : exc_tot= exc_spin up + exc_spin down = 2*exc_spin up but this applies to tau and rho (so it cancels) 1659 exc(:)=exc(:)+alpha*tau_updn(:,1)/rho_updn(:,1) 1660 else 1661 do ispden=1,nspden 1662 exc(:)=exc(:)+alpha*tau_updn(:,ispden)/(rho_updn(:,1)+rho_updn(:,2)) 1663 end do 1664 end if 1665 vxctau(:,:)=alpha 1666 case (32) 1667 alpha=0.01d0 1668 ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc) 1669 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 1670 if (nspden==1) then 1671 exc(:)=exc(:)+2.0d0*alpha*lrho_updn(:,1) 1672 vxcrho(:,1) =vxcrho(:,1)+2.0d0*alpha*lrho_updn(:,1) 1673 vxclrho(:,1)=alpha*2.0d0*rho_updn(:,1) 1674 else 1675 do ispden=1,nspden 1676 exc(:)=exc(:)+alpha*lrho_updn(:,ispden) 1677 vxcrho(:,ispden) =vxcrho(:,ispden)+alpha*(lrho_updn(:,1)+lrho_updn(:,2)) 1678 vxclrho(:,ispden)=alpha*(rho_updn(:,1)+rho_updn(:,2)) 1679 end do 1680 end if 1681 case (33) 1682 alpha=-0.010d0 1683 ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc) 1684 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 1685 if (nspden==1) then 1686 ! it should be : exc_tot= exc_spin up + exc_spin down = 2*exc_spin up but this applies to grho2 and rho 1687 ! (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) 1688 exc(:)=exc(:)+alpha*2.0d0*grho2_updn(:,1)/rho_updn(:,1) 1689 if(nvxcgrho==2)vxcgrho(:,1:2)=2.0d0*alpha 1690 if(nvxcgrho==3)vxcgrho(:,3)=2.0d0*alpha 1691 else 1692 exc(:)=exc(:)+alpha*grho2_updn(:,3)/(rho_updn(:,1)+rho_updn(:,2)) 1693 if(nvxcgrho==2)vxcgrho(:,1:2)=2.0d0*alpha 1694 if(nvxcgrho==3)vxcgrho(:,3)=2.0d0*alpha 1695 end if 1696 case (34) 1697 alpha=-0.010d0 1698 ! Compute first LDA XC (exc,vxc) and then add fake MGGA XC (exc,vxc) 1699 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc) 1700 if (nspden==1) then 1701 exc(:)=exc(:)+16.0d0*alpha*tau_updn(:,1) 1702 vxcrho(:,1)=vxcrho(:,1)+16.0d0*alpha*tau_updn(:,1) 1703 vxctau(:,1)=16.0d0*alpha*rho_updn(:,1) 1704 else 1705 do ispden=1,nspden 1706 exc(:)=exc(:)+8.0d0*alpha*tau_updn(:,ispden) 1707 vxcrho(:,ispden)=vxcrho(:,ispden)+8.0d0*alpha*(tau_updn(:,1)+tau_updn(:,2)) 1708 vxctau(:,ispden)=8.0d0*alpha*(rho_updn(:,1)+rho_updn(:,2)) 1709 end do 1710 end if 1711 end select 1712 1713 !>>>>> Hybrid PBE0 (1/4 and 1/3) 1714 else if(ixc>=41.and.ixc<=42) then 1715 ! Requires to evaluate exchange-correlation with PBE (optpbe=2) 1716 ! minus hyb_mixing*exchange with PBE (optpbe=-2) 1717 ndvxc_x=8 1718 ABI_ALLOCATE(exc_x,(npts)) 1719 ABI_ALLOCATE(vxcrho_x,(npts,nspden)) 1720 ABI_ALLOCATE(vxcgrho_x,(npts,nvxcgrho)) 1721 exc_x=zero;vxcrho_x=zero;vxcgrho_x=zero 1722 if (order**2 <= 1) then 1723 optpbe=2 !PBE exchange correlation 1724 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 1725 & dvxcdgr=vxcgrho,exexch=exexch_,grho2_updn=grho2_updn) 1726 optpbe=-2 !PBE exchange-only 1727 call xcpbe(exc_x,npts,nspden,optpbe,order,rho_updn,vxcrho_x,ndvxc,ngr2,nd2vxc,& 1728 & dvxcdgr=vxcgrho_x,exexch=exexch_,grho2_updn=grho2_updn) 1729 exc=exc-exc_x*hyb_mixing 1730 vxcrho=vxcrho-vxcrho_x*hyb_mixing 1731 vxcgrho=vxcgrho-vxcgrho_x*hyb_mixing 1732 else if (order /=3) then 1733 ABI_ALLOCATE(dvxc_x,(npts,ndvxc_x)) 1734 optpbe=2 !PBE exchange correlation 1735 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 1736 dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 1737 optpbe=-2 !PBE exchange-only 1738 call xcpbe(exc_x,npts,nspden,optpbe,order,rho_updn,vxcrho_x,ndvxc_x,ngr2,nd2vxc,& 1739 & dvxcdgr=vxcgrho_x,dvxci=dvxc_x,grho2_updn=grho2_updn) 1740 exc=exc-exc_x*hyb_mixing 1741 vxcrho=vxcrho-vxcrho_x*hyb_mixing 1742 vxcgrho=vxcgrho-vxcgrho_x*hyb_mixing 1743 dvxc(:,1:ndvxc_x)=dvxc(:,1:ndvxc_x)-dvxc_x(:,1:ndvxc_x)*hyb_mixing 1744 ABI_DEALLOCATE(dvxc_x) 1745 else if (order ==3) then 1746 ! The size of exchange-correlation with PBE (optpbe=2) 1747 ! is the one which defines the size for ndvxc. 1748 ABI_ALLOCATE(dvxc_x,(npts,ndvxc_x)) 1749 ABI_ALLOCATE(d2vxc_x,(npts,nd2vxc)) 1750 optpbe=2 !PBE exchange correlation 1751 call xcpbe(exc,npts,nspden,optpbe,order,rho_updn,vxcrho,ndvxc,ngr2,nd2vxc,& 1752 & d2vxci=d2vxc,dvxcdgr=vxcgrho,dvxci=dvxc,grho2_updn=grho2_updn) 1753 optpbe=-2 !PBE exchange-only 1754 call xcpbe(exc_x,npts,nspden,optpbe,order,rho_updn,vxcrho_x,ndvxc_x,ngr2,nd2vxc,& 1755 & d2vxci=d2vxc_x,dvxcdgr=vxcgrho_x,dvxci=dvxc_x,grho2_updn=grho2_updn) 1756 exc=exc-exc_x*hyb_mixing 1757 vxcrho=vxcrho-vxcrho_x*hyb_mixing 1758 vxcgrho=vxcgrho-vxcgrho_x*hyb_mixing 1759 d2vxc=d2vxc-d2vxc_x*hyb_mixing 1760 dvxc(:,1:ndvxc_x)=dvxc(:,1:ndvxc_x)-dvxc_x(:,1:ndvxc_x)*hyb_mixing 1761 ABI_DEALLOCATE(dvxc_x) 1762 ABI_DEALLOCATE(d2vxc_x) 1763 end if 1764 ABI_DEALLOCATE(exc_x) 1765 ABI_DEALLOCATE(vxcrho_x) 1766 ABI_DEALLOCATE(vxcgrho_x) 1767 1768 !>>>>> Ichimaru,Iyetomi,Tanaka, XC at finite temp (e- gaz) 1769 else if (ixc==50) then 1770 if (order**2 <= 1) then 1771 call xciit(exc,fxcT,npts,order,rspts,el_temp,vxcrho(:,1)) 1772 else 1773 call xciit(exc,fxcT,npts,order,rspts,el_temp,vxcrho(:,1),dvxc) 1774 end if 1775 1776 !>>>>> GGA counterpart of the B3LYP functional 1777 else if(ixc==1402000) then 1778 ! Requires to evaluate exchange-correlation 1779 ! with 5/4 B3LYP - 1/4 B3LYPc, where 1780 ! B3LYPc = (0.19 Ec VWN3 + 0.81 Ec LYP) 1781 1782 ! First evaluate B3LYP. 1783 if(present(xc_funcs))then 1784 if (abs(order)==1) then 1785 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1786 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,xc_functionals=xc_funcs) 1787 elseif (abs(order)==2) then 1788 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1789 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,xc_functionals=xc_funcs) 1790 else 1791 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1792 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,d2vxc=d2vxc,xc_functionals=xc_funcs) 1793 end if 1794 else 1795 if (abs(order)==1) then 1796 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1797 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho) 1798 elseif (abs(order)==2) then 1799 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1800 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc) 1801 else 1802 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1803 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,d2vxc=d2vxc) 1804 end if 1805 end if 1806 1807 ! Then renormalize B3LYP and subtract VWN3 contribution 1808 ABI_ALLOCATE(exc_c,(npts)) 1809 ABI_ALLOCATE(vxcrho_c,(npts,nspden)) 1810 if(order**2>1)then 1811 ABI_ALLOCATE(dvxc_c,(npts,ndvxc)) 1812 end if 1813 if(order**2>4)then 1814 ABI_ALLOCATE(d2vxc_c,(npts,nd2vxc)) 1815 end if 1816 exc_c=zero;vxcrho_c=zero 1817 call libxc_functionals_init(-30,nspden,xc_functionals=xc_funcs_vwn3) 1818 if (order**2 <= 1) then 1819 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 1820 & vxcrho_c,xc_functionals=xc_funcs_vwn3) 1821 elseif (order**2 <= 4) then 1822 dvxc_c=zero 1823 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 1824 & vxcrho_c,dvxc=dvxc_c,xc_functionals=xc_funcs_vwn3) 1825 else 1826 dvxc_c=zero 1827 d2vxc_c=zero 1828 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 1829 & vxcrho_c,dvxc=dvxc_c,d2vxc=d2vxc,xc_functionals=xc_funcs_vwn3) 1830 end if 1831 exc=1.25d0*exc-quarter*0.19d0*exc_c 1832 vxcrho=1.25d0*vxcrho-quarter*0.19d0*vxcrho_c 1833 if(order**2>1)dvxc=1.25d0*dvxc-quarter*0.19d0*dvxc_c 1834 if(order**2>4)d2vxc=1.25d0*d2vxc-quarter*0.19d0*d2vxc_c 1835 call libxc_functionals_end(xc_functionals=xc_funcs_vwn3) 1836 1837 ! Then subtract LYP contribution 1838 call libxc_functionals_init(-131,nspden,xc_functionals=xc_funcs_lyp) 1839 if (order**2 <= 1) then 1840 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 1841 & vxcrho_c,grho2=grho2_updn,vxcgr=vxcgrho,xc_functionals=xc_funcs_lyp) 1842 elseif (order**2 <= 4) then 1843 dvxc_c=zero 1844 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 1845 & vxcrho_c,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc_c,xc_functionals=xc_funcs_lyp) 1846 else 1847 dvxc_c=zero 1848 d2vxc_c=zero 1849 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc_c,& 1850 & vxcrho_c,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc_c,d2vxc=d2vxc,xc_functionals=xc_funcs_lyp) 1851 end if 1852 exc=exc-quarter*0.81d0*exc_c 1853 vxcrho=vxcrho-quarter*0.81d0*vxcrho_c 1854 if(order**2>1)dvxc=dvxc-quarter*0.81d0*dvxc_c 1855 if(order**2>4)d2vxc=d2vxc-quarter*0.81d0*d2vxc_c 1856 call libxc_functionals_end(xc_functionals=xc_funcs_lyp) 1857 1858 ABI_DEALLOCATE(exc_c) 1859 ABI_DEALLOCATE(vxcrho_c) 1860 if(allocated(dvxc_c))then 1861 ABI_DEALLOCATE(dvxc_c) 1862 end if 1863 if(allocated(d2vxc_c))then 1864 ABI_DEALLOCATE(d2vxc_c) 1865 end if 1866 1867 !>>>>> All libXC functionals 1868 else if( ixc<0 ) then 1869 if (is_mgga) then 1870 if(present(xc_funcs))then 1871 if (PRESENT(xc_tb09_c)) then 1872 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1873 & vxcrho,grho2_updn,vxcgrho,lrho_updn,vxclrho,tau_updn,vxctau,xc_functionals=xc_funcs,xc_tb09_c=xc_tb09_c) 1874 else 1875 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1876 & vxcrho,grho2_updn,vxcgrho,lrho_updn,vxclrho,tau_updn,vxctau,xc_functionals=xc_funcs) 1877 end if 1878 else 1879 if (PRESENT(xc_tb09_c)) then 1880 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1881 & vxcrho,grho2_updn,vxcgrho,lrho_updn,vxclrho,tau_updn,vxctau,xc_tb09_c=xc_tb09_c) 1882 else 1883 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1884 & vxcrho,grho2_updn,vxcgrho,lrho_updn,vxclrho,tau_updn,vxctau) 1885 end if 1886 end if 1887 ixc1 = (-ixc)/1000 1888 ixc2 = (-ixc) - ixc1*1000 1889 if(ixc1==206 .or. ixc1==207 .or. ixc1==208 .or. ixc1==209 .or. & 1890 & ixc2==206 .or. ixc2==207 .or. ixc2==208 .or. ixc2==209 )then 1891 ! Assume that that type of mGGA can only be used with a LDA correlation (see doc) 1892 vxcgrho(:,:)=zero 1893 vxclrho(:,:)=zero 1894 vxctau(:,:)=zero 1895 end if 1896 elseif (is_gga) then 1897 if(present(xc_funcs))then 1898 if (order**2 <= 1) then 1899 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1900 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,xc_functionals=xc_funcs) 1901 else 1902 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1903 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc,xc_functionals=xc_funcs) 1904 end if 1905 else 1906 if (order**2 <= 1) then 1907 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1908 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho) 1909 else 1910 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1911 & vxcrho,grho2=grho2_updn,vxcgr=vxcgrho,dvxc=dvxc) 1912 end if 1913 end if 1914 else 1915 if(present(xc_funcs))then 1916 if (order**2 <= 1) then 1917 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1918 & vxcrho,xc_functionals=xc_funcs) 1919 elseif (order**2 <= 4) then 1920 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1921 & vxcrho,dvxc=dvxc,xc_functionals=xc_funcs) 1922 else 1923 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1924 & vxcrho,dvxc=dvxc,d2vxc=d2vxc,xc_functionals=xc_funcs) 1925 end if 1926 else 1927 if (order**2 <= 1) then 1928 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1929 & vxcrho) 1930 elseif (order**2 <= 4) then 1931 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1932 & vxcrho,dvxc=dvxc) 1933 else 1934 call libxc_functionals_getvxc(ndvxc,nd2vxc,npts,nspden,order,rho_updn,exc,& 1935 & vxcrho,dvxc=dvxc,d2vxc=d2vxc) 1936 end if 1937 end if 1938 end if 1939 1940 end if 1941 1942 ! ================================================= 1943 ! == Finalization == 1944 ! ================================================= 1945 !Deallocate arrays 1946 if(allocated(rhotot)) then 1947 ABI_DEALLOCATE(rhotot) 1948 end if 1949 if(allocated(rspts)) then 1950 ABI_DEALLOCATE(rspts) 1951 end if 1952 if(allocated(zeta)) then 1953 ABI_DEALLOCATE(zeta) 1954 end if 1955 1956 end subroutine drivexc
ABINIT/drivexc_main [ Functions ]
NAME
drivexc_main
FUNCTION
Driver of XC functionals. Optionally, deliver the XC kernel, or even the derivative of the XC kernel (the third derivative of the XC energy)
INPUTS
ixc=index of the XC functional mgga=1 if functional is metaGGA ndvxc=size of dvxc(npts,ndvxc) nd2vxc=size of d2vxc(npts,nd2vxc) ngr2=size of grho2(npts,ngr2) npts=number of real space points on which the density (and its gradients) is provided nspden=number of spin-density components (1 or 2) nvxcgrho=size of vxcgrho(npts,nvxcgrho) order=gives the maximal derivative of Exc computed. 1=usual value (return exc and vxc) 2=also computes the kernel (return exc,vxc,kxc) -2=like 2, except (to be described) 3=also computes the derivative of the kernel (return exc,vxc,kxc,k3xc) rho(npts,nspden)=the spin-up and spin-down densities If nspden=1, only the spin-up density must be given. In the calling routine, the spin-down density must be equal to the spin-up density5 and both are half the total density. If nspden=2, the spin-up and spin-down density must be given xclevel= XC functional level === Optional input arguments === [el_temp]= electronic temperature (to be used for finite temperature XC functionals) [exexch]=choice of <<<local>>> exact exchange. Active if exexch=3 (only for GGA, and NOT for libxc) [grho2(npts,ngr2)]=the square of the gradients of spin-up, spin-down, and total density If nspden=1, only the square of the gradient of the spin-up density must be given. In the calling routine, the square of the gradient of the spin-down density must be equal to the square of the gradient of the spin-up density, and both must be equal to one-quarter of the square of the gradient of the total density. If nspden=2, the square of the gradients of spin-up, spin-down, and total density must be given. Note that the square of the gradient of the total density is usually NOT related to the square of the gradient of the spin-up and spin-down densities, because the gradients are not usually aligned. This is not the case when nspden=1. [hyb_mixing]= mixing parameter for the native PBEx functionals (ixc=41 and 42) [lrho(npts,nspden)]=the Laplacian of spin-up and spin-down densities If nspden=1, only the spin-up Laplacian density must be given. In the calling routine, the spin-down Laplacian density must be equal to the spin-up Laplacian density, and both are half the total Laplacian density. If nspden=2, the Laplacian of spin-up and spin-down densities must be given [tau(npts,nspden)]=the spin-up and spin-down kinetic energy densities If nspden=1, only the spin-up kinetic energy density must be given. In the calling routine, the spin-down kinetic energy density must be equal to the spin-up kinetic energy density, and both are half the total kinetic energy density. If nspden=2, the spin-up and spin-down kinetic energy densities must be given [xc_funcs(2)]= <type(libxc_functional_type)>, optional : libxc XC functionals. [xc_tb09_c]=c parameter for the mgga TB09 functional, within libxc
OUTPUT
exc(npts)=exchange-correlation energy density (hartree) vxcrho(npts,nspden)= (d($\rho$*exc)/d($\rho_up$)) (hartree) and (d($\rho$*exc)/d($\rho_down$)) (hartree) === Optional output arguments === [dvxc]=partial second derivatives of the xc energy, only if abs(order)>1 In case of local energy functional (option=1,-1 or 3): dvxc(npts,1+nspden)= (Hartree*bohr^3) if(nspden=1 .and. order==2): dvxci(:,1)=dvxc/d$\rho$ , dvxc(:,2) empty if(nspden=1 .and. order==-2): also compute dvxci(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$ if(nspden=2): dvxc(:,1)=dvxc($\uparrow$)/d$\rho(\downarrow)$, dvxc(:,2)=dvxc($\uparrow$)/d$\rho(\downarrow)$, dvxc(:,3)=dvxc($\downarrow$)/d$\rho(\downarrow)$ In case of gradient corrected functional (option=2,-2, 4, -4, 5, 6, 7): dvxc(npts,15)= dvxc(:,1)= d2Ex/drho_up drho_up dvxc(:,2)= d2Ex/drho_dn drho_dn dvxc(:,3)= dEx/d(abs(grad(rho_up))) / abs(grad(rho_up)) dvxc(:,4)= dEx/d(abs(grad(rho_dn))) / abs(grad(rho_dn)) dvxc(:,5)= d2Ex/d(abs(grad(rho_up))) drho_up / abs(grad(rho_up)) dvxc(:,6)= d2Ex/d(abs(grad(rho_dn))) drho_dn / abs(grad(rho_dn)) dvxc(:,7)= 1/abs(grad(rho_up)) * d/d(abs(grad(rho_up)) (dEx/d(abs(grad(rho_up))) /abs(grad(rho_up))) dvxc(:,8)= 1/abs(grad(rho_dn)) * d/d(abs(grad(rho_dn)) (dEx/d(abs(grad(rho_dn))) /abs(grad(rho_dn))) dvxc(:,9)= d2Ec/drho_up drho_up dvxc(:,10)=d2Ec/drho_up drho_dn dvxc(:,11)=d2Ec/drho_dn drho_dn dvxc(:,12)=dEc/d(abs(grad(rho))) / abs(grad(rho)) dvxc(:,13)=d2Ec/d(abs(grad(rho))) drho_up / abs(grad(rho)) dvxc(:,14)=d2Ec/d(abs(grad(rho))) drho_dn / abs(grad(rho)) dvxc(:,15)=1/abs(grad(rho)) * d/d(abs(grad(rho)) (dEc/d(abs(grad(rho))) /abs(grad(rho))) [d2vxc]=second derivative of the XC potential=3rd order derivative of energy, only if abs(order)>2 (only available for LDA and nspden=1) if nspden=1 d2vxc(npts,1)=second derivative of the XC potential=3rd order derivative of energy if nspden=2 d2vxc(npts,1), d2vxc(npts,2), d2vxc(npts,3), d2vxc(npts,4) (3rd derivative of energy) [fxcT(npts)]=XC free energy of the electron gaz at finite temperature (to be used for plasma systems) [vxcgrho(npts,3)]= 1/$|grad \rho_up|$ (d($\rho$*exc)/d($|grad \rho_up|$)) (hartree) 1/$|grad \rho_dn|$ (d($\rho$*exc)/d($|grad \rho_dn|$)) (hartree) and 1/$|grad \rho|$ (d($\rho$*exc)/d($|grad \rho|$)) (hartree) (will be zero if a LDA functional is used) [vxclrho(npts,nspden)]=(only for meta-GGA, i.e. optional output)= (d($\rho$*exc)/d($\lrho_up$)) (hartree) and (d($\rho$*exc)/d($\lrho_down$)) (hartree) [vxctau(npts,nspden)]=(only for meta-GGA, i.e. optional output)= derivative of XC energy density with respect to kinetic energy density (depsxcdtau). (d($\rho$*exc)/d($\tau_up$)) (hartree) and (d($\rho$*exc)/d($\tau_down$)) (hartree)
PARENTS
m_pawxc,rhotoxc
CHILDREN
drivexc
SOURCE
176 subroutine drivexc_main(exc,ixc,mgga,ndvxc,nd2vxc,ngr2,npts,nspden,nvxcgrho,order,rho,vxcrho,xclevel, & 177 & dvxc,d2vxc,el_temp,exexch,fxcT,grho2,& ! Optional arguments 178 & hyb_mixing,lrho,tau,vxcgrho,vxclrho,vxctau,xc_funcs,xc_tb09_c) ! Optional arguments 179 180 181 !This section has been created automatically by the script Abilint (TD). 182 !Do not modify the following lines by hand. 183 #undef ABI_FUNC 184 #define ABI_FUNC 'drivexc_main' 185 !End of the abilint section 186 187 implicit none 188 189 !Arguments ------------------------------------ 190 !scalars 191 integer,intent(in) :: ixc,mgga,ndvxc,nd2vxc,ngr2,npts,nspden,nvxcgrho,order,xclevel 192 integer,intent(in),optional :: exexch 193 real(dp),intent(in),optional :: el_temp,hyb_mixing,xc_tb09_c 194 !arrays 195 real(dp),intent(in) :: rho(npts,nspden) 196 real(dp),intent(in),optional :: grho2(npts,ngr2),lrho(npts,nspden*mgga),tau(npts,nspden*mgga) 197 real(dp),intent(out) :: exc(npts),vxcrho(npts,nspden) 198 real(dp),intent(out),optional :: dvxc(npts,ndvxc),d2vxc(npts,nd2vxc),fxcT(:),vxcgrho(npts,nvxcgrho) 199 real(dp),intent(out),optional :: vxclrho(npts,nspden*mgga),vxctau(npts,nspden*mgga) 200 type(libxc_functional_type),intent(inout),optional :: xc_funcs(2) 201 202 !Local variables------------------------------- 203 !scalars 204 real(dp) :: hyb_mixing_,xc_tb09_c_ 205 logical :: is_gga 206 207 ! ************************************************************************* 208 209 !Checks input parameters 210 if (mgga==1) then 211 if (.not.present(lrho)) then 212 MSG_BUG('lrho arg must be present in case of mGGA!') 213 end if 214 if (.not.present(tau)) then 215 MSG_BUG('tau arg must be present in case of mGGA!') 216 end if 217 if (.not.present(vxclrho)) then 218 MSG_BUG('vxclrho arg must be present in case of mGGA!') 219 end if 220 if (.not.present(vxctau)) then 221 MSG_BUG('vxctau arg must be present in case of mGGA!') 222 end if 223 end if 224 if (present(fxcT)) then 225 if (.not.present(el_temp)) then 226 MSG_BUG('el_temp arg must be present together with fxcT!') 227 end if 228 end if 229 230 xc_tb09_c_=99.99_dp;if (present(xc_tb09_c)) xc_tb09_c_=xc_tb09_c 231 232 if(ixc==41)hyb_mixing_=quarter 233 if(ixc==42)hyb_mixing_=third 234 if (present(hyb_mixing)) hyb_mixing_=hyb_mixing 235 236 !>>>>> All libXC functionals 237 238 if (ixc<0) then 239 240 if (present(xc_funcs))then 241 is_gga=libxc_functionals_isgga(xc_functionals=xc_funcs) 242 else 243 is_gga=libxc_functionals_isgga() 244 end if 245 246 if (mgga==1) then 247 if (abs(xc_tb09_c_-99.99_dp)>tol12) then 248 if (present(xc_funcs)) then 249 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 250 & lrho_updn=lrho,vxclrho=vxclrho,tau_updn=tau,vxctau=vxctau, & 251 & xc_funcs=xc_funcs,xc_tb09_c=xc_tb09_c_) 252 else 253 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 254 & grho2_updn=grho2,vxcgrho=vxcgrho, & 255 & lrho_updn=lrho,vxclrho=vxclrho,tau_updn=tau,vxctau=vxctau, & 256 & xc_tb09_c=xc_tb09_c_) 257 end if 258 else 259 if (present(xc_funcs)) then 260 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 261 & grho2_updn=grho2,vxcgrho=vxcgrho, & 262 & lrho_updn=lrho,vxclrho=vxclrho,tau_updn=tau,vxctau=vxctau, & 263 & xc_funcs=xc_funcs) 264 else 265 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 266 & grho2_updn=grho2,vxcgrho=vxcgrho, & 267 & lrho_updn=lrho,vxclrho=vxclrho,tau_updn=tau,vxctau=vxctau) 268 end if 269 end if 270 else if (is_gga) then 271 if (order**2<=1) then 272 if (present(xc_funcs)) then 273 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 274 & grho2_updn=grho2,vxcgrho=vxcgrho,xc_funcs=xc_funcs) 275 else 276 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 277 & grho2_updn=grho2,vxcgrho=vxcgrho) 278 end if 279 else 280 if (present(xc_funcs)) then 281 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 282 & grho2_updn=grho2,vxcgrho=vxcgrho,dvxc=dvxc,xc_funcs=xc_funcs) 283 else 284 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 285 & grho2_updn=grho2,vxcgrho=vxcgrho,dvxc=dvxc) 286 end if 287 end if 288 else 289 if (order**2<=1) then 290 if (present(xc_funcs)) then 291 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 292 xc_funcs=xc_funcs) 293 else 294 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho) 295 end if 296 else if (order**2<=4) then 297 if (present(xc_funcs)) then 298 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 299 & dvxc=dvxc, xc_funcs=xc_funcs) 300 else 301 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 302 & dvxc=dvxc) 303 end if 304 else 305 if (present(xc_funcs)) then 306 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 307 & dvxc=dvxc,d2vxc=d2vxc, xc_funcs=xc_funcs) 308 else 309 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 310 & dvxc=dvxc,d2vxc=d2vxc) 311 end if 312 end if 313 end if 314 315 else 316 317 ! Cases with gradient 318 if (xclevel==2)then 319 if (order**2<=1.or.ixc==16.or.ixc==17.or.ixc==26.or.ixc==27) then 320 if (ixc/=13) then 321 if (present(exexch)) then 322 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 323 & hyb_mixing=hyb_mixing_,grho2_updn=grho2,vxcgrho=vxcgrho, & 324 & exexch=exexch) 325 else 326 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 327 & hyb_mixing=hyb_mixing_,grho2_updn=grho2,vxcgrho=vxcgrho) 328 end if 329 else 330 if (present(exexch)) then 331 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 332 & grho2_updn=grho2, & 333 & exexch=exexch) 334 else 335 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 336 & grho2_updn=grho2) 337 end if 338 end if 339 else if (order/=3) then 340 if (ixc/=13) then 341 if (present(exexch)) then 342 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 343 & hyb_mixing=hyb_mixing_,dvxc=dvxc,grho2_updn=grho2,vxcgrho=vxcgrho, & 344 & exexch=exexch) 345 else 346 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 347 & hyb_mixing=hyb_mixing_,dvxc=dvxc,grho2_updn=grho2,vxcgrho=vxcgrho) 348 end if 349 else 350 if (present(exexch)) then 351 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 352 & dvxc=dvxc,grho2_updn=grho2, & 353 & exexch=exexch) 354 else 355 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 356 & dvxc=dvxc,grho2_updn=grho2) 357 end if 358 end if 359 else if (order==3) then 360 if (ixc/=13) then 361 if (present(exexch)) then 362 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 363 & dvxc=dvxc,d2vxc=d2vxc,grho2_updn=grho2,vxcgrho=vxcgrho, & 364 & hyb_mixing=hyb_mixing_,exexch=exexch) 365 else 366 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 367 & hyb_mixing=hyb_mixing_,dvxc=dvxc,d2vxc=d2vxc,grho2_updn=grho2,vxcgrho=vxcgrho) 368 end if 369 else 370 if (present(exexch)) then 371 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 372 & dvxc=dvxc,d2vxc=d2vxc,grho2_updn=grho2, & 373 & exexch=exexch) 374 else 375 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 376 & dvxc=dvxc,d2vxc=d2vxc,grho2_updn=grho2) 377 end if 378 end if 379 end if 380 381 382 ! Cases without gradient 383 else 384 if (order**2<=1) then 385 if (ixc>=31.and.ixc<=34) then !fake mgga functionals for testing purpose only (based on LDA functional) 386 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 387 & grho2_updn=grho2,vxcgrho=vxcgrho, & 388 & lrho_updn=lrho,vxclrho=vxclrho,tau_updn=tau,vxctau=vxctau) 389 else 390 if (present(fxcT)) then 391 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho,fxcT=fxcT,el_temp=el_temp) 392 else 393 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho) 394 end if 395 end if 396 else if (order==3.and.(ixc==3.or.ixc>=7.and.ixc<=10)) then 397 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 398 & dvxc=dvxc,d2vxc=d2vxc) 399 else 400 if (present(fxcT)) then 401 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 402 & dvxc=dvxc,fxcT=fxcT,el_temp=el_temp) 403 else 404 call drivexc(exc,ixc,npts,nspden,order,rho,vxcrho,ndvxc,ngr2,nd2vxc,nvxcgrho, & 405 & dvxc=dvxc) 406 end if 407 end if 408 end if 409 410 end if 411 412 end subroutine drivexc_main
ABINIT/echo_xc_name [ 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
PARENTS
driver
CHILDREN
wrtout
SOURCE
433 subroutine echo_xc_name (ixc) 434 435 436 !This section has been created automatically by the script Abilint (TD). 437 !Do not modify the following lines by hand. 438 #undef ABI_FUNC 439 #define ABI_FUNC 'echo_xc_name' 440 !End of the abilint section 441 442 implicit none 443 444 !Arguments ------------------------------- 445 integer, intent(in) :: ixc 446 447 !Local variables ------------------------- 448 integer :: l_citation 449 character(len=500) :: message, citation 450 451 ! ********************************************************************* 452 453 message ='' 454 citation ='' 455 456 !normal case (not libxc) 457 if (ixc >= 0) then 458 459 select case (ixc) 460 case (0) 461 message = 'No xc applied (usually for testing) - ixc=0' 462 citation = '' 463 ! LDA,LSD 464 case (1) 465 message = 'LDA: new Teter (4/93) with spin-polarized option - ixc=1' 466 citation = 'S. Goedecker, M. Teter, J. Huetter, PRB 54, 1703 (1996)' ! [[cite:Goedecker1996]] 467 case (2) 468 message = 'LDA: Perdew-Zunger-Ceperley-Alder - ixc=2' 469 citation = 'J.P.Perdew and A.Zunger, PRB 23, 5048 (1981) ' ! [[cite:Perdew1981]] 470 case (3) 471 message = 'LDA: old Teter (4/91) fit to Ceperley-Alder data - ixc=3' 472 citation = '' 473 case (4) 474 message = 'LDA: Wigner - ixc=4' 475 citation = 'E.P.Wigner, Trans. Faraday Soc. 34, 67 (1938)' ! [[cite:Wigner1938]] 476 case (5) 477 message = 'LDA: Hedin-Lundqvist - ixc=5' 478 citation = 'L.Hedin and B.I.Lundqvist, J. Phys. C4, 2064 (1971)' ! [[cite:Hedin1971]] 479 case (6) 480 message = 'LDA: "X-alpha" xc - ixc=6' 481 citation = 'Slater J. C., Phys. Rev. 81, 385 (1951)' ! [[cite:Slater1951]] 482 case (7) 483 message = 'LDA: Perdew-Wang 92 LSD fit to Ceperley-Alder data - ixc=7' 484 citation = 'J.P.Perdew and Y.Wang, PRB 45, 13244 (1992)' ! [[cite:Perdew1992a]] 485 case (8) 486 message = 'LDA: Perdew-Wang 92 LSD , exchange-only - ixc=8' 487 citation = 'J.P.Perdew and Y.Wang, PRB 45, 13244 (1992)' ! [[cite:Perdew1992a]] 488 case (9) 489 message = 'LDA: Perdew-Wang 92 Ex+Ec_RPA energy - ixc=9' 490 citation = 'J.P.Perdew and Y.Wang, PRB 45, 13244 (1992)' ! [[cite:Perdew1992]] 491 case (10) 492 message = 'LDA: RPA LSD energy (only the energy !!) - ixc=10' 493 citation = '' 494 ! GGA 495 case (11) 496 message = 'GGA: Perdew-Burke-Ernzerhof functional - ixc=11' 497 citation = 'J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996)' ! [[cite:Perdew1996]] 498 case (12) 499 message = 'GGA: x-only Perdew-Burke-Ernzerhof functional - ixc=12' 500 citation = 'J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996)' ! [[cite:Perdew1996]] 501 case (13) 502 message = 'GGA: LDA (ixc==7) energy, and the xc _potential_ is given by van Leeuwen-Baerends GGA - ixc=13' 503 citation = 'R. van Leeuwen and E. J. Baerends PRA 49, 2421 (1994)' ! [[cite:VanLeeuwen1994]] 504 case (14) 505 message = 'GGA: revPBE functional - ixc=14' 506 citation = 'Zhang and Yang, PRL 80, 890 (1998)' ! [[cite:Zhang1998]] 507 case (15) 508 message = 'GGA: RPBE functional - ixc=15' 509 citation = 'Hammer, L. B. Hansen, and J. K. Norskov, PRB 59, 7413 (1999)' ! [[cite:Hammer1999]] 510 case (16) 511 message = 'GGA: HCTH93 functional - ixc=16' 512 citation = 'F.A. Hamprecht, A.J. Cohen, D.J. Tozer, N.C. Handy, JCP 109, 6264 (1998)' ! [[cite:Hamprecht1998]] 513 case (17) 514 message = 'GGA: HCTH120 functional - ixc=17' 515 citation = 'A.D. Boese, N.L. Doltsinis, N.C. Handy, and M. Sprik, JCP 112, 1670 (2000)' ! [[cite:Boese2000]] 516 case (23) 517 message = 'GGA: Wu Cohen functional - ixc=23' 518 citation = 'Z. Wu and R. E. Cohen, PRB 73, 235116 (2006)' ! [[cite:Wu2006]] 519 case (24) 520 message = 'GGA: C09x exchange functional - ixc=24' 521 citation = 'Valentino R. Cooper, PRB 81, 161104(R) (2010)' ! [[cite:Cooper2010]] 522 case (26) 523 message = 'GGA: HCTH147 functional - ixc=26' 524 citation = 'A.D. Boese, N.L. Doltsinis, N.C. Handy, and M. Sprik, JCP 112, 1670 (2000)' ! [[cite:Boese2000]] 525 case (27) 526 message = 'GGA: HCTH407 functional - ixc=27' 527 citation = 'A.D. Boese, and N.C. Handy, JCP 114, 5497 (2001)' ! [[cite:Boese2001]] 528 ! Fermi-Amaldi 529 case (20) 530 message = 'Fermi-Amaldi correction - ixc=20' 531 citation = '' 532 case (21) 533 message = 'Fermi-Amaldi correction with LDA(ixc=1) kernel - ixc=21' 534 citation = '' 535 case (22) 536 message = 'Fermi-Amaldi correction with hybrid BPG kernel - ixc=22' 537 citation = '' 538 case (31) 539 message = 'Meta-GGA fake1 - ixc=31' 540 citation = '' 541 case (32) 542 message = 'Meta-GGA fake2 - ixc=32' 543 citation = '' 544 case (33) 545 message = 'Meta-GGA fake3 - ixc=33' 546 citation = '' 547 case (34) 548 message = 'Meta-GGA fake4 - ixc=34' 549 citation = '' 550 case (40) 551 message = 'Hartree-Fock with mixing coefficient alpha=1' 552 citation = '' 553 case (41) 554 message = 'PBE0 with alpha=0.25' 555 citation = '' 556 case (42) 557 message = 'modified PBE0 with alpha=0.33' 558 citation = '' 559 case (50) 560 message = 'LDA at finite T Ichimaru-Iyetomy-Tanaka - ixc=50' 561 citation = 'Ichimaru S., Iyetomi H., Tanaka S., Phys. Rep. 149, 91-205 (1987) ' ! [[cite:Ichimaru1987]] 562 case default 563 write(message,'(a,i0)')" echo_xc_name does not know how to handle ixc = ",ixc 564 MSG_WARNING(message) 565 end select 566 567 message = " Exchange-correlation functional for the present dataset will be:" // ch10 & 568 & // " " // trim(message) 569 570 l_citation=len_trim(citation) 571 citation = " Citation for XC functional:" // ch10 // " " // trim(citation) 572 573 call wrtout(ab_out,message,'COLL') 574 call wrtout(std_out,message,'COLL') 575 576 if(l_citation/=0)then 577 call wrtout(ab_out,citation,'COLL') 578 call wrtout(std_out,citation,'COLL') 579 end if 580 581 message =' ' 582 call wrtout(ab_out,message,'COLL') 583 call wrtout(std_out,message,'COLL') 584 585 end if ! end libxc if 586 587 end subroutine echo_xc_name
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-2018 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 .
PARENTS
CHILDREN
SOURCE
22 #if defined HAVE_CONFIG_H 23 #include "config.h" 24 #endif 25 26 #include "abi_common.h" 27 28 module m_drivexc 29 30 use defs_basis 31 use m_abicore 32 use m_errors 33 use libxc_functionals 34 35 use m_numeric_tools, only : invcb 36 use m_xciit, only : xciit 37 use m_xcpbe, only : xcpbe 38 use m_xchcth, only : xchcth 39 use m_xclda, only : xcpzca, xcspol, xctetr, xcwign, xchelu, xcxalp, xclb 40 41 implicit none 42 43 private
ABINIT/mkdenpos [ Functions ]
NAME
mkdenpos
FUNCTION
Make a ground-state 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).
PARENTS
bethe_salpeter,m_pawxc,mkcore_wvl,posdoppler,poslifetime,posratecore psolver_rhohxc,rhohxcpositron,rhotoxc,wvl_initro
CHILDREN
SOURCE
968 subroutine mkdenpos(iwarn,nfft,nspden,option,rhonow,xc_denpos) 969 970 971 !This section has been created automatically by the script Abilint (TD). 972 !Do not modify the following lines by hand. 973 #undef ABI_FUNC 974 #define ABI_FUNC 'mkdenpos' 975 !End of the abilint section 976 977 implicit none 978 979 !Arguments ------------------------------------ 980 !scalars 981 integer,intent(in) :: nfft,nspden,option 982 integer,intent(inout) :: iwarn 983 real(dp),intent(in) :: xc_denpos 984 !arrays 985 real(dp),intent(inout) :: rhonow(nfft,nspden) 986 987 !Local variables------------------------------- 988 !scalars 989 integer :: ifft,ispden,numneg 990 real(dp) :: rhotmp,worst 991 character(len=500) :: message 992 !arrays 993 real(dp) :: rho(2) 994 995 ! ************************************************************************* 996 997 numneg=0 998 worst=zero 999 1000 if(nspden==1)then 1001 1002 ! Non spin-polarized 1003 !$OMP PARALLEL DO PRIVATE(ifft,rhotmp) REDUCTION(MIN:worst) REDUCTION(+:numneg) SHARED(nfft,rhonow) 1004 do ifft=1,nfft 1005 rhotmp=rhonow(ifft,1) 1006 if(rhotmp<xc_denpos)then 1007 if(rhotmp<-xc_denpos)then 1008 ! This case is probably beyond machine precision considerations 1009 worst=min(worst,rhotmp) 1010 numneg=numneg+1 1011 end if 1012 rhonow(ifft,1)=xc_denpos 1013 end if 1014 end do 1015 else if (nspden==2) then 1016 1017 ! Spin-polarized 1018 1019 ! rhonow is stored as (up,dn) 1020 if (option==0) then 1021 1022 !$OMP PARALLEL DO PRIVATE(ifft,ispden,rho,rhotmp) REDUCTION(MIN:worst) REDUCTION(+:numneg) & 1023 !$OMP&SHARED(nfft,nspden,rhonow) 1024 do ifft=1,nfft 1025 ! For polarized case, rho(1) is spin-up density, rho(2) is spin-down density 1026 rho(1)=rhonow(ifft,1) 1027 rho(2)=rhonow(ifft,2) 1028 do ispden=1,nspden 1029 if (rho(ispden)<xc_denpos) then 1030 if (rho(ispden)<-xc_denpos) then 1031 ! This case is probably beyond machine precision considerations 1032 worst=min(worst,rho(ispden)) 1033 numneg=numneg+1 1034 end if 1035 rhonow(ifft,ispden)=xc_denpos 1036 end if 1037 end do 1038 end do 1039 1040 ! rhonow is stored as (up+dn,up) 1041 else if (option==1) then 1042 1043 !$OMP PARALLEL DO PRIVATE(ifft,ispden,rho,rhotmp) & 1044 !$OMP&REDUCTION(MIN:worst) REDUCTION(+:numneg) & 1045 !$OMP&SHARED(nfft,nspden,rhonow) 1046 do ifft=1,nfft 1047 ! For polarized case, rho(1) is spin-up density, rho(2) is spin-down density 1048 rho(1)=rhonow(ifft,2) 1049 rho(2)=rhonow(ifft,1)-rho(1) 1050 do ispden=1,nspden 1051 if (rho(ispden)<xc_denpos) then 1052 if (rho(ispden)<-xc_denpos) then 1053 ! This case is probably beyond machine precision considerations 1054 worst=min(worst,rho(ispden)) 1055 numneg=numneg+1 1056 end if 1057 rho(ispden)=xc_denpos 1058 rhonow(ifft,1)=rho(1)+rho(2) 1059 rhonow(ifft,2)=rho(1) 1060 end if 1061 end do 1062 end do 1063 1064 end if ! option 1065 1066 else 1067 MSG_BUG('nspden>2 not allowed !') 1068 end if ! End choice between non-spin polarized and spin-polarized. 1069 1070 if (numneg>0) then 1071 if (iwarn==0) then 1072 write(message,'(a,i0,a,a,a,es10.2,a,e10.2,a,a,a,a)')& 1073 & 'Density went too small (lower than xc_denpos) at ',numneg,' points',ch10,& 1074 & 'and was set to xc_denpos = ',xc_denpos,'. Lowest was ',worst,'.',ch10,& 1075 & 'Likely due to too low boxcut or too low ecut for',' pseudopotential core charge.' 1076 MSG_WARNING(message) 1077 end if 1078 iwarn=iwarn+1 1079 end if 1080 1081 end subroutine mkdenpos
ABINIT/size_dvxc [ Functions ]
NAME
size_dvxc
FUNCTION
Give the size of the array dvxc(npts,ndvxc) and the second dimension of the d2vxc(npts,nd2vxc) needed for the allocations depending on the routine which is called from the drivexc routine
INPUTS
[add_tfw]= optional flag controling the addition of Weiszacker gradient correction to Thomas-Fermi XC energy 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) [xc_funcs(2)]= <type(libxc_functional_type)>
OUTPUT
ndvxc size of the array dvxc(npts,ndvxc) for allocation ngr2 size of the array grho2_updn(npts,ngr2) for allocation nd2vxc size of the array d2vxc(npts,nd2vxc) for allocation nvxcdgr size of the array dvxcdgr(npts,nvxcdgr) for allocation
PARENTS
m_pawxc,rhotoxc
CHILDREN
SOURCE
702 subroutine size_dvxc(ixc,ndvxc,ngr2,nd2vxc,nspden,nvxcdgr,order,& 703 & add_tfw,xc_funcs) ! Optional 704 705 706 !This section has been created automatically by the script Abilint (TD). 707 !Do not modify the following lines by hand. 708 #undef ABI_FUNC 709 #define ABI_FUNC 'size_dvxc' 710 !End of the abilint section 711 712 implicit none 713 714 !Arguments---------------------- 715 integer, intent(in) :: ixc,nspden,order 716 integer, intent(out) :: ndvxc,nd2vxc,ngr2,nvxcdgr 717 logical, intent(in), optional :: add_tfw 718 type(libxc_functional_type),intent(in),optional :: xc_funcs(2) 719 720 !Local variables---------------- 721 logical :: add_tfw_,isgga,ismgga,is_hybrid 722 723 ! ************************************************************************* 724 725 add_tfw_=.false.;if (present(add_tfw)) add_tfw_=add_tfw 726 isgga=.false. ; ismgga=.false. ; is_hybrid=.false. 727 if(ixc<0)then 728 if(present(xc_funcs))then 729 isgga=libxc_functionals_isgga(xc_functionals=xc_funcs) 730 ismgga=libxc_functionals_ismgga(xc_functionals=xc_funcs) 731 is_hybrid=libxc_functionals_is_hybrid(xc_functionals=xc_funcs) 732 else 733 isgga=libxc_functionals_isgga() 734 ismgga=libxc_functionals_ismgga() 735 is_hybrid=libxc_functionals_is_hybrid() 736 end if 737 end if 738 739 ngr2=0 740 nvxcdgr=0 741 ndvxc=0 742 nd2vxc=0 743 744 !Dimension for the gradient of the density (only allocated for GGA or mGGA) 745 if ((ixc>=11.and.ixc<=17).or.(ixc>=23.and.ixc<=24).or.ixc==26.or.ixc==27.or. & 746 & (ixc>=31.and.ixc<=34).or.(ixc==41.or.ixc==42).or.ixc==1402000.or.(add_tfw_)) ngr2=2*min(nspden,2)-1 747 if (ixc<0.and.isgga.or.ismgga.or.is_hybrid) ngr2=2*min(nspden,2)-1 748 749 !A-Only Exc and Vxc 750 !======================================================================================= 751 if (order**2 <= 1) then 752 if (((ixc>=11 .and. ixc<=15) .and. ixc/=13) .or. (ixc>=23 .and. ixc<=24) .or. & 753 & (ixc==41 .or. ixc==42) .or. ixc==1402000) nvxcdgr=3 754 if (ixc==16.or.ixc==17.or.ixc==26.or.ixc==27) nvxcdgr=2 755 if (ixc<0) nvxcdgr=3 756 if (ixc>=31 .and. ixc<=34) nvxcdgr=3 !Native fake metaGGA functionals (for testing purpose only) 757 if (add_tfw_) nvxcdgr=3 758 759 ! B- Exc+Vxc and other derivatives 760 ! ======================================================================================= 761 else 762 763 ! Definition of ndvxc and nvxcdgr, 2nd dimension of the arrays of 2nd-order derivatives 764 ! ------------------------------------------------------------------------------------- 765 if (ixc==1 .or. ixc==21 .or. ixc==22 .or. (ixc>=7 .and. ixc<=10) .or. ixc==13) then 766 ! Routine xcspol: new Teter fit (4/93) to Ceperley-Alder data, with spin-pol option routine xcspol 767 ! Routine xcpbe, with different options (optpbe) and orders (order) 768 ndvxc=min(nspden,2)+1 769 else if (ixc>=2 .and. ixc<=6) then 770 ! Perdew-Zunger fit to Ceperly-Alder data (no spin-pol) !routine xcpzca 771 ! Teter fit (4/91) to Ceperley-Alder values (no spin-pol) !routine xctetr 772 ! Wigner xc (no spin-pol) !routine xcwign 773 ! Hedin-Lundqvist xc (no spin-pol) !routine xchelu 774 ! X-alpha (no spin-pol) !routine xcxalp 775 ndvxc=1 776 else if (ixc==12 .or. ixc==24) then 777 ! Routine xcpbe, with optpbe=-2 and different orders (order) 778 ndvxc=8 779 nvxcdgr=3 780 else if ((ixc>=11 .and. ixc<=15 .and. ixc/=13) .or. ixc==23 .or. ixc==41 .or. ixc==42) then 781 ! Routine xcpbe, with different options (optpbe) and orders (order) 782 ndvxc=15 783 nvxcdgr=3 784 else if(ixc==16 .or. ixc==17 .or. ixc==26 .or. ixc==27 ) then 785 nvxcdgr=2 786 else if (ixc==50) then 787 ndvxc=1 ! IIT xc (no spin-pol) 788 else if (ixc==1402000) then 789 ndvxc=15 790 nvxcdgr=3 791 else if (ixc<0) then 792 if(libxc_functionals_isgga() .or. libxc_functionals_ismgga() .or. libxc_functionals_is_hybrid()) then 793 ndvxc=15 794 else 795 ndvxc=3 796 end if 797 nvxcdgr=3 798 end if 799 if (add_tfw_) nvxcdgr=3 800 801 ! Definition of nd2vxc, 2nd dimension of the array of 3rd-order derivatives 802 ! ------------------------------------------------------------------------------------- 803 if (order==3) then 804 if (ixc==3) nd2vxc=1 ! Non spin polarized LDA case 805 if ((ixc>=7 .and. ixc<=10) .or. (ixc==13)) nd2vxc=3*min(nspden,2)-2 806 ! Following line to be corrected when the calculation of d2vxcar is implemented for these functionals 807 if ((ixc>=11 .and. ixc<=15 .and. ixc/=13) .or. (ixc==23.and.ixc<=24) .or. (ixc==41.or.ixc==42)) nd2vxc=1 808 if(ixc==1402000)nd2vxc=3*min(nspden,2)-2 809 if ((ixc<0.and.(.not.(libxc_functionals_isgga().or.libxc_functionals_ismgga().or.libxc_functionals_is_hybrid() )))) & 810 & nd2vxc=3*min(nspden,2)-2 811 end if 812 813 end if 814 815 end subroutine size_dvxc
ABINIT/xcmult [ 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.
PARENTS
m_pawxc,rhotoxc
CHILDREN
SOURCE
868 subroutine xcmult (depsxc,nfft,ngrad,nspden,nspgrad,rhonow) 869 870 871 !This section has been created automatically by the script Abilint (TD). 872 !Do not modify the following lines by hand. 873 #undef ABI_FUNC 874 #define ABI_FUNC 'xcmult' 875 !End of the abilint section 876 877 implicit none 878 879 !Arguments ------------------------------------ 880 !scalars 881 integer,intent(in) :: nfft,ngrad,nspden,nspgrad 882 !arrays 883 real(dp),intent(in) :: depsxc(nfft,nspgrad) 884 real(dp),intent(inout) :: rhonow(nfft,nspden,ngrad*ngrad) 885 886 !Local variables------------------------------- 887 !scalars 888 integer :: idir,ifft 889 real(dp) :: rho_tot,rho_up 890 891 ! ************************************************************************* 892 893 do idir=1,3 894 895 if(nspden==1)then 896 !$OMP PARALLEL DO PRIVATE(ifft) SHARED(depsxc,idir,nfft,rhonow) 897 do ifft=1,nfft 898 rhonow(ifft,1,1+idir)=rhonow(ifft,1,1+idir)*depsxc(ifft,2) 899 end do 900 901 else 902 903 ! In the spin-polarized case, there are more factors to take into account 904 !$OMP PARALLEL DO PRIVATE(ifft,rho_tot,rho_up) SHARED(depsxc,idir,nfft,rhonow) 905 do ifft=1,nfft 906 rho_tot=rhonow(ifft,1,1+idir) 907 rho_up =rhonow(ifft,2,1+idir) 908 rhonow(ifft,1,1+idir)=rho_up *depsxc(ifft,3) + rho_tot*depsxc(ifft,5) 909 rhonow(ifft,2,1+idir)=(rho_tot-rho_up)*depsxc(ifft,4)+ rho_tot*depsxc(ifft,5) 910 end do 911 912 end if ! nspden==1 913 914 end do ! End loop on directions 915 916 end subroutine xcmult