TABLE OF CONTENTS
m_xc_tb09/m_xc_tb09 [ Modules ]
NAME
m_xc_tb09
FUNCTION
This module contains routine(s) related to exchange-correlation Tran-Blaha 2009 functional (modified Becke-Johnson)
COPYRIGHT
Copyright (C) 2023-2024 ABINIT group (MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 MODULE m_xc_tb09 24 25 use defs_basis 26 use m_abicore 27 use m_xmpi 28 use m_errors 29 use m_time, only : timab 30 use defs_abitypes, only : MPI_type 31 32 use m_xctk, only : xcden 33 use m_drivexc, only : mkdenpos 34 use m_geometry, only : metric 35 36 use m_pawang, only : pawang_type 37 use m_pawrad, only : pawrad_type,pawrad_deducer0,nderiv_gen,simp_gen 38 use m_pawtab, only : pawtab_type 39 use m_pawrhoij, only : pawrhoij_type 40 use m_paw_denpot, only : pawdensities 41 use m_paral_atom, only : get_my_atmtab,free_my_atmtab 42 43 use libxc_functionals, only : libxc_functional_type,libxc_functionals_set_c_tb09, & 44 & libxc_functionals_is_tb09,libxc_functionals_ixc 45 46 implicit none 47 48 private 49 50 !public procedure(s) 51 public :: xc_tb09_update_c ! Compute the C parameter for the TB09 functional (NCPP+PAW) 52 53 !Private variables 54 real(dp),save :: grho_over_rho_pw ! Contribution to TB09 c from the plane-wave pseudo-density 55 real(dp),save :: grho_over_rho_paw ! Contribution to TB09 c from the on-site PAW densities 56 57 CONTAINS !========================================================================================
m_xc_tb09/xc_tb09_update_c [ Functions ]
[ Top ] [ m_xc_tb09 ] [ Functions ]
NAME
xc_tb09_update_c
FUNCTION
This routine computes from the density the C parameter of the Tran-Blaha 2009 XC functional (modified Becke-Johnson) C = 1/V Int[|Grad(Rho)|/Rho] + Sum_PAW_aug_regions[ Natom/V [|Grad(Rho)|/Rho] ]
INPUTS
intxc = 1 if the XC functional has to be interpolated on a more refined mesh ixc= choice of exchange-correlation scheme mpi_enreg= information about MPI parallelization natom= total number of atoms in cell nfft= (effective) number of FFT grid points (for this processor) ngfft(18)= contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nhat(nfft,nspden*nhatdim)= -PAW only- compensation density nhatdim= -PAW only- 0 if nhat array is not used ; 1 otherwise nhatgr(nfft,xcdata%nspden,3*nhatgrdim)= -PAW only- cartesian gradients of compensation density nhatgrdim= -PAW only- 0 if nhatgr array is not used ; 1 otherwise nspden= number of spin-density components ntypat= number of types of atoms in unit cell. n3xccc= dimension of the xccc3d array (0 or nfft or cplx*nfft). pawang <type(pawang_type)> =paw angular mesh and related data pawrad(ntypat) <type(pawrad_type)> =paw radial mesh and related data pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom) pawtab(ntypat) <type(pawtab_type)>= paw tabulated starting data pawxcdev= choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments) rhor(nfft,nspden)= electron density in real space in electrons/bohr**3 rprimd(3,3)=dimensional primitive translations in real space (bohr) usepaw = 1 for PAW computation, 0 else xccc3d(n3xccc)=3D core electron density for XC core correction (bohr^-3) xc_denpos= lowest allowed density ----- Optional arguments ----- [computation_type]= 'only_pw' : update only plane-wave contribution 'only_paw': update only PAW on-site contribution 'all' (default): compute all contributions [mpi_atmtab(:)]= indexes of the atoms treated by current proc [comm_atom]= MPI communicator over atoms [xc_funcs(2)]= <type(libxc_functional_type)>=libxc XC functionals (if not the global one)
OUTPUT
No output The C parameter is directly set in the xc_functional datastructure
SOURCE
110 subroutine xc_tb09_update_c(intxc,ixc,mpi_enreg,natom,nfft,ngfft,nhat,nhatdim, & 111 & nhatgr,nhatgrdim,nspden,ntypat,n3xccc,pawang,pawrad, & 112 & pawrhoij,pawtab,pawxcdev,rhor,rprimd,usepaw,xccc3d,xc_denpos, & 113 & computation_type,mpi_atmtab,comm_atom,xc_funcs) ! Optional arguments 114 115 !Arguments ------------------------------------ 116 !scalars 117 integer,intent(in) :: intxc,ixc,n3xccc,nfft,natom,nhatdim,nhatgrdim 118 integer,intent(in) :: nspden,ntypat,pawxcdev,usepaw 119 integer,intent(in),optional :: comm_atom 120 real(dp),intent(in),optional :: xc_denpos 121 character(len=*),intent(in),optional :: computation_type 122 type(pawang_type),intent(in) :: pawang 123 type(MPI_type),intent(in) :: mpi_enreg 124 !arrays 125 integer,intent(in) :: ngfft(18) 126 integer,intent(in),optional,target :: mpi_atmtab(:) 127 real(dp),intent(in) :: nhat(nfft,nspden*nhatdim),nhatgr(nfft,nspden,3*nhatgrdim) 128 real(dp),intent(in) :: rhor(nfft,nspden) 129 real(dp),intent(in) :: rprimd(3,3),xccc3d(n3xccc) 130 type(libxc_functional_type),intent(inout),optional :: xc_funcs(2) 131 type(pawrad_type),intent(in) :: pawrad(ntypat*usepaw) 132 type(pawrhoij_type),intent(in) :: pawrhoij(:) 133 type(pawtab_type),intent(in),target :: pawtab(ntypat*usepaw) 134 135 !Local variables------------------------------- 136 !scalars 137 integer :: cplex,iatom,iatom_tot,ierr,ii,ilm,iloop,ir,itypat,ixc_from_lib,ipts 138 integer :: ishift,iwarn,lm_size,lm_size_eff,mesh_size,my_comm_atom,my_natom,ngrad 139 integer :: nfftot,npts,nspden1,nzlmopt,option,opt_compch,opt_dens,pawprtvol,usecore 140 integer :: usenhat,usexcnhat 141 real(dp) :: compch_dum,factor,grho_over_rho,grho2,rho,sumg,ucvol,xc_c_tb09 142 logical :: my_atmtab_allocated,paral_atom,test_tb09,use_exact_nhat_gradient 143 logical :: compute_pw,compute_paw 144 character(len=500) :: msg 145 !arrays 146 integer,pointer :: my_atmtab(:) 147 logical,allocatable :: lmselect(:),lmselect_tmp(:) 148 real(dp) :: gmet(3,3),gprimd(3,3),qphon(3),rmet(3,3) 149 !real(dp) :: tsec(2) 150 real(dp),allocatable :: dylmdr(:,:,:),ff(:) 151 real(dp),allocatable :: rhotot(:),grhotot(:,:),rhonow(:,:,:) 152 real(dp),allocatable :: rhoxc(:,:),drhoxc(:),dcorexc(:) 153 real(dp),allocatable :: rho1(:,:,:),trho1(:,:,:),nhat1(:,:,:) 154 real(dp), ABI_CONTIGUOUS pointer :: corexc(:) 155 156 ! ************************************************************************* 157 158 !call timab(xx,1,tsec) 159 160 !Type of calculation 161 compute_pw=.true. ; compute_paw=.true. 162 if (present(computation_type)) then 163 if (computation_type=='only_pw') compute_paw=.false. 164 if (computation_type=='only_paw') compute_pw=.false. 165 end if 166 if ((.not.compute_pw).and.(.not.compute_paw)) return 167 168 !Nothing to do if XC functional is not TB09 169 if (present(xc_funcs)) then 170 test_tb09=libxc_functionals_is_tb09(xc_functionals=xc_funcs) 171 else 172 test_tb09=libxc_functionals_is_tb09() 173 end if 174 if (.not.test_tb09) return 175 176 !Check options 177 if(ixc<0) then 178 if (present(xc_funcs)) then 179 ixc_from_lib=libxc_functionals_ixc(xc_functionals=xc_funcs) 180 else 181 ixc_from_lib=libxc_functionals_ixc() 182 end if 183 if (ixc /= ixc_from_lib) then 184 write(msg, '(a,i0,a,a,i0)')& 185 & 'The value of ixc specified in input, ixc = ',ixc,ch10,& 186 & 'differs from the one used to initialize the functional ',ixc_from_lib 187 ABI_BUG(msg) 188 end if 189 end if 190 if (usepaw==1.and.compute_paw) then 191 if (pawxcdev/=0) then 192 msg='TB09 (mBJ) XC functional only available with pawxcdev=0!' 193 ABI_BUG(msg) 194 end if 195 if(pawang%angl_size==0) then 196 msg='pawang%angl_size=0!' 197 ABI_BUG(msg) 198 end if 199 if(.not.allocated(pawang%ylmr)) then 200 msg='pawang%ylmr must be allocated!' 201 ABI_BUG(msg) 202 end if 203 if(.not.allocated(pawang%ylmrgr)) then 204 msg='pawang%ylmrgr must be allocated!' 205 ABI_BUG(msg) 206 end if 207 if (size(pawrhoij)/=mpi_enreg%my_natom) then 208 msg='Invalid pawrhoij size!' 209 ABI_BUG(msg) 210 end if 211 end if 212 213 !Initialize cell data 214 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 215 216 !Other initializations 217 usexcnhat=0; if (usepaw==1) usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat) 218 219 !Initialization of Int[|Grad(Rho)|/Rho] contributions (plane-wave and on-site if PAW) 220 if (compute_pw) grho_over_rho_pw=zero 221 if (compute_paw) grho_over_rho_paw=zero 222 223 ! ======================================================================== 224 ! =========== Plane wave contribution to Int[|Grad(Rho)|/Rho] ============ 225 226 if (compute_pw) then 227 228 ! Total density used in XC (with/without core density and compensation density) 229 ABI_MALLOC(rhotot,(nfft)) 230 rhotot(1:nfft)=rhor(1:nfft,1) 231 if (usexcnhat==0.and.nhatdim==1) rhotot(1:nfft)=rhotot(1:nfft)-nhat(1:nfft,1) 232 if (n3xccc>0) rhotot(1:nfft)=rhotot(1:nfft)+xccc3d(1:nfft) 233 234 ! When possible, it is better to use the exact (analytical) gradient of nhat 235 ! In that case, we substract nhat from the density now 236 ! If not possible, the gradient will be computed in Fourier space (see call to xcden) 237 use_exact_nhat_gradient=(nhatdim==1.and.nhatgrdim==1.and.usexcnhat==1.and.intxc==0) 238 if (use_exact_nhat_gradient) rhotot(1:nfft)=rhotot(1:nfft)-nhat(1:nfft,1) 239 240 ! Loop on unshifted or shifted grids 241 do ishift=0,intxc 242 243 ! Set up density on unshifted or shifted grid (will be in rhonow(:,:,1)), 244 ! as well as the gradient of the density, also on the unshifted 245 ! or shifted grid (will be in rhonow(:,:,2:4)), if needed. 246 ABI_MALLOC(rhonow,(nfft,1,4)) 247 nspden1=1 ; cplex=1 ; ngrad=2 ; qphon(:)=zero 248 call xcden(cplex,gprimd,ishift,mpi_enreg,nfft,ngfft,ngrad,nspden1,qphon,rhotot,rhonow) 249 250 ! PAW: if they were previously substracted, re-add nhat and its "exact" gradient 251 if (use_exact_nhat_gradient) then 252 rhonow(1:nfft,1,1)=rhonow(1:nfft,1,1)+nhat(1:nfft,1) 253 do ii=1,3 254 rhonow(1:nfft,1,ii+1)=rhonow(1:nfft,1,ii+1)+nhatgr(1:nfft,1,ii) 255 end do 256 end if 257 258 ! Make the density positive everywhere (but do not care about gradients) 259 nspden1=1 ; iwarn=1 ; option=1 260 call mkdenpos(iwarn,nfft,nspden1,option,rhonow(:,1,1),xc_denpos) 261 262 ! Accumulate integral of |Grad_rho|/Rho (to be used for TB09 XC) 263 do ipts=1,nfft 264 rho=rhonow(ipts,1,1) 265 if (abs(rho)>tol10) then 266 grho2=rhonow(ipts,1,2)**2+rhonow(ipts,1,3)**2+rhonow(ipts,1,4)**2 267 grho_over_rho_pw=grho_over_rho_pw+sqrt(grho2)/rho 268 end if 269 end do 270 271 ABI_FREE(rhonow) 272 273 ! End loop on unshifted or shifted grids 274 end do 275 276 ! Normalize the integral 277 nfftot=ngfft(1)*ngfft(2)*ngfft(3) 278 grho_over_rho_pw=grho_over_rho_pw*ucvol/dble(nfftot)/dble(intxc+1) 279 280 ! Reduction in case of FFT distribution 281 if (mpi_enreg%nproc_fft>1) then 282 call xmpi_sum(grho_over_rho_pw,mpi_enreg%comm_fft,ierr) 283 end if 284 285 ABI_FREE(rhotot) 286 287 end if ! compute_pw 288 289 ! ======================================================================== 290 ! ========== PAW on-site contributions to Int[|Grad(Rho)|/Rho] =========== 291 292 if (usepaw==1.and.compute_paw) then 293 294 ! Set up parallelism over atoms 295 my_natom=mpi_enreg%my_natom 296 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 297 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 298 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 299 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom) 300 301 usecore=1 ; nzlmopt=-1 302 npts=pawang%angl_size 303 304 ! Compute Ylm gradients 305 ! dYlm/dr_i = { dYlm/dr_i^hat - r_i^hat * Sum_j[dYlm/dr_j^hat r_j^hat] } * (1/r) 306 ABI_MALLOC(dylmdr,(3,npts,pawang%ylm_size)) 307 do ilm=1,pawang%ylm_size 308 do ipts=1,npts 309 factor=sum(pawang%ylmrgr(1:3,ilm,ipts)*pawang%anginit(1:3,ipts)) 310 dylmdr(1:3,ipts,ilm)=pawang%ylmrgr(1:3,ilm,ipts)-factor*pawang%anginit(1:3,ipts) 311 end do 312 end do 313 314 ! Loop on atom sites 315 do iatom=1,my_natom 316 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 317 318 itypat=pawrhoij(iatom)%itypat 319 mesh_size=pawtab(itypat)%mesh_size 320 lm_size=pawtab(itypat)%lcut_size**2 321 lm_size_eff=min(lm_size,pawang%ylm_size) 322 ABI_MALLOC(ff,(mesh_size)) 323 324 ! Compute on-site densities 325 ABI_MALLOC(rho1 ,(mesh_size,lm_size,1)) 326 ABI_MALLOC(trho1,(mesh_size,lm_size,1)) 327 ABI_MALLOC(nhat1,(mesh_size,lm_size,usexcnhat)) 328 ABI_MALLOC(lmselect,(lm_size)) 329 ABI_MALLOC(lmselect_tmp,(lm_size)) 330 lmselect_tmp(:)=.true. 331 cplex=1 ; nspden1=1 ; opt_compch=0 ; opt_dens=1-usexcnhat ; pawprtvol=0 332 call pawdensities(compch_dum,cplex,iatom_tot,lmselect_tmp,lmselect,& 333 & lm_size,nhat1,nspden1,nzlmopt,opt_compch,opt_dens,-1,0,pawang,pawprtvol,& 334 & pawrad(itypat),pawrhoij(iatom),pawtab(itypat),rho1,trho1) 335 ABI_FREE(lmselect_tmp) 336 337 ! Allocation of temporary densities 338 ABI_MALLOC(rhotot,(mesh_size)) 339 ABI_MALLOC(grhotot,(mesh_size,3)) 340 ABI_MALLOC(rhoxc,(mesh_size,lm_size)) 341 ABI_MALLOC(drhoxc,(mesh_size)) 342 ABI_MALLOC(dcorexc,(mesh_size)) 343 344 ! Loop over AE and PS contributions 345 do iloop=1,2 346 if (iloop==1) then 347 rhoxc(:,:)=rho1(:,:,1) 348 corexc => pawtab(itypat)%coredens(:) 349 usenhat=0 350 else 351 rhoxc(:,:)=trho1(:,:,1) 352 corexc => pawtab(itypat)%tcoredens(:,1) 353 usenhat=usexcnhat 354 end if 355 356 ! Add hat density if needed 357 if (usenhat>0) rhoxc(:,:)=rhoxc(:,:)+nhat1(:,:,1) 358 359 ! Need derivative of core density 360 if (usecore==1) then 361 call nderiv_gen(dcorexc,corexc,pawrad(itypat)) 362 end if 363 364 ! Do loop on the angular part (theta,phi) 365 do ipts=1,npts 366 367 ! Compute the density and its gradient for this (theta,phi) 368 rhotot(:)=zero ; grhotot(:,:)=zero 369 do ilm=1,lm_size_eff 370 if (lmselect(ilm)) then 371 !Density 372 rhotot(:)=rhotot(:)+rhoxc(:,ilm)*pawang%ylmr(ilm,ipts) 373 !Gradient 374 ff(1:mesh_size)=rhoxc(1:mesh_size,ilm) 375 call nderiv_gen(drhoxc,ff,pawrad(itypat)) 376 ff(2:mesh_size)=ff(2:mesh_size)/pawrad(itypat)%rad(2:mesh_size) 377 call pawrad_deducer0(ff,mesh_size,pawrad(itypat)) 378 do ii=1,3 379 grhotot(1:mesh_size,ii)=grhotot(1:mesh_size,ii) & 380 & +drhoxc(1:mesh_size)*pawang%ylmr(ilm,ipts)*pawang%anginit(ii,ipts) & 381 & +ff(1:mesh_size)*dylmdr(ii,ipts,ilm) 382 end do 383 end if 384 end do 385 if (usecore==1) then 386 rhotot(:)=rhotot(:)+corexc(:) 387 do ii=1,3 388 grhotot(:,ii)=grhotot(:,ii)+dcorexc(:)*pawang%anginit(ii,ipts) 389 end do 390 end if 391 392 ! Make the density positive everywhere (but do not care about gradients) 393 nspden1=1 ; iwarn=1 394 call mkdenpos(iwarn,mesh_size,nspden1,0,rhotot,xc_denpos) 395 396 ! Accumulate integral of |Grad_rho|/Rho 397 do ir=1,mesh_size 398 rho=rhotot(ir) 399 if (abs(rho)>tol10) then 400 grho2=grhotot(ir,1)**2+grhotot(ir,2)**2+grhotot(ir,3)**2 401 ff(ir)=sqrt(grho2)/rho 402 end if 403 end do 404 ff(1:mesh_size)=ff(1:mesh_size)*pawrad(itypat)%rad(1:mesh_size)**2 405 call simp_gen(sumg,ff,pawrad(itypat)) 406 407 ! Add the contribution of the atom 408 if (iloop==1) then 409 grho_over_rho_paw=grho_over_rho_paw+sumg*four_pi*pawang%angwgth(ipts) 410 else 411 grho_over_rho_paw=grho_over_rho_paw-sumg*four_pi*pawang%angwgth(ipts) 412 end if 413 414 end do ! npts 415 end do ! AE/PS loop 416 417 ! Deallocate temporary memory 418 ABI_FREE(ff) 419 ABI_FREE(rho1) 420 ABI_FREE(trho1) 421 ABI_FREE(nhat1) 422 ABI_FREE(lmselect) 423 ABI_FREE(rhotot) 424 ABI_FREE(grhotot) 425 ABI_FREE(rhoxc) 426 ABI_FREE(drhoxc) 427 ABI_FREE(dcorexc) 428 429 end do ! atom sites 430 431 ! Deallocate temporary memory 432 ABI_FREE(dylmdr) 433 434 ! Reduction in case of parallelism 435 if (paral_atom) then 436 call xmpi_sum(grho_over_rho_paw,my_comm_atom,ierr) 437 end if 438 439 ! Destroy atom table used for parallelism 440 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 441 442 end if ! PAW and compute_paw 443 444 ! ======================================================================== 445 ! ============ Assemble contributions and update mBJ C value ============ 446 447 grho_over_rho = (grho_over_rho_pw + grho_over_rho_paw)/ucvol 448 xc_c_tb09 = -0.012_dp + 1.023_dp * sqrt(grho_over_rho) 449 if (xc_c_tb09<one) xc_c_tb09=one 450 if (xc_c_tb09>two) xc_c_tb09=two 451 452 write(msg,'(2a,f10.5,a)') ch10, & 453 & ' In the mGGA functional TB09, c (computed value) = ',xc_c_tb09,ch10 454 call wrtout(std_out,msg,'COLL') 455 456 if (present(xc_funcs)) then 457 call libxc_functionals_set_c_tb09(xc_c_tb09,xc_functionals=xc_funcs) 458 else 459 call libxc_functionals_set_c_tb09(xc_c_tb09) 460 end if 461 462 !call timab(xx,2,tsec) 463 464 end subroutine xc_tb09_update_c