TABLE OF CONTENTS
ABINIT/paw_gencond [ 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
PARENTS
bethe_salpeter,gstate,respfn,screening,sigma,wfk_analyze
CHILDREN
SOURCE
649 subroutine paw_gencond(Dtset,gnt_option,mode,call_pawinit) 650 651 use defs_basis 652 use m_errors 653 654 use defs_abitypes, only : dataset_type 655 656 !This section has been created automatically by the script Abilint (TD). 657 !Do not modify the following lines by hand. 658 #undef ABI_FUNC 659 #define ABI_FUNC 'paw_gencond' 660 !End of the abilint section 661 662 implicit none 663 664 !Arguments ------------------------------------ 665 integer,intent(in) :: gnt_option 666 logical,intent(out) :: call_pawinit 667 character(len=*),intent(in) :: mode 668 type(dataset_type),intent(in) :: Dtset 669 670 !Local variables------------------------------- 671 !scalars 672 integer,save :: gencond(9)=(/-1,-1,-1,-1,-1,-1,-1,-1,-1/) 673 674 ! ********************************************************************* 675 676 call_pawinit = .False. 677 select case (mode) 678 case ("test") 679 680 if (gencond(1)/=Dtset%pawlcutd .or.gencond(2)/=Dtset%pawlmix .or.& 681 & gencond(3)/=Dtset%pawnphi .or.gencond(4)/=Dtset%pawntheta.or.& 682 & gencond(5)/=Dtset%pawspnorb.or.gencond(6)/=Dtset%pawxcdev.or.& 683 & gencond(7)/=Dtset%nsym .or.gencond(8)/=gnt_option.or.& 684 & gencond(9)/=Dtset%usepotzero) call_pawinit = .True. 685 686 case ("save") 687 ! Update internal values 688 gencond(1)=Dtset%pawlcutd ; gencond(2)=Dtset%pawlmix 689 gencond(3)=Dtset%pawnphi ; gencond(4)=Dtset%pawntheta 690 gencond(5)=Dtset%pawspnorb; gencond(6)=Dtset%pawxcdev 691 gencond(7)=Dtset%nsym ; gencond(8)=gnt_option 692 gencond(9)=Dtset%usepotzero 693 694 case ("reset") 695 gencond = -1 696 697 case default 698 MSG_BUG("Wrong value for mode: "//trim(mode)) 699 end select 700 701 end subroutine paw_gencond
ABINIT/pawinit [ 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)
7-Compute Ex-correlation energy for the core density
COPYRIGHT
Copyright (C) 1998-2018 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
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
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)
xclevel=XC functional level (1=LDA, 2=GGA)
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
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)
PARENTS
bethe_salpeter,gstate,respfn,screening,sigma,wfk_analyze
CHILDREN
SOURCE
118 #if defined HAVE_CONFIG_H 119 #include "config.h" 120 #endif 121 122 #include "abi_common.h" 123 124 subroutine pawinit(gnt_option,gsqcut_eff,hyb_range_fock,lcutdens,lmix,mpsang,nphi,nsym,ntheta,& 125 & pawang,pawrad,pawspnorb,pawtab,pawxcdev,xclevel,usepotzero) 126 127 use defs_basis 128 use m_errors 129 use m_profiling_abi 130 use m_splines 131 132 use m_pawpsp, only : pawpsp_nl 133 use m_paw_atom,only : atompaw_shpfun 134 use m_pawang, only : pawang_type, pawang_init, pawang_free 135 use m_pawrad, only : pawrad_type, simp_gen, nderiv_gen, poisson 136 use m_pawtab, only : pawtab_type 137 use m_paw_numeric, only : paw_derfc 138 139 !This section has been created automatically by the script Abilint (TD). 140 !Do not modify the following lines by hand. 141 #undef ABI_FUNC 142 #define ABI_FUNC 'pawinit' 143 use interfaces_18_timing 144 !End of the abilint section 145 146 implicit none 147 148 !Arguments --------------------------------------------- 149 !scalars 150 integer,intent(in) :: gnt_option,lcutdens,lmix,mpsang,nphi,nsym,ntheta 151 integer,intent(in) :: pawspnorb,pawxcdev,xclevel,usepotzero 152 real(dp),intent(in) :: gsqcut_eff,hyb_range_fock 153 type(pawang_type),intent(inout) :: pawang 154 !arrays 155 type(pawrad_type),intent(in) :: pawrad(:) 156 type(pawtab_type),target,intent(inout) :: pawtab(:) 157 158 !Local variables ------------------------------ 159 !scalars 160 integer,parameter :: mqgrid_shp_default=300 161 integer :: basis_size,i0lm,i0ln,ij_size,il,ilm,ilmn,iln,iloop,iq,isel,isel1 162 integer :: itypat,j0lm,j0lmn,j0ln,jl,jlm,jlmn,jln,klm,klm1 163 integer :: klmn,klmn1,kln,kln1,l_size,ll,lm0,lmax,lmax1,lmin,lmin1,lmn2_size 164 integer :: lmn_size,lmnmix,mesh_size,meshsz,mm,ntypat,usexcnhat,use_ls_ylm,use_ylm 165 real(dp) :: dq,gnrm,intg,ql,ql1,rg,rg1,vh1,yp1,ypn 166 character(len=500) :: message 167 !arrays 168 integer,allocatable :: indl(:,:),klm_diag(:),kmix_tmp(:) 169 integer,ABI_CONTIGUOUS pointer :: indlmn(:,:) 170 real(dp) :: tsec(2) 171 real(dp),allocatable :: ff(:),gg(:),hh(:),indklmn_(:,:),intvhatl(:) 172 real(dp),allocatable :: rad(:),rgl(:,:),vhatijl(:,:),vhatl(:),work(:) 173 real(dp),pointer :: eijkl(:,:) 174 175 !************************************************************************ 176 177 DBG_ENTER("COLL") 178 179 call timab(553,1,tsec) 180 181 ntypat=size(pawtab) 182 if (size(pawrad)/=ntypat) then 183 MSG_BUG('pawrad and pawtab should have the same size!') 184 end if 185 186 ! Immediately set the value of usepotzero 187 ! it will be used later on in this subroutine 188 pawtab%usepotzero=usepotzero 189 190 !================================================== 191 !1- INITIALIZE DATA RELATED TO ANGULAR MESH 192 !* ANGULAR GRID 193 !* REAL SPHERICAL HARMONICS 194 !* REAL GAUNT COEFFICIENTS 195 196 use_ylm=0;if (pawxcdev==0) use_ylm=1 197 use_ls_ylm=0;if (pawspnorb>0) use_ls_ylm=1 198 call pawang_free(pawang) 199 call pawang_init(pawang,gnt_option,mpsang-1,nphi,nsym,ntheta,pawxcdev,use_ls_ylm,use_ylm,xclevel) 200 201 usexcnhat=maxval(pawtab(1:ntypat)%usexcnhat) 202 203 !******************* 204 !Loop on atom types 205 !******************* 206 do itypat=1,ntypat 207 mesh_size=pawtab(itypat)%mesh_size 208 l_size=pawtab(itypat)%l_size 209 lmn_size=pawtab(itypat)%lmn_size 210 lmn2_size=pawtab(itypat)%lmn2_size 211 basis_size=pawtab(itypat)%basis_size 212 ij_size=pawtab(itypat)%ij_size 213 indlmn => pawtab(itypat)%indlmn(:,:) 214 ABI_ALLOCATE(indklmn_,(8,lmn2_size)) 215 ABI_ALLOCATE(klm_diag,(lmn2_size)) 216 ABI_ALLOCATE(ff,(mesh_size)) 217 ABI_ALLOCATE(gg,(mesh_size)) 218 ABI_ALLOCATE(hh,(mesh_size)) 219 ABI_ALLOCATE(rad,(mesh_size)) 220 rad(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size) 221 222 if (pawtab(itypat)%usexcnhat/=usexcnhat) then 223 write(message, '(7a)' )& 224 & 'You cannot simultaneously use atomic data with different',ch10,& 225 & 'formulation of XC [using compensation charge in XC or not] !',ch10,& 226 & 'Action: change at least one of your atomic data (psp) file',ch10,& 227 & ' or use usexcnhat keyword in input file.' 228 MSG_ERROR(message) 229 end if 230 231 ! ================================================== 232 ! 2- TABULATE SHAPE FUNCTION 233 234 ! Allocated shape function 235 if (pawtab(itypat)%shape_type/=-1) then 236 if (allocated(pawtab(itypat)%shapefunc)) then 237 ABI_DEALLOCATE(pawtab(itypat)%shapefunc) 238 end if 239 ABI_ALLOCATE(pawtab(itypat)%shapefunc,(mesh_size,l_size)) 240 else if (.not.allocated(pawtab(itypat)%shapefunc)) then 241 message='shapefunc should be allocated with shape_type=-1' 242 MSG_ERROR(message) 243 end if 244 if (allocated(pawtab(itypat)%gnorm)) then 245 ABI_DEALLOCATE(pawtab(itypat)%gnorm) 246 end if 247 ABI_ALLOCATE(pawtab(itypat)%gnorm,(l_size)) 248 249 ! Compute shape function 250 do il=1,l_size 251 ll=il-1 252 call atompaw_shpfun(ll,pawrad(itypat),gnrm,pawtab(itypat),ff) 253 pawtab(itypat)%shapefunc(1:mesh_size,il)=ff(1:mesh_size) 254 pawtab(itypat)%gnorm(il)=gnrm 255 end do 256 ! In case of numerical shape function, compute some derivatives 257 if (pawtab(itypat)%shape_type==-1) then 258 if (allocated(pawtab(itypat)%dshpfunc)) then 259 ABI_DEALLOCATE(pawtab(itypat)%dshpfunc) 260 end if 261 ABI_ALLOCATE(pawtab(itypat)%dshpfunc,(mesh_size,l_size,4)) 262 ABI_ALLOCATE(work,(mesh_size)) 263 do il=1,l_size 264 call nderiv_gen(pawtab(itypat)%dshpfunc(:,il,1),pawtab(itypat)%shapefunc(:,il),pawrad(itypat)) 265 yp1=pawtab(itypat)%dshpfunc(1,il,1);ypn=pawtab(itypat)%dshpfunc(mesh_size,il,1) 266 call spline(rad,pawtab(itypat)%shapefunc(:,il),mesh_size,yp1,ypn,pawtab(itypat)%dshpfunc(:,il,2)) 267 yp1=pawtab(itypat)%dshpfunc(1,il,2);ypn=pawtab(itypat)%dshpfunc(mesh_size,il,2) 268 call spline(rad,pawtab(itypat)%dshpfunc(:,il,1),mesh_size,yp1,ypn,pawtab(itypat)%dshpfunc(:,il,3)) 269 yp1=pawtab(itypat)%dshpfunc(1,il,3);ypn=pawtab(itypat)%dshpfunc(mesh_size,il,3) 270 call spline(rad,pawtab(itypat)%dshpfunc(:,il,2),mesh_size,yp1,ypn,pawtab(itypat)%dshpfunc(:,il,4)) 271 end do 272 ABI_DEALLOCATE(work) 273 end if 274 275 ! In some cases, has to store radial mesh for shape function in pawtab variable 276 if (pawtab(itypat)%shape_type==-1) then 277 if (allocated(pawtab(itypat)%rad_for_spline)) then 278 ABI_DEALLOCATE(pawtab(itypat)%rad_for_spline) 279 end if 280 ABI_ALLOCATE(pawtab(itypat)%rad_for_spline,(mesh_size)) 281 pawtab(itypat)%rad_for_spline(1:mesh_size)=pawrad(itypat)%rad(1:mesh_size) 282 end if 283 284 ! In some cases, has to store shape function in reciprocal space 285 if (pawtab(itypat)%has_shapefncg>0) then 286 if (gsqcut_eff<tol8) then 287 message='Computation of shapefncg only possible when gsqcut>0!' 288 MSG_BUG(message) 289 end if 290 pawtab(itypat)%mqgrid_shp=mqgrid_shp_default 291 if (allocated(pawtab(itypat)%shapefncg)) then 292 ABI_DEALLOCATE(pawtab(itypat)%shapefncg) 293 end if 294 if (allocated(pawtab(itypat)%qgrid_shp)) then 295 ABI_DEALLOCATE(pawtab(itypat)%qgrid_shp) 296 end if 297 ABI_ALLOCATE(pawtab(itypat)%shapefncg,(pawtab(itypat)%mqgrid_shp,2,l_size)) 298 ABI_ALLOCATE(pawtab(itypat)%qgrid_shp,(pawtab(itypat)%mqgrid_shp)) 299 dq=1.1_dp*sqrt(gsqcut_eff)/dble(pawtab(itypat)%mqgrid_shp-1) 300 do iq=1,pawtab(itypat)%mqgrid_shp 301 pawtab(itypat)%qgrid_shp(iq)=dble(iq-1)*dq 302 end do 303 ABI_ALLOCATE(indl,(6,l_size)) 304 ABI_ALLOCATE(rgl,(mesh_size,il)) 305 do il=1,l_size 306 indl(:,il)=0;indl(1,il)=il-1;indl(5,il)=il 307 rgl(1:mesh_size,il)=rad(1:mesh_size)*pawtab(itypat)%shapefunc(1:mesh_size,il) 308 end do 309 call pawpsp_nl(pawtab(itypat)%shapefncg,indl,l_size,l_size,& 310 & pawtab(itypat)%mqgrid_shp,pawtab(itypat)%qgrid_shp,pawrad(itypat),rgl) 311 pawtab(itypat)%shapefncg=four_pi*pawtab(itypat)%shapefncg 312 ABI_DEALLOCATE(indl) 313 ABI_DEALLOCATE(rgl) 314 else 315 pawtab(itypat)%mqgrid_shp=0 316 end if 317 318 ! ================================================== 319 ! 3- COMPUTE indklmn INDEXES GIVING klm, kln, abs(il-jl) and (il+jl), ilmn and jlmn 320 ! for each klmn=(ilmn,jlmn) 321 322 if (allocated(pawtab(itypat)%indklmn)) then 323 ABI_DEALLOCATE(pawtab(itypat)%indklmn) 324 end if 325 ABI_ALLOCATE(pawtab(itypat)%indklmn,(8,lmn2_size)) 326 327 klm_diag=0 328 do jlmn=1,lmn_size 329 jl= indlmn(1,jlmn);jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn) 330 j0lmn=jlmn*(jlmn-1)/2 331 j0lm =jlm *(jlm -1)/2 332 j0ln =jln *(jln -1)/2 333 do ilmn=1,jlmn 334 il= indlmn(1,ilmn);ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn) 335 klmn=j0lmn+ilmn 336 if (ilm<=jlm) then 337 indklmn_(1,klmn)=j0lm+ilm 338 else 339 i0lm=ilm*(ilm-1)/2 340 indklmn_(1,klmn)=i0lm+jlm 341 end if 342 if (iln<=jln) then 343 indklmn_(2,klmn)=j0ln+iln 344 else 345 i0ln=iln*(iln-1)/2 346 indklmn_(2,klmn)=i0ln+jln 347 end if 348 indklmn_(3,klmn)=min(abs(il-jl),lcutdens) 349 indklmn_(4,klmn)=min(il+jl,lcutdens) 350 indklmn_(5,klmn)=ilm 351 indklmn_(6,klmn)=jlm 352 indklmn_(7,klmn)=ilmn 353 indklmn_(8,klmn)=jlmn 354 pawtab(itypat)%indklmn(:,klmn)=indklmn_(:,klmn) 355 if (ilm==jlm) klm_diag(klmn)=1 356 end do 357 end do 358 359 ! ================================================== 360 ! 4- COMPUTE various FACTORS/SIZES (depending on (l,m,n)) 361 362 pawtab(itypat)%lcut_size=min(l_size,lcutdens+1) 363 364 if (allocated(pawtab(itypat)%dltij)) then 365 ABI_DEALLOCATE(pawtab(itypat)%dltij) 366 end if 367 ABI_ALLOCATE(pawtab(itypat)%dltij,(lmn2_size)) 368 pawtab(itypat)%dltij(:)=two 369 do ilmn=1,lmn_size 370 pawtab(itypat)%dltij(ilmn*(ilmn+1)/2)=one 371 end do 372 373 lmnmix=zero 374 ABI_ALLOCATE(kmix_tmp,(lmn2_size)) 375 do jlmn=1,lmn_size 376 jl=indlmn(1,jlmn) 377 if (jl<=lmix) then 378 j0lmn=jlmn*(jlmn-1)/2 379 do ilmn=1,jlmn 380 il=indlmn(1,ilmn) 381 if (il<=lmix) then 382 lmnmix=lmnmix+1 383 kmix_tmp(lmnmix)=j0lmn+ilmn 384 end if 385 end do 386 end if 387 end do 388 if (allocated(pawtab(itypat)%kmix)) then 389 ABI_DEALLOCATE(pawtab(itypat)%kmix) 390 end if 391 ABI_ALLOCATE(pawtab(itypat)%kmix,(lmnmix)) 392 pawtab(itypat)%lmnmix_sz=lmnmix 393 pawtab(itypat)%kmix(1:lmnmix)=kmix_tmp(1:lmnmix) 394 ABI_DEALLOCATE(kmix_tmp) 395 396 ! ================================================== 397 ! 5- COMPUTE Qijl TERMS AND Sij MATRIX 398 399 ! Store some usefull quantities 400 if (allocated(pawtab(itypat)%phiphj)) then 401 ABI_DEALLOCATE(pawtab(itypat)%phiphj) 402 end if 403 if (allocated(pawtab(itypat)%tphitphj)) then 404 ABI_DEALLOCATE(pawtab(itypat)%tphitphj) 405 end if 406 ABI_ALLOCATE(pawtab(itypat)%phiphj,(mesh_size,ij_size)) 407 ABI_ALLOCATE(pawtab(itypat)%tphitphj,(mesh_size,ij_size)) 408 do jln=1,basis_size 409 j0ln=jln*(jln-1)/2 410 do iln=1,jln 411 kln=j0ln+iln 412 pawtab(itypat)%phiphj (1:mesh_size,kln)=pawtab(itypat)%phi (1:mesh_size,iln)& 413 & *pawtab(itypat)%phi (1:mesh_size,jln) 414 pawtab(itypat)%tphitphj(1:mesh_size,kln)=pawtab(itypat)%tphi(1:mesh_size,iln)& 415 & *pawtab(itypat)%tphi(1:mesh_size,jln) 416 end do 417 end do 418 419 ! Compute q_ijL and S_ij=q_ij0 420 if (allocated(pawtab(itypat)%qijl)) then 421 ABI_DEALLOCATE(pawtab(itypat)%qijl) 422 end if 423 if (allocated(pawtab(itypat)%sij)) then 424 ABI_DEALLOCATE(pawtab(itypat)%sij) 425 end if 426 ABI_ALLOCATE(pawtab(itypat)%qijl,(l_size*l_size,lmn2_size)) 427 ABI_ALLOCATE(pawtab(itypat)%sij,(lmn2_size)) 428 pawtab(itypat)%qijl=zero 429 pawtab(itypat)%sij=zero 430 do klmn=1,lmn2_size 431 klm=indklmn_(1,klmn);kln=indklmn_(2,klmn) 432 lmin=indklmn_(3,klmn);lmax=indklmn_(4,klmn) 433 do ll=lmin,lmax,2 434 lm0=ll*ll+ll+1;ff(1)=zero 435 ff(2:mesh_size)=(pawtab(itypat)%phiphj (2:mesh_size,kln)& 436 & -pawtab(itypat)%tphitphj(2:mesh_size,kln))& 437 & *rad(2:mesh_size)**ll 438 call simp_gen(intg,ff,pawrad(itypat)) 439 do mm=-ll,ll 440 isel=pawang%gntselect(lm0+mm,klm) 441 if (isel>0) pawtab(itypat)%qijl(lm0+mm,klmn)=intg*pawang%realgnt(isel) 442 end do 443 end do 444 if (klm_diag(klmn)==1) pawtab(itypat)%sij(klmn)= & 445 & pawtab(itypat)%qijl(1,klmn)*sqrt(four_pi) 446 end do 447 448 ! ================================================== 449 ! 6- COMPUTE Eijkl TERMS (Hartree) 450 ! Compute eventually short-range screened version of Eijkl (Fock) 451 452 if (allocated(pawtab(itypat)%eijkl)) then 453 ABI_DEALLOCATE(pawtab(itypat)%eijkl) 454 end if 455 ABI_ALLOCATE(pawtab(itypat)%eijkl,(lmn2_size,lmn2_size)) 456 if (abs(hyb_range_fock)>tol8) then 457 if (allocated(pawtab(itypat)%eijkl_sr)) then 458 ABI_DEALLOCATE(pawtab(itypat)%eijkl_sr) 459 end if 460 ABI_ALLOCATE(pawtab(itypat)%eijkl_sr,(lmn2_size,lmn2_size)) 461 end if 462 463 ! First loop is for eijkl (Hartree) 464 ! 2nd loop is for eijkl_sr (short-range screened Fock exchange) 465 do iloop=1,2 466 if (iloop==2.and.abs(hyb_range_fock)<=tol8) cycle 467 if (iloop==1) eijkl => pawtab(itypat)%eijkl 468 if (iloop==2) eijkl => pawtab(itypat)%eijkl_sr 469 470 ! Compute: 471 ! vhatL(r) according to eq. (A14) in Holzwarth et al., PRB 55, 2005 (1997) 472 ! intvhatL=$\int_{0}^{r_c}{vhatL(r) shapefunc_L(r) r^2\,dr}$ 473 ! vhatijL =$\int_{0}^{r_c}{vhatL(r) \tilde{\phi}_i \tilde{\phi}_j \,dr}$ 474 ! ----------------------------------------------------------------- 475 ABI_ALLOCATE(vhatl,(mesh_size)) 476 ABI_ALLOCATE(vhatijl,(lmn2_size,l_size)) 477 ABI_ALLOCATE(intvhatl,(l_size)) 478 intvhatl(:)=zero;vhatl(:)=zero;vhatijl(:,:)=zero 479 do il=1,l_size 480 vhatl(1)=zero;ff(1)=zero 481 ff(2:mesh_size)=pawtab(itypat)%shapefunc(2:mesh_size,il)*rad(2:mesh_size)**2 482 if (iloop==1) call poisson(ff,il-1,pawrad(itypat),vhatl) 483 if (iloop==2) call poisson(ff,il-1,pawrad(itypat),vhatl,screened_sr_separation=hyb_range_fock) 484 vhatl(2:mesh_size)=two*vhatl(2:mesh_size)/rad(2:mesh_size) 485 gg(1:mesh_size)=vhatl(1:mesh_size)*ff(1:mesh_size) 486 call simp_gen(intvhatl(il),gg,pawrad(itypat)) 487 do klmn=1,lmn2_size 488 kln=indklmn_(2,klmn) 489 hh(1:mesh_size)=vhatl(1:mesh_size)*pawtab(itypat)%tphitphj(1:mesh_size,kln) 490 call simp_gen(vhatijl(klmn,il),hh,pawrad(itypat)) 491 end do 492 end do 493 ABI_DEALLOCATE(vhatl) 494 495 ! Compute: 496 ! eijkl=$ vh1_ijkl - Vhatijkl - Bijkl - Cijkl$ 497 ! With: 498 ! $vh1_ijkl =\sum_{L,m} {vh1*Gaunt(i,j,Lm)*Gaunt(k,l,Lm)}$ 499 ! $Vhat_ijkl=\sum_{L,m} {vhatijL*Gaunt(i,j,Lm)*q_klL}$ 500 ! $B_ijkl =\sum_{L,m} {vhatijL*Gaunt(k,l,Lm)*q_ijL}$ 501 ! $C_ijkl =\sum_{L,m} {intvhatL*q_ijL*q_klL}$ 502 ! and: 503 ! vh1 according to eq. (A17) in Holzwarth et al., PRB 55, 2005 (1997) 504 ! Warning: compute only eijkl for (i,j)<=(k,l) 505 ! ----------------------------------------------------------------- 506 eijkl(:,:)=zero 507 meshsz=pawrad(itypat)%int_meshsz;if (mesh_size>meshsz) ff(meshsz+1:mesh_size)=zero 508 do klmn=1,lmn2_size 509 klm=indklmn_(1,klmn);kln=indklmn_(2,klmn) 510 lmin=indklmn_(3,klmn);lmax=indklmn_(4,klmn) 511 do ll=lmin,lmax,2 512 lm0=ll*ll+ll+1 513 ff(1:meshsz)=pawtab(itypat)%phiphj (1:meshsz,kln) 514 if (iloop==1) call poisson(ff,ll,pawrad(itypat),gg) 515 if (iloop==2) call poisson(ff,ll,pawrad(itypat),gg,screened_sr_separation=hyb_range_fock) 516 ff(1:meshsz)=pawtab(itypat)%tphitphj(1:meshsz,kln) 517 if (iloop==1) call poisson(ff,ll,pawrad(itypat),hh) 518 if (iloop==2) call poisson(ff,ll,pawrad(itypat),hh,screened_sr_separation=hyb_range_fock) 519 do klmn1=klmn,lmn2_size 520 klm1=indklmn_(1,klmn1);kln1=indklmn_(2,klmn1) 521 lmin1=indklmn_(3,klmn1);lmax1=indklmn_(4,klmn1) 522 vh1=zero 523 if ((ll.ge.lmin1).and.(ll.le.lmax1)) then 524 ff(1)=zero 525 ff(2:meshsz)=(pawtab(itypat)%phiphj (2:meshsz,kln1)*gg(2:meshsz)& 526 & -pawtab(itypat)%tphitphj(2:meshsz,kln1)*hh(2:meshsz))& 527 & *two/rad(2:meshsz) 528 call simp_gen(vh1,ff,pawrad(itypat)) 529 end if 530 do mm=-ll,ll 531 isel =pawang%gntselect(lm0+mm,klm) 532 isel1=pawang%gntselect(lm0+mm,klm1) 533 if (isel>0.and.isel1>0) then 534 rg =pawang%realgnt(isel) 535 rg1=pawang%realgnt(isel1) 536 ql =pawtab(itypat)%qijl(lm0+mm,klmn) 537 ql1=pawtab(itypat)%qijl(lm0+mm,klmn1) 538 eijkl(klmn,klmn1)=eijkl(klmn,klmn1)& 539 & +( vh1 *rg *rg1& ! vh1_ijkl 540 & - vhatijl(klmn ,ll+1)*rg *ql1& ! Vhat_ijkl 541 & - vhatijl(klmn1,ll+1)*rg1*ql & ! B_ijkl 542 & - intvhatl(ll+1) *ql *ql1& ! C_ijkl 543 & )*two_pi 544 end if 545 end do 546 end do 547 end do 548 end do 549 ABI_DEALLOCATE(vhatijl) 550 ABI_DEALLOCATE(intvhatl) 551 end do ! iloop 552 553 ! ================================================== 554 ! 7- COMPUTE gamma_ij TERMS 555 ! corrections to get the background right 556 557 if (pawtab(itypat)%usepotzero==1) then 558 if (allocated(pawtab(itypat)%gammaij)) then 559 ABI_DEALLOCATE(pawtab(itypat)%gammaij) 560 end if 561 ABI_ALLOCATE(pawtab(itypat)%gammaij,(lmn2_size)) 562 ABI_ALLOCATE(work,(mesh_size)) 563 do klmn=1,lmn2_size 564 if (klm_diag(klmn)==1) then 565 kln=indklmn_(2,klmn) 566 ff(1)=zero 567 ff(2:mesh_size)=pawtab(itypat)%phiphj(2:mesh_size,kln)-pawtab(itypat)%tphitphj(2:mesh_size,kln) 568 ! First, compute q_ij^00 569 call simp_gen(intg,ff,pawrad(itypat)) 570 ! Second, compute phi^2 - tphi^2 - 4pi*r^2 q_ij^00 g_0(r) 571 ff(2:mesh_size)= ff(2:mesh_size) - intg*pawtab(itypat)%shapefunc(2:mesh_size,1)*rad(2:mesh_size)**2 572 call poisson(ff,0,pawrad(itypat),work) 573 ! work is r*vh; should be then multiplied by r to prepare the integration in the sphere 574 work(1)=zero ; work(2:mesh_size)=work(2:mesh_size)*rad(2:mesh_size) 575 ! Third, average it over the sphere 576 call simp_gen(intg,work,pawrad(itypat)) 577 ! Finally, store it in pawtab%gammaij 578 pawtab(itypat)%gammaij(klmn)=intg*four_pi 579 else 580 pawtab(itypat)%gammaij(klmn)=zero 581 end if 582 end do 583 ABI_DEALLOCATE(work) 584 end if 585 586 ! *********************** 587 ! End Loop on atom types 588 ! *********************** 589 ABI_DEALLOCATE(ff) 590 ABI_DEALLOCATE(gg) 591 ABI_DEALLOCATE(hh) 592 ABI_DEALLOCATE(indklmn_) 593 ABI_DEALLOCATE(klm_diag) 594 ABI_DEALLOCATE(rad) 595 end do 596 597 call timab(553,2,tsec) 598 599 DBG_EXIT("COLL") 600 601 end subroutine pawinit