TABLE OF CONTENTS
ABINIT/initylmg [ Functions ]
NAME
initylmg
FUNCTION
Calculate the real spherical harmonics Ylm (and gradients) over a set of (reciprocal space) (k+G) vectors
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
gprimd(3,3)=dimensional reciprocal space primitive translations (b^-1) kg(3,mpw)=integer coordinates of G vectors in basis sphere kptns(3,nkpt)=k points in terms of reciprocal translations mkmem =number of k points treated by this node mpi_enreg=information about MPI parallelization mpsang=1+maximum angular momentum for nonlocal pseudopotential mpw =maximum number of planewaves in basis sphere (large number) nband(nkpt*nsppol)=number of bands at each k point nkpt =number of k points npwarr(nkpt)=array holding npw for each k point nsppol=1 for unpolarized, 2 for polarized optder= 0=compute Ylm(K) 1=compute Ylm(K) and dYlm/dKi 2=compute Ylm(K), dYlm/dKi and d2Ylm/dKidKj -1=compute only dYlm/dKi rprimd(3,3)=dimensional primitive translations in real space (bohr)
OUTPUT
if (optder>=0)
ylm(mpw*mkmem,mpsang*mpsang) = real spherical harmonics
for each G and k point
if (optder>=1 or optder==-1)
ylm_gr(mpw*mkmem,1:3,mpsang*mpsang)= gradients of real
spherical harmonics wrt (G+k) in reduced coordinates
if (optder>=2)
ylm_gr(mpw*mkmem,4:9,mpsang*mpsang)= second gradients of
real spherical harmonics wrt (G+k) in reduced coordinates
NOTES
Remember the expression of complex spherical harmonics:
$Y_{lm}(%theta ,%phi)=sqrt{{(2l+1) over (4 %pi)}
{fact(l-m) over fact(l+m)} } P_l^m(cos(%theta))
func e^{i m %phi}$
Remember the expression of real spherical harmonics as
linear combination of complex spherical harmonics:
$Yr_{lm}(%theta ,%phi)=(Re{Y_{l-m}}+(-1)^m Re{Y_{lm}})/sqrt{2}
$Yr_{l-m}(%theta ,%phi)=(Im{Y_{l-m}}-(-1)^m Im{Y_{lm}})/sqrt{2}
PARENTS
debug_tools,dfpt_looppert,dfpt_nstpaw,dfptnl_loop,forstr,gstate ks_ddiago,m_cut3d,m_pawpwij,m_shirley,m_wfd,mover,nonlop_test partial_dos_fractions,respfn,scfcv
CHILDREN
plm_coeff
SOURCE
70 #if defined HAVE_CONFIG_H 71 #include "config.h" 72 #endif 73 74 #include "abi_common.h" 75 76 subroutine initylmg(gprimd,kg,kptns,mkmem,mpi_enreg,mpsang,mpw,& 77 & nband,nkpt,npwarr,nsppol,optder,rprimd,ylm,ylm_gr) 78 79 use defs_basis 80 use defs_abitypes 81 use m_profiling_abi 82 use m_errors 83 use m_xmpi 84 85 use m_paw_sphharm, only : ass_leg_pol, plm_dtheta, plm_dphi, plm_coeff 86 87 !This section has been created automatically by the script Abilint (TD). 88 !Do not modify the following lines by hand. 89 #undef ABI_FUNC 90 #define ABI_FUNC 'initylmg' 91 use interfaces_32_util 92 !End of the abilint section 93 94 implicit none 95 96 !Arguments ------------------------------------ 97 !scalars 98 integer,intent(in) :: mkmem,mpsang,mpw,nkpt,nsppol,optder 99 type(MPI_type),intent(in) :: mpi_enreg 100 !arrays 101 integer,intent(in) :: kg(3,mpw*mkmem),nband(nkpt*nsppol) 102 integer,intent(in) :: npwarr(nkpt) 103 real(dp),intent(in) :: gprimd(3,3),kptns(3,nkpt),rprimd(3,3) 104 real(dp),intent(out) :: ylm(mpw*mkmem,mpsang*mpsang) 105 real(dp),intent(out) :: ylm_gr(mpw*mkmem,3+6*(optder/2),mpsang*mpsang) 106 107 !Local variables ------------------------------ 108 !scalars 109 integer :: dimgr,ia,ib,ii,ikg,ikpt,ilang,ipw 110 integer :: jj,kk,l0,ll 111 integer :: me_distrb,mm,npw_k 112 real(dp),parameter :: tol=1.d-10 113 real(dp) :: cphi,ctheta,fact,onem,rr,sphi,stheta,work1,work2 114 real(dp) :: xx,ylmcst,ylmcst2 115 real(dp) :: yy,zz 116 !character(len=500) :: message 117 !arrays 118 integer,parameter :: alpha(6)=(/1,2,3,3,3,2/) 119 integer,parameter :: beta(6)=(/1,2,3,2,1,1/) 120 integer,allocatable :: kg_k(:,:) 121 real(dp) :: dphi(3),dtheta(3),iphase(mpsang-1),kpg(3) 122 real(dp) :: rphase(mpsang-1) 123 real(dp),allocatable :: blm(:,:) 124 real(dp),allocatable :: ylmgr2_cart(:,:,:),ylmgr2_tmp(:,:) 125 real(dp),allocatable :: ylmgr_cart(:,:) 126 real(dp),allocatable :: ylmgr_red(:,:) 127 128 !***************************************************************** 129 130 !Begin executable 131 me_distrb=mpi_enreg%me_kpt 132 !Initialisation of spherical harmonics (and gradients) 133 if (optder>=0) ylm(:,:) =zero 134 if (optder/=0) ylm_gr(:,:,:)=zero 135 dimgr=3+6*(optder/2) 136 137 !Allocate some memory 138 if (optder/=0) then 139 ABI_ALLOCATE(ylmgr_cart,(3,2)) 140 end if 141 if (optder/=0.and.optder/=2) then 142 ABI_ALLOCATE(ylmgr_red,(3,2)) 143 end if 144 if (optder==2) then 145 ABI_ALLOCATE(ylmgr2_cart,(3,3,2)) 146 ABI_ALLOCATE(ylmgr2_tmp,(3,3)) 147 ABI_ALLOCATE(ylmgr_red,(6,2)) 148 ABI_ALLOCATE(blm,(5,mpsang*mpsang)) 149 end if 150 151 !Loop over k-points: 152 ikg=0 153 do ikpt=1,nkpt 154 155 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband(ikpt),-1,me_distrb)) cycle 156 157 158 ! Get k+G-vectors, for this k-point: 159 npw_k=npwarr(ikpt) 160 ABI_ALLOCATE(kg_k,(3,npw_k)) 161 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 162 163 ! Special case for l=0 164 if (optder>=0) ylm(1+ikg:npw_k+ikg,1)=1._dp/sqrt(four_pi) 165 if (optder/=0) ylm_gr(1+ikg:npw_k+ikg,1:dimgr,1)=zero 166 167 if (mpsang>1) then 168 ! Loop over all k+G 169 do ipw=1,npw_k 170 171 ! Load k+G 172 kpg(1)=kptns(1,ikpt)+real(kg_k(1,ipw),dp) 173 kpg(2)=kptns(2,ikpt)+real(kg_k(2,ipw),dp) 174 kpg(3)=kptns(3,ikpt)+real(kg_k(3,ipw),dp) 175 176 ! Calculate module of k+G 177 xx=gprimd(1,1)*kpg(1)+gprimd(1,2)*kpg(2)+gprimd(1,3)*kpg(3) 178 yy=gprimd(2,1)*kpg(1)+gprimd(2,2)*kpg(2)+gprimd(2,3)*kpg(3) 179 zz=gprimd(3,1)*kpg(1)+gprimd(3,2)*kpg(2)+gprimd(3,3)*kpg(3) 180 rr=sqrt(xx**2+yy**2+zz**2) 181 182 ! Continue only for k+G<>0 183 if (rr>tol) then 184 185 ! Determine theta and phi 186 cphi=one 187 sphi=zero 188 ctheta=zz/rr 189 stheta=sqrt(abs((one-ctheta)*(one+ctheta))) 190 if (stheta>tol) then 191 cphi=xx/(rr*stheta) 192 sphi=yy/(rr*stheta) 193 end if 194 do mm=1,mpsang-1 195 rphase(mm)=dreal(dcmplx(cphi,sphi)**mm) 196 iphase(mm)=aimag(dcmplx(cphi,sphi)**mm) 197 end do 198 199 ! Determine gradients of theta and phi 200 if (optder/=0) then 201 dtheta(1)=ctheta*cphi 202 dtheta(2)=ctheta*sphi 203 dtheta(3)=-stheta 204 dphi(1)=-sphi 205 dphi(2)=cphi 206 dphi(3)=zero 207 end if 208 209 ! COMPUTE Ylm(K) 210 ! ============================================ 211 if (optder>=0) then 212 ! Loop over angular momentum l 213 do ilang=2,mpsang 214 ll=ilang-1 215 l0=ll**2+ll+1 216 fact=1._dp/real(ll*(ll+1),dp) 217 ylmcst=sqrt(real(2*ll+1,dp)/four_pi) 218 ! Special case m=0 219 ylm(ikg+ipw,l0)=ylmcst*ass_leg_pol(ll,0,ctheta) 220 ! Compute for m>0 221 onem=one 222 do mm=1,ll 223 onem=-onem 224 work1=ylmcst*sqrt(fact)*onem*ass_leg_pol(ll,mm,ctheta)*sqrt(2._dp) 225 ylm(ikg+ipw,l0+mm)=work1*rphase(mm) 226 ylm(ikg+ipw,l0-mm)=work1*iphase(mm) 227 if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp) 228 end do ! End loop over m 229 end do ! End loop over l 230 end if 231 232 ! COMPUTE dYlm/dKi 233 ! ============================================ 234 if (optder/=0) then 235 ! Loop over angular momentum l 236 do ilang=2,mpsang 237 ll=ilang-1 238 l0=ll**2+ll+1 239 fact=1._dp/real(ll*(ll+1),dp) 240 ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/rr 241 ! === Special case m=0 === 242 ! 1-compute gradients in cartesian coordinates 243 work1=ylmcst*plm_dtheta(ll,0,ctheta) 244 ylmgr_cart(1:3,1)=work1*dtheta(1:3) 245 ! 2-Transfer gradients into reduced coordinates 246 do ii=1,3 247 ylmgr_red(ii,1)=(rprimd(1,ii)*ylmgr_cart(1,1)+& 248 & rprimd(2,ii)*ylmgr_cart(2,1)+& 249 & rprimd(3,ii)*ylmgr_cart(3,1)) 250 end do 251 ! 3-Store gradients 252 ylm_gr(ikg+ipw,1:3,l0) =ylmgr_red(1:3,1) 253 ! === Compute for m>0 === 254 onem=one 255 do mm=1,ll 256 onem=-onem 257 ! 1-compute gradients in cartesian coordinates 258 work1=ylmcst*sqrt(fact)*onem*plm_dtheta(ll,mm,ctheta)*sqrt(2._dp) 259 work2=ylmcst*sqrt(fact)*onem*plm_dphi (ll,mm,ctheta)*sqrt(2._dp) 260 ylmgr_cart(1:3,1)=rphase(mm)*work1*dtheta(1:3)-iphase(mm)*work2*dphi(1:3) 261 ylmgr_cart(1:3,2)=iphase(mm)*work1*dtheta(1:3)+rphase(mm)*work2*dphi(1:3) 262 ! 2-Transfer gradients into reduced coordinates 263 do kk=1,2 264 do ii=1,3 265 ylmgr_red(ii,kk)=(rprimd(1,ii)*ylmgr_cart(1,kk)+& 266 & rprimd(2,ii)*ylmgr_cart(2,kk)+& 267 & rprimd(3,ii)*ylmgr_cart(3,kk)) 268 end do 269 end do 270 ! 3-Store gradients 271 ylm_gr(ikg+ipw,1:3,l0+mm) =ylmgr_red(1:3,1) 272 ylm_gr(ikg+ipw,1:3,l0-mm) =ylmgr_red(1:3,2) 273 if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp) 274 end do ! End loop over m 275 end do ! End loop over l 276 end if 277 278 ! COMPUTE d2Ylm/dKidKj 279 ! ============================================ 280 if (optder==2) then 281 call plm_coeff(blm,mpsang,ctheta) 282 ! Loop over angular momentum l 283 do ilang=2,mpsang 284 ll=ilang-1 285 l0=ll**2+ll+1 286 fact=1._dp/real(ll*(ll+1),dp) 287 ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/(rr**2) 288 ! === Special case m=0 === 289 ! 1-compute gradients in cartesian coordinates 290 ylmgr2_cart(1,1,1)=ylmcst*(-blm(3,l0)*sphi*sphi+blm(4,l0)*cphi*cphi) 291 ylmgr2_cart(2,2,1)=ylmcst*(-blm(3,l0)*cphi*cphi+blm(4,l0)*sphi*sphi) 292 ylmgr2_cart(3,3,1)=ylmcst*blm(1,l0) 293 ylmgr2_cart(3,1,1)=ylmcst*blm(2,l0)*cphi 294 ylmgr2_cart(3,2,1)=ylmcst*blm(2,l0)*sphi 295 ylmgr2_cart(2,1,1)=ylmcst*(blm(3,l0)+blm(4,l0))*sphi*cphi 296 ylmgr2_cart(1,3,1)=ylmgr2_cart(3,1,1) 297 ylmgr2_cart(1,2,1)=ylmgr2_cart(2,1,1) 298 ylmgr2_cart(2,3,1)=ylmgr2_cart(3,2,1) 299 ! 2-Transfer gradients into reduced coordinates 300 do jj=1,3 301 do ii=1,3 302 ylmgr2_tmp(ii,jj)=(rprimd(1,jj)*ylmgr2_cart(1,ii,1)+& 303 & rprimd(2,jj)*ylmgr2_cart(2,ii,1)+& 304 & rprimd(3,jj)*ylmgr2_cart(3,ii,1)) 305 end do 306 end do 307 do ii=1,6 308 ia=alpha(ii);ib=beta(ii) 309 ylmgr_red(ii,1)=(rprimd(1,ia)*ylmgr2_tmp(1,ib)+& 310 & rprimd(2,ia)*ylmgr2_tmp(2,ib)+& 311 & rprimd(3,ia)*ylmgr2_tmp(3,ib)) 312 end do 313 ylm_gr(ikg+ipw,4:9,l0) =ylmgr_red(1:6,1) 314 ! === Compute for m>0 === 315 onem=one 316 do mm=1,ll 317 onem=-onem;ylmcst2=ylmcst*sqrt(fact)*sqrt(two) 318 ylmgr2_cart(1,1,1)=ylmcst2*((-blm(3,l0+mm)*sphi*sphi+blm(4,l0+mm)*cphi*cphi)*rphase(mm)-& 319 & blm(5,l0+mm)*2.d0*cphi*sphi*mm*iphase(mm)) 320 ylmgr2_cart(1,1,2)=ylmcst2*((-blm(3,l0+mm)*sphi*sphi+blm(4,l0+mm)*cphi*cphi)*iphase(mm)+& 321 & blm(5,l0+mm)*2.d0*cphi*sphi*mm*rphase(mm)) 322 ylmgr2_cart(2,2,1)=ylmcst2*((-blm(3,l0+mm)*cphi*cphi+blm(4,l0+mm)*sphi*sphi)*rphase(mm)+& 323 & blm(5,l0+mm)*2.d0*cphi*sphi*mm*iphase(mm)) 324 ylmgr2_cart(2,2,2)=ylmcst2*((-blm(3,l0+mm)*cphi*cphi+blm(4,l0+mm)*sphi*sphi)*iphase(mm)-& 325 & blm(5,l0+mm)*2.d0*cphi*sphi*mm*rphase(mm)) 326 ylmgr2_cart(3,3,1)=ylmcst2*blm(1,l0+mm)*rphase(mm) 327 ylmgr2_cart(3,3,2)=ylmcst2*blm(1,l0+mm)*iphase(mm) 328 ylmgr2_cart(3,1,1)=ylmcst2*(blm(2,l0+mm)*cphi*rphase(mm)-& 329 & mm*iphase(mm)*sphi*onem*plm_dtheta(ll,mm,ctheta)) 330 ylmgr2_cart(3,1,2)=ylmcst2*(blm(2,l0+mm)*cphi*iphase(mm)+& 331 & mm*rphase(mm)*sphi*onem*plm_dtheta(ll,mm,ctheta)) 332 ylmgr2_cart(3,2,1)=ylmcst2*(blm(2,l0+mm)*sphi*rphase(mm)+& 333 & mm*iphase(mm)*cphi*onem*plm_dtheta(ll,mm,ctheta)) 334 ylmgr2_cart(3,2,2)=ylmcst2*(blm(2,l0+mm)*sphi*iphase(mm)-& 335 & mm*rphase(mm)*cphi*onem*plm_dtheta(ll,mm,ctheta)) 336 ylmgr2_cart(2,1,1)=ylmcst2*((blm(3,l0+mm)+blm(4,l0+mm))*sphi*cphi*rphase(mm)-& 337 & blm(5,l0+mm)*(sphi*sphi-cphi*cphi)*mm*iphase(mm)) 338 ylmgr2_cart(2,1,2)=ylmcst2*((blm(3,l0+mm)+blm(4,l0+mm))*sphi*cphi*iphase(mm)+& 339 & blm(5,l0+mm)*(sphi*sphi-cphi*cphi)*mm*rphase(mm)) 340 ylmgr2_cart(1,3,:)=ylmgr2_cart(3,1,:) 341 ylmgr2_cart(1,2,:)=ylmgr2_cart(2,1,:) 342 ylmgr2_cart(2,3,:)=ylmgr2_cart(3,2,:) 343 ! 2-Transfer gradients into reduced coordinates 344 do kk=1,2 345 do jj=1,3 346 do ii=1,3 347 ylmgr2_tmp(ii,jj)=(rprimd(1,jj)*ylmgr2_cart(1,ii,kk)+& 348 & rprimd(2,jj)*ylmgr2_cart(2,ii,kk)+& 349 & rprimd(3,jj)*ylmgr2_cart(3,ii,kk)) 350 end do 351 end do 352 do ii=1,6 353 ia=alpha(ii);ib=beta(ii) 354 ylmgr_red(ii,kk)=(rprimd(1,ia)*ylmgr2_tmp(1,ib)+& 355 & rprimd(2,ia)*ylmgr2_tmp(2,ib)+& 356 & rprimd(3,ia)*ylmgr2_tmp(3,ib)) 357 end do 358 end do 359 ylm_gr(ikg+ipw,4:9,l0+mm) =ylmgr_red(1:6,1) 360 ylm_gr(ikg+ipw,4:9,l0-mm) =ylmgr_red(1:6,2) 361 if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp) 362 end do ! End loop over m 363 end do ! End loop over l 364 end if 365 366 ! End condition r<>0 367 end if 368 369 ! End loop over k+G 370 end do 371 372 ! End condition l<>0 373 end if 374 375 ABI_DEALLOCATE(kg_k) 376 377 ikg=ikg+npw_k 378 end do ! End Loop over k-points 379 380 !Release the temporary memory 381 !Allocate some memory 382 if (optder/=0) then 383 ABI_DEALLOCATE(ylmgr_cart) 384 end if 385 if (optder/=0.and.optder/=2) then 386 ABI_DEALLOCATE(ylmgr_red) 387 end if 388 if (optder==2) then 389 ABI_DEALLOCATE(ylmgr2_cart) 390 ABI_DEALLOCATE(ylmgr2_tmp) 391 ABI_DEALLOCATE(ylmgr_red) 392 ABI_DEALLOCATE(blm) 393 end if 394 395 end subroutine initylmg