TABLE OF CONTENTS
ABINIT/m_paw_init [ Modules ]
NAME
m_paw_init
FUNCTION
This module contains routines related tp PAW calculations initialization.
COPYRIGHT
Copyright (C) 2018-2022 ABINIT group (FJ, 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 MODULE m_paw_init 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 use m_splines 28 use m_dtset 29 30 use m_time, only : timab 31 use m_pawpsp, only : pawpsp_nl 32 use m_paw_atom,only : atompaw_shpfun 33 use m_pawang, only : pawang_type, pawang_init, pawang_free 34 use m_pawrad, only : pawrad_type, simp_gen, nderiv_gen, poisson, pawrad_deducer0 35 use m_pawtab, only : pawtab_type 36 use m_pawxc, only : pawxc_get_usekden,pawxc_get_uselaplacian,pawxc_get_xclevel 37 use m_paw_numeric, only : paw_derfc 38 39 implicit none 40 41 private 42 43 !public procedures. 44 public :: pawinit ! Initialize some tabulated data for PAW calculations 45 public :: paw_gencond ! Test whether we have to call pawinit to regenerate tabulated data. 46 47 CONTAINS !========================================================================================
m_paw_init/paw_gencond [ Functions ]
[ Top ] [ m_paw_init ] [ Functions ]
NAME
paw_gencond
FUNCTION
This routine tests whether we have to call pawinit in one of the optdriver routines since important values have changed wrt to the previous dataset. The function uses an internal array to store of the previous values Usage example: call paw_gencond(Dtset,gnt_option,"test",call_pawinit) if (psp_gencond == 1 .or. call_pawinit) then call pawinit(...) call paw_gencond(Dtset, gnt_option, "save", call_pawinit) end if where psp_gencond is the value returned by pspini. INPUT Dtset<type(dataset_type)>=all input variables for this dataset gnt_option=flag activated if pawang%gntselect and pawang%realgnt have to be allocated also determine the size of these pointers mode= "test" to test if pawinit must be called "save" to update the internal variables. "reset" to reset the internal variables
OUTPUT
call_pawinit=True if pawinit must be called. Meaninfull only if mode=="test"
SIDE EFFECTS
mode=="save" updates the internal variables. "reset" reset the internal variables to -1
SOURCE
735 subroutine paw_gencond(Dtset,gnt_option,mode,call_pawinit) 736 737 !Arguments ------------------------------------ 738 integer,intent(in) :: gnt_option 739 logical,intent(out) :: call_pawinit 740 character(len=*),intent(in) :: mode 741 type(dataset_type),intent(in) :: Dtset 742 743 !Local variables------------------------------- 744 !scalars 745 integer,save :: gencond(10)=(/-1,-1,-1,-1,-1,-1,-1,-1,-1,-1/) 746 747 ! ********************************************************************* 748 749 call_pawinit = .False. 750 select case (mode) 751 case ("test") 752 753 if (gencond(1)/=Dtset%pawlcutd .or.gencond(2) /=Dtset%pawlmix .or.& 754 & gencond(3)/=Dtset%pawnphi .or.gencond(4) /=Dtset%pawntheta.or.& 755 & gencond(5)/=Dtset%pawspnorb .or.gencond(6) /=Dtset%pawxcdev .or.& 756 & gencond(7)/=Dtset%nsym .or.gencond(8) /=gnt_option .or.& 757 & gencond(9)/=Dtset%usepotzero.or.gencond(10)/=Dtset%usekden) call_pawinit = .True. 758 759 case ("save") 760 ! Update internal values 761 gencond(1)=Dtset%pawlcutd ; gencond(2) =Dtset%pawlmix 762 gencond(3)=Dtset%pawnphi ; gencond(4) =Dtset%pawntheta 763 gencond(5)=Dtset%pawspnorb ; gencond(6) =Dtset%pawxcdev 764 gencond(7)=Dtset%nsym ; gencond(8) =gnt_option 765 gencond(9)=Dtset%usepotzero; gencond(10)=Dtset%usekden 766 767 case ("reset") 768 gencond = -1 769 770 case default 771 ABI_BUG("Wrong value for mode: "//trim(mode)) 772 end select 773 774 end subroutine paw_gencond
m_paw_init/pawinit [ Functions ]
[ Top ] [ m_paw_init ] [ Functions ]
NAME
pawinit
FUNCTION
Initialize some starting values of several arrays used in PAW calculations. 1-Initialize data related to angular mesh 2-Tabulate normalized shape function g(r) 3-Compute indklmn indexes giving some l,m,n,lp,mp,np info from klmn=[(l,m,n),(lp,mp,np)] 4-Compute various factors/sizes (depending on (l,m,n)) 5-Compute $q_ijL=\displaystyle \int_{0}^{r_c}{(\phi_i\phi_j-\widetilde{\phi_i}\widetilde{\phi_j}) r^l\,dr} Gaunt(l_i m_i,l_j m_j,l m))$ $S_ij=\displaystyle \sqrt{4 \pi} q_ij0$ 6-Compute $e_ijkl= vh1_ijkl - Vhatijkl - Bijkl - Cijkl$ With: $vh1_ijkl =\sum_{L,m} {vh1*Gaunt(i,j,Lm)*Gaunt(k,l,Lm)}$ $Vhat_ijkl=\sum_{L,m} {vhatijL*Gaunt(i,j,Lm)*q_klL}$ $B_ijkl =\sum_{L,m} {vhatijL*Gaunt(k,l,Lm)*q_ijL}$ $C_ijkl =\sum_{L,m} {intvhatL*q_ijL*q_klL}$ and: vh1 according to eq. (A17) in Holzwarth et al., PRB 55, 2005 (1997) [[cite:Holzwarth1997]] 7-Compute Ex-correlation energy for the core density
COPYRIGHT
Copyright (C) 1998-2022 ABINIT group (FJ, 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
effmass_free=effective mass for electrons (1. in common case) gnt_option=flag activated if pawang%gntselect and pawang%realgnt have to be allocated also determine the size of these pointers gsqcut_shp=effective cut-off to determine shape functions in reciprocal space hyb_range_fock=range coefficient for screened hybrid XC functionals ixc=choice of exchange-correlation functional lcutdens=max. l for densities/potentials moments computations lmix=max. l for which spherical terms will be mixed durinf SCF cycle mpsang=1+maximum angular momentum nphi="phi" dimension of paw angular mesh nsym=Number of symmetry elements in space group ntheta="theta" dimension of paw angular mesh pawang <type(pawang_type)>=paw angular mesh and related data %lmax=Maximum value of angular momentum l+1 %gntselect((2*l_max-1)**2,l_max**2,l_max**2)= selection rules for Gaunt coefficients pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data: %mesh_size=Dimension of radial mesh %rad(mesh_size)=The coordinates of all the points of the radial mesh %radfact(mesh_size)=Factor used to compute radial integrals pawspnorb=flag: 1 if spin-orbit coupling is activated pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data: %basis_size=Number of elements for the PAW nl basis %l_size=Maximum value of l+1 leading to non zero Gaunt coeffs %lmn_size=Number of (l,m,n) elements for the PAW basis %lmn2_size=lmn_size*(lmn_size+1)/2 %dltij(lmn2_size)=factors used to compute sums over (ilmn,jlmn) %phi(mesh_size,basis_size)=PAW all electron wavefunctions %rshp=shape function radius (radius for compensation charge) %shape_type=Radial shape function type %shape_alpha=Alpha parameters in Bessel shape function %shape_lambda=Lambda parameter in gaussian shape function %shape_q=Q parameters in Bessel shape function %shape_sigma=Sigma parameter in gaussian shape function %tphi(mesh_size,basis_size)=PAW atomic pseudowavefunctions pawxcdev=Choice of XC development (0=no dev. (use of angular mesh) ; 1=dev. on moments) usekden= 1 is kinetic energy density has to be computed, 0 otherwise
OUTPUT
pawang %gntselect(l_size_max**2,l_max**2*(l_max**2+1)/2)=selection rules for Gaunt coefficients %l_max=maximum value of angular momentum l+1 %l_size_max=maximum value of angular momentum l_size=2*l_max-1 %nsym=number of symmetry elements in space group %ngnt=number of non-zero Gaunt coefficients %realgnt(pawang%ngnt)=non-zero real Gaunt coefficients === only if pawxcdev==1 == %anginit(3,angl_size)=for each point of the angular mesh, gives the coordinates of the corresponding point on an unitary sphere %angl_size=dimension of paw angular mesh (angl_size=ntheta*nphi) %angwgth(angl_size)=for each point of the angular mesh, gives the weight of the corresponding point on an unitary sphere %ntheta, nphi=dimensions of paw angular mesh %ylmr(l_size_max**2,angl_size)=real Ylm calculated in real space %ylmrgr(1:3,l_size_max**2,angl_size)=first gradients of real Ylm calculated in real space === only if pawspnorb==1 == %ls_ylm(2,l_max**2,l_max**2,4)=LS operator in the real spherical harmonics basis %use_ls_ylm=flag activated if ls_ylm is allocated %usespnorb=flag activated for spin-orbit coupling pawtab(ntypat) <type(pawtab_type)>=paw tabulated data read at start: %lcut_size_=max. value of l+1 leading to non zero Gaunt coeffs modified by lcutdens %lmnmix_sz=number of (lmn,lmn_prime) verifying l<=lmix and l_prime<=lmix %mqgrid_shp=number of points in reciprocal space for shape function %indklmn(8,lmn2_size)=array giving klm, kln, abs(il-jl) and (il+jl), ilmn and jlmn for each klmn=(ilmn,jlmn) %dshpfunc(mesh_size,l_size,4)=derivatives of shape function (used only for numerical shape functions) %eijkl(lmn2_size,lmn2_size)=part of the Dij that depends only from the projected occupation coeffs %exccore=Exchange-correlation energy for the core density %gnorm(l_size)=normalization factor of radial shape function %phiphj(:,:)=useful product Phi(:,i)*Phi(:,j) %qgrid_shp(mqgrid_shp)=points in reciprocal space for shape function %qijl(l_size**2,lmn2_size)=moments of the difference charge density between AE and PS partial wave %rad_for_spline(mesh_size)=radial grid used for spline (copy of pawrad%rad) %shapefunc(mesh_size,l_size)=normalized radial shape function %shapefncg(mqgrid_shp,l_size)=normalized radial shape function in reciprocal space %sij(lmn2_size)=nonlocal part of the overlap operator %tphitphj(:,:)=useful product tPhi(:,i)*tPhi(:,j)
SOURCE
165 subroutine pawinit(effmass_free,gnt_option,gsqcut_eff,hyb_range_fock,lcutdens,lmix,mpsang,& 166 & nphi,nsym,ntheta,pawang,pawrad,pawspnorb,pawtab,pawxcdev,ixc,usepotzero) 167 168 !Arguments --------------------------------------------- 169 !scalars 170 integer,intent(in) :: gnt_option,ixc,lcutdens,lmix,mpsang,nphi,nsym,ntheta 171 integer,intent(in) :: pawspnorb,pawxcdev,usepotzero 172 real(dp),intent(in) :: effmass_free,gsqcut_eff,hyb_range_fock 173 type(pawang_type),intent(inout) :: pawang 174 !arrays 175 type(pawrad_type),intent(in) :: pawrad(:) 176 type(pawtab_type),target,intent(inout) :: pawtab(:) 177 178 !Local variables ------------------------------ 179 !scalars 180 integer,parameter :: mqgrid_shp_default=300 181 integer :: basis_size,i0lm,i0ln,ij_size,il,ilm,ilmn,iln,iloop,iq,isel,isel1 182 integer :: itypat,j0lm,j0lmn,j0ln,jl,jlm,jlmn,jln,klm,klm1 183 integer :: klmn,klmn1,kln,kln1,l_size,ll,lm0,lmax,lmax1,lmin,lmin1,lmn2_size 184 integer :: lmn_size,lmnmix,mesh_size,meshsz,mm,nabgnt_option,ngrad2_ylm,ntypat,pw_mesh_size 185 integer :: usexcnhat,use_angular_grid,use_ls_ylm,use_ylm,usekden 186 real(dp) :: dq,gnrm,intg,ql,ql1,rg,rg1,vh1,yp1,ypn 187 character(len=500) :: message 188 !arrays 189 integer,allocatable :: indl(:,:),klm_diag(:),kmix_tmp(:) 190 integer, ABI_CONTIGUOUS pointer :: indlmn(:,:) 191 real(dp) :: tsec(2) 192 real(dp),allocatable :: der(:),ff(:),gg(:),hh(:),indklmn_(:,:),intvhatl(:) 193 real(dp),allocatable :: rad(:),rgl(:,:),vhatijl(:,:),vhatl(:),work(:) 194 real(dp),pointer :: eijkl(:,:) 195 196 !************************************************************************ 197 198 DBG_ENTER("COLL") 199 200 call timab(553,1,tsec) 201 202 ntypat=size(pawtab) 203 if (size(pawrad)/=ntypat) then 204 ABI_BUG('pawrad and pawtab should have the same size!') 205 end if 206 207 !Immediately set the value of usepotzero 208 !it will be used later on in this subroutine 209 pawtab%usepotzero=usepotzero 210 211 usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat) 212 213 !================================================== 214 !1- INITIALIZE DATA RELATED TO ANGULAR MESH 215 !* ANGULAR GRID 216 !* REAL SPHERICAL HARMONICS 217 !* REAL GAUNT COEFFICIENTS 218 219 usekden=pawxc_get_usekden(ixc) 220 nabgnt_option=0;if (usekden>0) nabgnt_option=1 ! If kin. ene. density is used, need nabla Gaunt coeffs 221 use_angular_grid=0;if (pawxcdev==0) use_angular_grid=1 222 use_ylm=0;if (pawxcdev==0) use_ylm=1 223 use_ls_ylm=0;if (pawspnorb>0) use_ls_ylm=1 224 ngrad2_ylm=0;if (pawxc_get_xclevel(ixc)>=2) ngrad2_ylm=1 225 if (pawxc_get_uselaplacian(ixc)>0) ngrad2_ylm=2 226 call pawang_free(pawang) 227 call pawang_init(pawang,gnt_option,nabgnt_option,mpsang-1,nphi,ntheta,nsym,ngrad2_ylm,& 228 & use_angular_grid,use_ylm,use_ls_ylm) 229 230 !******************* 231 !Loop on atom types 232 !******************* 233 do itypat=1,ntypat 234 mesh_size=pawtab(itypat)%mesh_size 235 l_size=pawtab(itypat)%l_size 236 lmn_size=pawtab(itypat)%lmn_size 237 lmn2_size=pawtab(itypat)%lmn2_size 238 basis_size=pawtab(itypat)%basis_size 239 ij_size=pawtab(itypat)%ij_size 240 indlmn => pawtab(itypat)%indlmn(:,:) 241 ABI_MALLOC(indklmn_,(8,lmn2_size)) 242 ABI_MALLOC(klm_diag,(lmn2_size)) 243 ABI_MALLOC(ff,(mesh_size)) 244 ABI_MALLOC(gg,(mesh_size)) 245 ABI_MALLOC(hh,(mesh_size)) 246 ABI_MALLOC(rad,(mesh_size)) 247 rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size) 248 249 if (pawtab(itypat)%usexcnhat/=usexcnhat) then 250 write(message, '(7a)' )& 251 & 'You cannot simultaneously use atomic data with different',ch10,& 252 & 'formulation of XC [using compensation charge in XC or not] !',ch10,& 253 & 'Action: change at least one of your atomic data (psp) file',ch10,& 254 & ' or use usexcnhat keyword in input file.' 255 ABI_ERROR(message) 256 end if 257 258 ! ================================================== 259 ! 2- TABULATE SHAPE FUNCTION 260 261 ! Allocated shape function 262 if (pawtab(itypat)%shape_type/=-1) then 263 if (allocated(pawtab(itypat)%shapefunc)) then 264 ABI_FREE(pawtab(itypat)%shapefunc) 265 end if 266 ABI_MALLOC(pawtab(itypat)%shapefunc,(mesh_size,l_size)) 267 else if (.not.allocated(pawtab(itypat)%shapefunc)) then 268 message='shapefunc should be allocated with shape_type=-1' 269 ABI_ERROR(message) 270 end if 271 if (allocated(pawtab(itypat)%gnorm)) then 272 ABI_FREE(pawtab(itypat)%gnorm) 273 end if 274 ABI_MALLOC(pawtab(itypat)%gnorm,(l_size)) 275 276 ! Compute shape function 277 do il=1,l_size 278 ll=il-1 279 call atompaw_shpfun(ll,pawrad(itypat),gnrm,pawtab(itypat),ff) 280 pawtab(itypat)%shapefunc(1:mesh_size,il)=ff(1:mesh_size) 281 pawtab(itypat)%gnorm(il)=gnrm 282 end do 283 ! In case of numerical shape function, compute some derivatives 284 if (pawtab(itypat)%shape_type==-1) then 285 if (allocated(pawtab(itypat)%dshpfunc)) then 286 ABI_FREE(pawtab(itypat)%dshpfunc) 287 end if 288 ABI_MALLOC(pawtab(itypat)%dshpfunc,(mesh_size,l_size,4)) 289 ABI_MALLOC(work,(mesh_size)) 290 do il=1,l_size 291 call nderiv_gen(pawtab(itypat)%dshpfunc(:,il,1),pawtab(itypat)%shapefunc(:,il),pawrad(itypat)) 292 yp1=pawtab(itypat)%dshpfunc(1,il,1);ypn=pawtab(itypat)%dshpfunc(mesh_size,il,1) 293 call spline(rad,pawtab(itypat)%shapefunc(:,il),mesh_size,yp1,ypn,pawtab(itypat)%dshpfunc(:,il,2)) 294 yp1=pawtab(itypat)%dshpfunc(1,il,2);ypn=pawtab(itypat)%dshpfunc(mesh_size,il,2) 295 call spline(rad,pawtab(itypat)%dshpfunc(:,il,1),mesh_size,yp1,ypn,pawtab(itypat)%dshpfunc(:,il,3)) 296 yp1=pawtab(itypat)%dshpfunc(1,il,3);ypn=pawtab(itypat)%dshpfunc(mesh_size,il,3) 297 call spline(rad,pawtab(itypat)%dshpfunc(:,il,2),mesh_size,yp1,ypn,pawtab(itypat)%dshpfunc(:,il,4)) 298 end do 299 ABI_FREE(work) 300 end if 301 302 ! In some cases, has to store radial mesh for shape function in pawtab variable 303 if (pawtab(itypat)%shape_type==-1) then 304 if (allocated(pawtab(itypat)%rad_for_spline)) then 305 ABI_FREE(pawtab(itypat)%rad_for_spline) 306 end if 307 ABI_MALLOC(pawtab(itypat)%rad_for_spline,(mesh_size)) 308 pawtab(itypat)%rad_for_spline(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size) 309 end if 310 311 ! In some cases, has to store shape function in reciprocal space 312 if (pawtab(itypat)%has_shapefncg>0) then 313 if (gsqcut_eff<tol8) then 314 message='Computation of shapefncg only possible when gsqcut>0!' 315 ABI_BUG(message) 316 end if 317 pawtab(itypat)%mqgrid_shp=mqgrid_shp_default 318 if (allocated(pawtab(itypat)%shapefncg)) then 319 ABI_FREE(pawtab(itypat)%shapefncg) 320 end if 321 if (allocated(pawtab(itypat)%qgrid_shp)) then 322 ABI_FREE(pawtab(itypat)%qgrid_shp) 323 end if 324 ABI_MALLOC(pawtab(itypat)%shapefncg,(pawtab(itypat)%mqgrid_shp,2,l_size)) 325 ABI_MALLOC(pawtab(itypat)%qgrid_shp,(pawtab(itypat)%mqgrid_shp)) 326 dq=1.1_dp*sqrt(gsqcut_eff)/dble(pawtab(itypat)%mqgrid_shp-1) 327 do iq=1,pawtab(itypat)%mqgrid_shp 328 pawtab(itypat)%qgrid_shp(iq)=dble(iq-1)*dq 329 end do 330 ABI_MALLOC(indl,(6,l_size)) 331 ABI_MALLOC(rgl,(mesh_size,il)) 332 do il=1,l_size 333 indl(:,il)=0;indl(1,il)=il-1;indl(5,il)=il 334 rgl(1:mesh_size,il)=rad(1:mesh_size)*pawtab(itypat)%shapefunc(1:mesh_size,il) 335 end do 336 call pawpsp_nl(pawtab(itypat)%shapefncg,indl,l_size,l_size,& 337 & pawtab(itypat)%mqgrid_shp,pawtab(itypat)%qgrid_shp,pawrad(itypat),rgl) 338 pawtab(itypat)%shapefncg=four_pi*pawtab(itypat)%shapefncg 339 ABI_FREE(indl) 340 ABI_FREE(rgl) 341 else 342 pawtab(itypat)%mqgrid_shp=0 343 end if 344 345 ! ================================================== 346 ! 3- COMPUTE indklmn INDEXES GIVING klm, kln, abs(il-jl) and (il+jl), ilmn and jlmn 347 ! for each klmn=(ilmn,jlmn) 348 349 if (allocated(pawtab(itypat)%indklmn)) then 350 ABI_FREE(pawtab(itypat)%indklmn) 351 end if 352 ABI_MALLOC(pawtab(itypat)%indklmn,(8,lmn2_size)) 353 354 klm_diag=0 355 do jlmn=1,lmn_size 356 jl= indlmn(1,jlmn);jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 357 j0lmn=jlmn*(jlmn-1)/2 358 j0lm =jlm *(jlm -1)/2 359 j0ln =jln *(jln -1)/2 360 do ilmn=1,jlmn 361 il= indlmn(1,ilmn);ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 362 klmn=j0lmn+ilmn 363 if (ilm<=jlm) then 364 indklmn_(1,klmn)=j0lm+ilm 365 else 366 i0lm=ilm*(ilm-1)/2 367 indklmn_(1,klmn)=i0lm+jlm 368 end if 369 if (iln<=jln) then 370 indklmn_(2,klmn)=j0ln+iln 371 else 372 i0ln=iln*(iln-1)/2 373 indklmn_(2,klmn)=i0ln+jln 374 end if 375 indklmn_(3,klmn)=min(abs(il-jl),lcutdens) 376 indklmn_(4,klmn)=min(il+jl,lcutdens) 377 indklmn_(5,klmn)=ilm 378 indklmn_(6,klmn)=jlm 379 indklmn_(7,klmn)=ilmn 380 indklmn_(8,klmn)=jlmn 381 pawtab(itypat)%indklmn(:,klmn)=indklmn_(:,klmn) 382 if (ilm==jlm) klm_diag(klmn)=1 383 end do 384 end do 385 386 ! ================================================== 387 ! 4- COMPUTE various FACTORS/SIZES (depending on (l,m,n)) 388 389 pawtab(itypat)%usespnorb=pawspnorb 390 pawtab(itypat)%lcut_size=min(l_size,lcutdens+1) 391 392 if (allocated(pawtab(itypat)%dltij)) then 393 ABI_FREE(pawtab(itypat)%dltij) 394 end if 395 ABI_MALLOC(pawtab(itypat)%dltij,(lmn2_size)) 396 pawtab(itypat)%dltij(:)=two 397 do ilmn=1,lmn_size 398 pawtab(itypat)%dltij(ilmn*(ilmn+1)/2)=one 399 end do 400 401 lmnmix=zero 402 ABI_MALLOC(kmix_tmp,(lmn2_size)) 403 do jlmn=1,lmn_size 404 jl=indlmn(1,jlmn) 405 if (jl<=lmix) then 406 j0lmn=jlmn*(jlmn-1)/2 407 do ilmn=1,jlmn 408 il=indlmn(1,ilmn) 409 if (il<=lmix) then 410 lmnmix=lmnmix+1 411 kmix_tmp(lmnmix)=j0lmn+ilmn 412 end if 413 end do 414 end if 415 end do 416 if (allocated(pawtab(itypat)%kmix)) then 417 ABI_FREE(pawtab(itypat)%kmix) 418 end if 419 ABI_MALLOC(pawtab(itypat)%kmix,(lmnmix)) 420 pawtab(itypat)%lmnmix_sz=lmnmix 421 pawtab(itypat)%kmix(1:lmnmix)=kmix_tmp(1:lmnmix) 422 ABI_FREE(kmix_tmp) 423 424 ! ================================================== 425 ! 5- STORE SOME USEFUL QUANTITIES FROM PARTIAL WAVES 426 427 if (allocated(pawtab(itypat)%phiphj)) then 428 ABI_FREE(pawtab(itypat)%phiphj) 429 end if 430 if (allocated(pawtab(itypat)%tphitphj)) then 431 ABI_FREE(pawtab(itypat)%tphitphj) 432 end if 433 ABI_MALLOC(pawtab(itypat)%phiphj,(mesh_size,ij_size)) 434 ABI_MALLOC(pawtab(itypat)%tphitphj,(mesh_size,ij_size)) 435 do jln=1,basis_size 436 j0ln=jln*(jln-1)/2 437 do iln=1,jln 438 kln=j0ln+iln 439 pawtab(itypat)%phiphj(1:mesh_size,kln)=pawtab(itypat)%phi(1:mesh_size,iln)& 440 & *pawtab(itypat)%phi(1:mesh_size,jln) 441 pawtab(itypat)%tphitphj(1:mesh_size,kln)=pawtab(itypat)%tphi(1:mesh_size,iln)& 442 & *pawtab(itypat)%tphi(1:mesh_size,jln) 443 end do 444 end do 445 446 if (usekden==1) then 447 pw_mesh_size=pawtab(itypat)%partialwave_mesh_size 448 if (allocated(pawtab(itypat)%nablaphi)) then 449 ABI_FREE(pawtab(itypat)%nablaphi) 450 end if 451 ABI_MALLOC(pawtab(itypat)%nablaphi,(pw_mesh_size,basis_size)) 452 if (allocated(pawtab(itypat)%tnablaphi)) then 453 ABI_FREE(pawtab(itypat)%tnablaphi) 454 end if 455 ABI_MALLOC(pawtab(itypat)%tnablaphi,(pw_mesh_size,basis_size)) 456 ABI_MALLOC(der,(pw_mesh_size)) 457 do iln=1,basis_size 458 call nderiv_gen(der,pawtab(itypat)%phi(1:pw_mesh_size,iln),pawrad(itypat)) 459 pawtab(itypat)%nablaphi(2:pw_mesh_size,iln)=der(2:pw_mesh_size) & 460 & -pawtab(itypat)%phi(2:pw_mesh_size,iln)/pawrad(itypat)%rad(2:pw_mesh_size) 461 call nderiv_gen(der,pawtab(itypat)%tphi(1:pw_mesh_size,iln),pawrad(itypat)) 462 pawtab(itypat)%tnablaphi(2:pw_mesh_size,iln)=der(2:pw_mesh_size) & 463 & -pawtab(itypat)%tphi(2:pw_mesh_size,iln)/pawrad(itypat)%rad(2:pw_mesh_size) 464 call pawrad_deducer0(pawtab(itypat)%nablaphi(1:pw_mesh_size,iln),pw_mesh_size,pawrad(itypat)) 465 call pawrad_deducer0(pawtab(itypat)%tnablaphi(1:pw_mesh_size,iln),pw_mesh_size,pawrad(itypat)) 466 end do 467 ABI_FREE(der) 468 pawtab(itypat)%has_nablaphi=2 469 end if 470 471 ! ================================================== 472 ! 6- COMPUTE Qijl TERMS AND Sij MATRIX 473 474 ! Store some usefull quantities 475 if (allocated(pawtab(itypat)%phiphj)) then 476 ABI_FREE(pawtab(itypat)%phiphj) 477 end if 478 if (allocated(pawtab(itypat)%tphitphj)) then 479 ABI_FREE(pawtab(itypat)%tphitphj) 480 end if 481 ABI_MALLOC(pawtab(itypat)%phiphj,(mesh_size,ij_size)) 482 ABI_MALLOC(pawtab(itypat)%tphitphj,(mesh_size,ij_size)) 483 do jln=1,basis_size 484 j0ln=jln*(jln-1)/2 485 do iln=1,jln 486 kln=j0ln+iln 487 pawtab(itypat)%phiphj (1:mesh_size,kln)=pawtab(itypat)%phi (1:mesh_size,iln)& 488 & *pawtab(itypat)%phi (1:mesh_size,jln) 489 pawtab(itypat)%tphitphj(1:mesh_size,kln)=pawtab(itypat)%tphi(1:mesh_size,iln)& 490 & *pawtab(itypat)%tphi(1:mesh_size,jln) 491 end do 492 end do 493 494 ! Compute q_ijL and S_ij=q_ij0 495 if (allocated(pawtab(itypat)%qijl)) then 496 ABI_FREE(pawtab(itypat)%qijl) 497 end if 498 if (allocated(pawtab(itypat)%sij)) then 499 ABI_FREE(pawtab(itypat)%sij) 500 end if 501 ABI_MALLOC(pawtab(itypat)%qijl,(l_size*l_size,lmn2_size)) 502 ABI_MALLOC(pawtab(itypat)%sij,(lmn2_size)) 503 pawtab(itypat)%qijl=zero 504 pawtab(itypat)%sij=zero 505 do klmn=1,lmn2_size 506 klm=indklmn_(1,klmn);kln=indklmn_(2,klmn) 507 lmin=indklmn_(3,klmn);lmax=indklmn_(4,klmn) 508 do ll=lmin,lmax,2 509 lm0=ll*ll+ll+1;ff(1)=zero 510 ff(2:mesh_size)=(pawtab(itypat)%phiphj (2:mesh_size,kln)& 511 & -pawtab(itypat)%tphitphj(2:mesh_size,kln))& 512 & *rad(2:mesh_size)**ll 513 call simp_gen(intg,ff,pawrad(itypat)) 514 do mm=-ll,ll 515 isel=pawang%gntselect(lm0+mm,klm) 516 if (isel>0) pawtab(itypat)%qijl(lm0+mm,klmn)=intg*pawang%realgnt(isel) 517 end do 518 end do 519 if (klm_diag(klmn)==1) pawtab(itypat)%sij(klmn)= & 520 & pawtab(itypat)%qijl(1,klmn)*sqrt(four_pi) 521 end do 522 523 ! ================================================== 524 ! 6- COMPUTE Eijkl TERMS (Hartree) 525 ! Compute eventually short-range screened version of Eijkl (Fock) 526 527 if (allocated(pawtab(itypat)%eijkl)) then 528 ABI_FREE(pawtab(itypat)%eijkl) 529 end if 530 ABI_MALLOC(pawtab(itypat)%eijkl,(lmn2_size,lmn2_size)) 531 if (abs(hyb_range_fock)>tol8) then 532 if (allocated(pawtab(itypat)%eijkl_sr)) then 533 ABI_FREE(pawtab(itypat)%eijkl_sr) 534 end if 535 ABI_MALLOC(pawtab(itypat)%eijkl_sr,(lmn2_size,lmn2_size)) 536 end if 537 538 ! First loop is for eijkl (Hartree) 539 ! 2nd loop is for eijkl_sr (short-range screened Fock exchange) 540 do iloop=1,2 541 if (iloop==2.and.abs(hyb_range_fock)<=tol8) cycle 542 if (iloop==1) eijkl => pawtab(itypat)%eijkl 543 if (iloop==2) eijkl => pawtab(itypat)%eijkl_sr 544 545 ! Compute: 546 ! vhatL(r) according to eq. (A14) in Holzwarth et al., PRB 55, 2005 (1997) [[cite:Holzwarth1997]] 547 ! intvhatL=$\int_{0}^{r_c}{vhatL(r) shapefunc_L(r) r^2\,dr}$ 548 ! vhatijL =$\int_{0}^{r_c}{vhatL(r) \tilde{\phi}_i \tilde{\phi}_j \,dr}$ 549 ! ----------------------------------------------------------------- 550 ABI_MALLOC(vhatl,(mesh_size)) 551 ABI_MALLOC(vhatijl,(lmn2_size,l_size)) 552 ABI_MALLOC(intvhatl,(l_size)) 553 intvhatl(:)=zero;vhatl(:)=zero;vhatijl(:,:)=zero 554 do il=1,l_size 555 vhatl(1)=zero;ff(1)=zero 556 ff(2:mesh_size)=pawtab(itypat)%shapefunc(2:mesh_size,il)*rad(2:mesh_size)**2 557 if (iloop==1) call poisson(ff,il-1,pawrad(itypat),vhatl) 558 if (iloop==2) call poisson(ff,il-1,pawrad(itypat),vhatl,screened_sr_separation=hyb_range_fock) 559 vhatl(2:mesh_size)=two*vhatl(2:mesh_size)/rad(2:mesh_size) 560 gg(1:mesh_size)=vhatl(1:mesh_size)*ff(1:mesh_size) 561 call simp_gen(intvhatl(il),gg,pawrad(itypat)) 562 do klmn=1,lmn2_size 563 kln=indklmn_(2,klmn) 564 hh(1:mesh_size)=vhatl(1:mesh_size)*pawtab(itypat)%tphitphj(1:mesh_size,kln) 565 call simp_gen(vhatijl(klmn,il),hh,pawrad(itypat)) 566 end do 567 end do 568 ABI_FREE(vhatl) 569 570 ! Compute: 571 ! eijkl=$ vh1_ijkl - Vhatijkl - Bijkl - Cijkl$ 572 ! With: 573 ! $vh1_ijkl =\sum_{L,m} {vh1*Gaunt(i,j,Lm)*Gaunt(k,l,Lm)}$ 574 ! $Vhat_ijkl=\sum_{L,m} {vhatijL*Gaunt(i,j,Lm)*q_klL}$ 575 ! $B_ijkl =\sum_{L,m} {vhatijL*Gaunt(k,l,Lm)*q_ijL}$ 576 ! $C_ijkl =\sum_{L,m} {intvhatL*q_ijL*q_klL}$ 577 ! and: 578 ! vh1 according to eq. (A17) in Holzwarth et al., PRB 55, 2005 (1997) [[cite:Holzwarth1997]] 579 ! Warning: compute only eijkl for (i,j)<=(k,l) 580 ! ----------------------------------------------------------------- 581 eijkl(:,:)=zero 582 meshsz=pawrad(itypat)%int_meshsz;if (mesh_size>meshsz) ff(meshsz+1:mesh_size)=zero 583 do klmn=1,lmn2_size 584 klm=indklmn_(1,klmn);kln=indklmn_(2,klmn) 585 lmin=indklmn_(3,klmn);lmax=indklmn_(4,klmn) 586 do ll=lmin,lmax,2 587 lm0=ll*ll+ll+1 588 ff(1:meshsz)=pawtab(itypat)%phiphj (1:meshsz,kln) 589 if (iloop==1) call poisson(ff,ll,pawrad(itypat),gg) 590 if (iloop==2) call poisson(ff,ll,pawrad(itypat),gg,screened_sr_separation=hyb_range_fock) 591 ff(1:meshsz)=pawtab(itypat)%tphitphj(1:meshsz,kln) 592 if (iloop==1) call poisson(ff,ll,pawrad(itypat),hh) 593 if (iloop==2) call poisson(ff,ll,pawrad(itypat),hh,screened_sr_separation=hyb_range_fock) 594 do klmn1=klmn,lmn2_size 595 klm1=indklmn_(1,klmn1);kln1=indklmn_(2,klmn1) 596 lmin1=indklmn_(3,klmn1);lmax1=indklmn_(4,klmn1) 597 vh1=zero 598 if ((ll.ge.lmin1).and.(ll.le.lmax1)) then 599 ff(1)=zero 600 ff(2:meshsz)=(pawtab(itypat)%phiphj (2:meshsz,kln1)*gg(2:meshsz)& 601 & -pawtab(itypat)%tphitphj(2:meshsz,kln1)*hh(2:meshsz))& 602 & *two/rad(2:meshsz) 603 call simp_gen(vh1,ff,pawrad(itypat)) 604 end if 605 do mm=-ll,ll 606 isel =pawang%gntselect(lm0+mm,klm) 607 isel1=pawang%gntselect(lm0+mm,klm1) 608 if (isel>0.and.isel1>0) then 609 rg =pawang%realgnt(isel) 610 rg1=pawang%realgnt(isel1) 611 ql =pawtab(itypat)%qijl(lm0+mm,klmn) 612 ql1=pawtab(itypat)%qijl(lm0+mm,klmn1) 613 eijkl(klmn,klmn1)=eijkl(klmn,klmn1)& 614 & +( vh1 *rg *rg1& ! vh1_ijkl 615 & - vhatijl(klmn ,ll+1)*rg *ql1& ! Vhat_ijkl 616 & - vhatijl(klmn1,ll+1)*rg1*ql & ! B_ijkl 617 & - intvhatl(ll+1) *ql *ql1& ! C_ijkl 618 & )*two_pi 619 end if 620 end do 621 end do 622 end do 623 end do 624 ABI_FREE(vhatijl) 625 ABI_FREE(intvhatl) 626 end do ! iloop 627 628 ! ================================================== 629 ! 7- COMPUTE gamma_ij TERMS 630 ! Corrections to get the background right 631 632 if (pawtab(itypat)%usepotzero==1) then 633 if (allocated(pawtab(itypat)%gammaij)) then 634 ABI_FREE(pawtab(itypat)%gammaij) 635 end if 636 ABI_MALLOC(pawtab(itypat)%gammaij,(lmn2_size)) 637 ABI_MALLOC(work,(mesh_size)) 638 do klmn=1,lmn2_size 639 if (klm_diag(klmn)==1) then 640 kln=indklmn_(2,klmn) 641 ff(1)=zero 642 ff(2:mesh_size)=pawtab(itypat)%phiphj(2:mesh_size,kln)-pawtab(itypat)%tphitphj(2:mesh_size,kln) 643 ! First, compute q_ij^00 644 call simp_gen(intg,ff,pawrad(itypat)) 645 ! Second, compute phi^2 - tphi^2 - 4pi*r^2 q_ij^00 g_0(r) 646 ff(2:mesh_size)= ff(2:mesh_size) - intg*pawtab(itypat)%shapefunc(2:mesh_size,1)*rad(2:mesh_size)**2 647 call poisson(ff,0,pawrad(itypat),work) 648 ! work is r*vh; should be then multiplied by r to prepare the integration in the sphere 649 work(1)=zero ; work(2:mesh_size)=work(2:mesh_size)*rad(2:mesh_size) 650 ! Third, average it over the sphere 651 call simp_gen(intg,work,pawrad(itypat)) 652 ! Finally, store it in pawtab%gammaij 653 pawtab(itypat)%gammaij(klmn)=intg*four_pi 654 else 655 pawtab(itypat)%gammaij(klmn)=zero 656 end if 657 end do 658 ABI_FREE(work) 659 end if 660 661 ! ================================================== 662 ! 8- TAKE into account a modified effective mass for the electrons 663 664 if (abs(effmass_free-one)>tol8) then 665 if (pawtab(itypat)%has_kij/=2) then 666 message='we need kij and has_kij/=2!' 667 ABI_BUG(message) 668 end if 669 if (allocated(pawtab(itypat)%dij0)) then 670 pawtab(itypat)%dij0(1:lmn2_size)=pawtab(itypat)%dij0(1:lmn2_size)-pawtab(itypat)%kij(1:lmn2_size) 671 end if 672 pawtab(itypat)%kij(1:lmn2_size)=pawtab(itypat)%kij(1:lmn2_size)/effmass_free 673 if (allocated(pawtab(itypat)%dij0)) then 674 pawtab(itypat)%dij0(1:lmn2_size)=pawtab(itypat)%dij0(1:lmn2_size)+pawtab(itypat)%kij(1:lmn2_size) 675 end if 676 end if 677 678 ! *********************** 679 ! End Loop on atom types 680 ! *********************** 681 ABI_FREE(ff) 682 ABI_FREE(gg) 683 ABI_FREE(hh) 684 ABI_FREE(indklmn_) 685 ABI_FREE(klm_diag) 686 ABI_FREE(rad) 687 end do 688 689 call timab(553,2,tsec) 690 691 DBG_EXIT("COLL") 692 693 end subroutine pawinit