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
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}
SOURCE
94 subroutine initylmg(gprimd, kg, kptns, mkmem, mpi_enreg, mpsang, mpw, & 95 nband, nkpt, npwarr, nsppol, optder, rprimd, ylm, ylm_gr) 96 97 !Arguments ------------------------------------ 98 !scalars 99 integer,intent(in) :: mkmem,mpsang,mpw,nkpt,nsppol,optder 100 type(MPI_type),intent(in) :: mpi_enreg 101 !arrays 102 integer,intent(in) :: kg(3,mpw*mkmem),nband(nkpt*nsppol) 103 integer,intent(in) :: npwarr(nkpt) 104 real(dp),intent(in) :: gprimd(3,3),kptns(3,nkpt),rprimd(3,3) 105 real(dp),intent(out) :: ylm(mpw*mkmem,mpsang*mpsang) 106 real(dp),intent(out) :: ylm_gr(mpw*mkmem,3+6*(optder/2),mpsang*mpsang) 107 108 !Local variables ------------------------------ 109 !scalars 110 integer :: dimgr,ia,ib,ii,ikg,ikpt,ilang,ipw 111 integer :: jj,kk,l0,ll 112 integer :: me_distrb,mm,npw_k 113 real(dp),parameter :: tol=1.d-10 114 real(dp) :: cphi,ctheta,fact,onem,rr,sphi,stheta,work1,work2 115 real(dp) :: xx,ylmcst,ylmcst2 116 real(dp) :: yy,zz 117 !character(len=500) :: message 118 !arrays 119 integer,parameter :: alpha(6)=(/1,2,3,3,3,2/) 120 integer,parameter :: beta(6)=(/1,2,3,2,1,1/) 121 integer,allocatable :: kg_k(:,:) 122 real(dp) :: dphi(3),dtheta(3),iphase(mpsang-1),kpg(3) 123 real(dp) :: rphase(mpsang-1) 124 real(dp),allocatable :: blm(:,:) 125 real(dp),allocatable :: ylmgr2_cart(:,:,:),ylmgr2_tmp(:,:) 126 real(dp),allocatable :: ylmgr_cart(:,:) 127 real(dp),allocatable :: ylmgr_red(:,:) 128 129 !***************************************************************** 130 131 !Begin executable 132 me_distrb=mpi_enreg%me_kpt 133 !Initialisation of spherical harmonics (and gradients) 134 if (optder>=0) ylm(:,:) =zero 135 if (optder/=0) ylm_gr(:,:,:)=zero 136 dimgr=3+6*(optder/2) 137 138 !Allocate some memory 139 if (optder/=0) then 140 ABI_MALLOC(ylmgr_cart,(3,2)) 141 end if 142 if (optder/=0.and.optder/=2) then 143 ABI_MALLOC(ylmgr_red,(3,2)) 144 end if 145 if (optder==2) then 146 ABI_MALLOC(ylmgr2_cart,(3,3,2)) 147 ABI_MALLOC(ylmgr2_tmp,(3,3)) 148 ABI_MALLOC(ylmgr_red,(6,2)) 149 ABI_MALLOC(blm,(5,mpsang*mpsang)) 150 end if 151 152 !Loop over k-points: 153 ikg=0 154 do ikpt=1,nkpt 155 156 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband(ikpt),-1,me_distrb)) cycle 157 158 159 ! Get k+G-vectors, for this k-point: 160 npw_k=npwarr(ikpt) 161 ABI_MALLOC(kg_k,(3,npw_k)) 162 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 163 164 ! Special case for l=0 165 if (optder>=0) ylm(1+ikg:npw_k+ikg,1)=1._dp/sqrt(four_pi) 166 if (optder/=0) ylm_gr(1+ikg:npw_k+ikg,1:dimgr,1)=zero 167 168 if (mpsang>1) then 169 ! Loop over all k+G 170 do ipw=1,npw_k 171 172 ! Load k+G 173 kpg(1)=kptns(1,ikpt)+real(kg_k(1,ipw),dp) 174 kpg(2)=kptns(2,ikpt)+real(kg_k(2,ipw),dp) 175 kpg(3)=kptns(3,ikpt)+real(kg_k(3,ipw),dp) 176 177 ! Calculate module of k+G 178 xx=gprimd(1,1)*kpg(1)+gprimd(1,2)*kpg(2)+gprimd(1,3)*kpg(3) 179 yy=gprimd(2,1)*kpg(1)+gprimd(2,2)*kpg(2)+gprimd(2,3)*kpg(3) 180 zz=gprimd(3,1)*kpg(1)+gprimd(3,2)*kpg(2)+gprimd(3,3)*kpg(3) 181 rr=sqrt(xx**2+yy**2+zz**2) 182 183 ! Continue only for k+G<>0 184 if (rr>tol) then 185 186 ! Determine theta and phi 187 cphi=one 188 sphi=zero 189 ctheta=zz/rr 190 stheta=sqrt(abs((one-ctheta)*(one+ctheta))) 191 if (stheta>tol) then 192 cphi=xx/(rr*stheta) 193 sphi=yy/(rr*stheta) 194 end if 195 do mm=1,mpsang-1 196 rphase(mm)=dreal(dcmplx(cphi,sphi)**mm) 197 iphase(mm)=aimag(dcmplx(cphi,sphi)**mm) 198 end do 199 200 ! Determine gradients of theta and phi 201 if (optder/=0) then 202 dtheta(1)=ctheta*cphi 203 dtheta(2)=ctheta*sphi 204 dtheta(3)=-stheta 205 dphi(1)=-sphi 206 dphi(2)=cphi 207 dphi(3)=zero 208 end if 209 210 ! COMPUTE Ylm(K) 211 ! ============================================ 212 if (optder>=0) then 213 ! Loop over angular momentum l 214 do ilang=2,mpsang 215 ll=ilang-1 216 l0=ll**2+ll+1 217 fact=1._dp/real(ll*(ll+1),dp) 218 ylmcst=sqrt(real(2*ll+1,dp)/four_pi) 219 ! Special case m=0 220 ylm(ikg+ipw,l0)=ylmcst*ass_leg_pol(ll,0,ctheta) 221 ! Compute for m>0 222 onem=one 223 do mm=1,ll 224 onem=-onem 225 work1=ylmcst*sqrt(fact)*onem*ass_leg_pol(ll,mm,ctheta)*sqrt(2._dp) 226 ylm(ikg+ipw,l0+mm)=work1*rphase(mm) 227 ylm(ikg+ipw,l0-mm)=work1*iphase(mm) 228 if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp) 229 end do ! End loop over m 230 end do ! End loop over l 231 end if 232 233 ! COMPUTE dYlm/dKi 234 ! ============================================ 235 if (optder/=0) then 236 ! Loop over angular momentum l 237 do ilang=2,mpsang 238 ll=ilang-1 239 l0=ll**2+ll+1 240 fact=1._dp/real(ll*(ll+1),dp) 241 ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/rr 242 ! === Special case m=0 === 243 ! 1-compute gradients in cartesian coordinates 244 work1=ylmcst*plm_dtheta(ll,0,ctheta) 245 ylmgr_cart(1:3,1)=work1*dtheta(1:3) 246 ! 2-Transfer gradients into reduced coordinates 247 do ii=1,3 248 ylmgr_red(ii,1)=(rprimd(1,ii)*ylmgr_cart(1,1)+& 249 & rprimd(2,ii)*ylmgr_cart(2,1)+& 250 & rprimd(3,ii)*ylmgr_cart(3,1)) 251 end do 252 ! 3-Store gradients 253 ylm_gr(ikg+ipw,1:3,l0) =ylmgr_red(1:3,1) 254 ! === Compute for m>0 === 255 onem=one 256 do mm=1,ll 257 onem=-onem 258 ! 1-compute gradients in cartesian coordinates 259 work1=ylmcst*sqrt(fact)*onem*plm_dtheta(ll,mm,ctheta)*sqrt(2._dp) 260 work2=ylmcst*sqrt(fact)*onem*plm_dphi (ll,mm,ctheta)*sqrt(2._dp) 261 ylmgr_cart(1:3,1)=rphase(mm)*work1*dtheta(1:3)-iphase(mm)*work2*dphi(1:3) 262 ylmgr_cart(1:3,2)=iphase(mm)*work1*dtheta(1:3)+rphase(mm)*work2*dphi(1:3) 263 ! 2-Transfer gradients into reduced coordinates 264 do kk=1,2 265 do ii=1,3 266 ylmgr_red(ii,kk)=(rprimd(1,ii)*ylmgr_cart(1,kk)+& 267 & rprimd(2,ii)*ylmgr_cart(2,kk)+& 268 & rprimd(3,ii)*ylmgr_cart(3,kk)) 269 end do 270 end do 271 ! 3-Store gradients 272 ylm_gr(ikg+ipw,1:3,l0+mm) =ylmgr_red(1:3,1) 273 ylm_gr(ikg+ipw,1:3,l0-mm) =ylmgr_red(1:3,2) 274 if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp) 275 end do ! End loop over m 276 end do ! End loop over l 277 end if 278 279 ! COMPUTE d2Ylm/dKidKj 280 ! ============================================ 281 if (optder==2) then 282 call plm_coeff(blm,mpsang,ctheta) 283 ! Loop over angular momentum l 284 do ilang=2,mpsang 285 ll=ilang-1 286 l0=ll**2+ll+1 287 fact=1._dp/real(ll*(ll+1),dp) 288 ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/(rr**2) 289 ! === Special case m=0 === 290 ! 1-compute gradients in cartesian coordinates 291 ylmgr2_cart(1,1,1)=ylmcst*(-blm(3,l0)*sphi*sphi+blm(4,l0)*cphi*cphi) 292 ylmgr2_cart(2,2,1)=ylmcst*(-blm(3,l0)*cphi*cphi+blm(4,l0)*sphi*sphi) 293 ylmgr2_cart(3,3,1)=ylmcst*blm(1,l0) 294 ylmgr2_cart(3,1,1)=ylmcst*blm(2,l0)*cphi 295 ylmgr2_cart(3,2,1)=ylmcst*blm(2,l0)*sphi 296 ylmgr2_cart(2,1,1)=ylmcst*(blm(3,l0)+blm(4,l0))*sphi*cphi 297 ylmgr2_cart(1,3,1)=ylmgr2_cart(3,1,1) 298 ylmgr2_cart(1,2,1)=ylmgr2_cart(2,1,1) 299 ylmgr2_cart(2,3,1)=ylmgr2_cart(3,2,1) 300 ! 2-Transfer gradients into reduced coordinates 301 do jj=1,3 302 do ii=1,3 303 ylmgr2_tmp(ii,jj)=(rprimd(1,jj)*ylmgr2_cart(1,ii,1)+& 304 & rprimd(2,jj)*ylmgr2_cart(2,ii,1)+& 305 & rprimd(3,jj)*ylmgr2_cart(3,ii,1)) 306 end do 307 end do 308 do ii=1,6 309 ia=alpha(ii);ib=beta(ii) 310 ylmgr_red(ii,1)=(rprimd(1,ia)*ylmgr2_tmp(1,ib)+& 311 & rprimd(2,ia)*ylmgr2_tmp(2,ib)+& 312 & rprimd(3,ia)*ylmgr2_tmp(3,ib)) 313 end do 314 ylm_gr(ikg+ipw,4:9,l0) =ylmgr_red(1:6,1) 315 ! === Compute for m>0 === 316 onem=one 317 do mm=1,ll 318 onem=-onem;ylmcst2=ylmcst*sqrt(fact)*sqrt(two) 319 ylmgr2_cart(1,1,1)=ylmcst2*((-blm(3,l0+mm)*sphi*sphi+blm(4,l0+mm)*cphi*cphi)*rphase(mm)-& 320 & blm(5,l0+mm)*2.d0*cphi*sphi*mm*iphase(mm)) 321 ylmgr2_cart(1,1,2)=ylmcst2*((-blm(3,l0+mm)*sphi*sphi+blm(4,l0+mm)*cphi*cphi)*iphase(mm)+& 322 & blm(5,l0+mm)*2.d0*cphi*sphi*mm*rphase(mm)) 323 ylmgr2_cart(2,2,1)=ylmcst2*((-blm(3,l0+mm)*cphi*cphi+blm(4,l0+mm)*sphi*sphi)*rphase(mm)+& 324 & blm(5,l0+mm)*2.d0*cphi*sphi*mm*iphase(mm)) 325 ylmgr2_cart(2,2,2)=ylmcst2*((-blm(3,l0+mm)*cphi*cphi+blm(4,l0+mm)*sphi*sphi)*iphase(mm)-& 326 & blm(5,l0+mm)*2.d0*cphi*sphi*mm*rphase(mm)) 327 ylmgr2_cart(3,3,1)=ylmcst2*blm(1,l0+mm)*rphase(mm) 328 ylmgr2_cart(3,3,2)=ylmcst2*blm(1,l0+mm)*iphase(mm) 329 ylmgr2_cart(3,1,1)=ylmcst2*(blm(2,l0+mm)*cphi*rphase(mm)-& 330 & mm*iphase(mm)*sphi*onem*plm_dtheta(ll,mm,ctheta)) 331 ylmgr2_cart(3,1,2)=ylmcst2*(blm(2,l0+mm)*cphi*iphase(mm)+& 332 & mm*rphase(mm)*sphi*onem*plm_dtheta(ll,mm,ctheta)) 333 ylmgr2_cart(3,2,1)=ylmcst2*(blm(2,l0+mm)*sphi*rphase(mm)+& 334 & mm*iphase(mm)*cphi*onem*plm_dtheta(ll,mm,ctheta)) 335 ylmgr2_cart(3,2,2)=ylmcst2*(blm(2,l0+mm)*sphi*iphase(mm)-& 336 & mm*rphase(mm)*cphi*onem*plm_dtheta(ll,mm,ctheta)) 337 ylmgr2_cart(2,1,1)=ylmcst2*((blm(3,l0+mm)+blm(4,l0+mm))*sphi*cphi*rphase(mm)-& 338 & blm(5,l0+mm)*(sphi*sphi-cphi*cphi)*mm*iphase(mm)) 339 ylmgr2_cart(2,1,2)=ylmcst2*((blm(3,l0+mm)+blm(4,l0+mm))*sphi*cphi*iphase(mm)+& 340 & blm(5,l0+mm)*(sphi*sphi-cphi*cphi)*mm*rphase(mm)) 341 ylmgr2_cart(1,3,:)=ylmgr2_cart(3,1,:) 342 ylmgr2_cart(1,2,:)=ylmgr2_cart(2,1,:) 343 ylmgr2_cart(2,3,:)=ylmgr2_cart(3,2,:) 344 ! 2-Transfer gradients into reduced coordinates 345 do kk=1,2 346 do jj=1,3 347 do ii=1,3 348 ylmgr2_tmp(ii,jj)=(rprimd(1,jj)*ylmgr2_cart(1,ii,kk)+& 349 & rprimd(2,jj)*ylmgr2_cart(2,ii,kk)+& 350 & rprimd(3,jj)*ylmgr2_cart(3,ii,kk)) 351 end do 352 end do 353 do ii=1,6 354 ia=alpha(ii);ib=beta(ii) 355 ylmgr_red(ii,kk)=(rprimd(1,ia)*ylmgr2_tmp(1,ib)+& 356 & rprimd(2,ia)*ylmgr2_tmp(2,ib)+& 357 & rprimd(3,ia)*ylmgr2_tmp(3,ib)) 358 end do 359 end do 360 ylm_gr(ikg+ipw,4:9,l0+mm) =ylmgr_red(1:6,1) 361 ylm_gr(ikg+ipw,4:9,l0-mm) =ylmgr_red(1:6,2) 362 if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp) 363 end do ! End loop over m 364 end do ! End loop over l 365 end if 366 367 ! End condition r<>0 368 end if 369 370 ! End loop over k+G 371 end do 372 373 ! End condition l<>0 374 end if 375 376 ABI_FREE(kg_k) 377 378 ikg=ikg+npw_k 379 end do ! End Loop over k-points 380 381 !Release the temporary memory 382 !Allocate some memory 383 if (optder/=0) then 384 ABI_FREE(ylmgr_cart) 385 end if 386 if (optder/=0.and.optder/=2) then 387 ABI_FREE(ylmgr_red) 388 end if 389 if (optder==2) then 390 ABI_FREE(ylmgr2_cart) 391 ABI_FREE(ylmgr2_tmp) 392 ABI_FREE(ylmgr_red) 393 ABI_FREE(blm) 394 end if 395 396 end subroutine initylmg
ABINIT/initylmg_k [ Functions ]
NAME
initylmg_k
FUNCTION
Simplified interface to compute Ylm for a single k-point. See initylmg for the the meaning of the input variables.
SOURCE
409 subroutine initylmg_k(npw, mpsang, optder, rprimd, gprimd, kpt, kg, ylm, ylm_gr) 410 411 !Arguments ------------------------------------ 412 !scalars 413 integer,intent(in) :: mpsang, npw, optder 414 !arrays 415 real(dp),intent(in) :: kpt(3) 416 integer,intent(in) :: kg(3,npw) 417 real(dp),intent(in) :: gprimd(3,3), rprimd(3,3) 418 real(dp),intent(out) :: ylm(npw, mpsang*mpsang) 419 real(dp),intent(out) :: ylm_gr(npw, 3+6*(optder/2), mpsang*mpsang) 420 421 !Local variables ------------------------------ 422 !scalars 423 integer,parameter :: nkpt_1 = 1, nsppol_1 = 1 424 integer :: npwarr__(nkpt_1), nband__(nkpt_1 * nsppol_1) 425 real(dp) :: kptns__(3,nkpt_1) 426 type(MPI_type) :: seq_mpi_enreg 427 428 !***************************************************************** 429 430 call initmpi_seq(seq_mpi_enreg) 431 npwarr__(1) = npw 432 kptns__(:,1) = kpt 433 nband__ = 1 ! Not used in sequential 434 435 call initylmg(gprimd, kg, kptns__, nkpt_1, seq_mpi_enreg, mpsang, npw, & 436 nband__, nkpt_1, npwarr__, nsppol_1, optder, rprimd, ylm, ylm_gr) 437 438 call destroy_mpi_enreg(seq_mpi_enreg) 439 440 end subroutine initylmg_k
ABINIT/m_initylmg [ Modules ]
NAME
m_initylmg
FUNCTION
Calculate the real spherical harmonics Ylm (and gradients) over a set of (reciprocal space) (k+G) vectors
COPYRIGHT
Copyright (C) 1998-2024 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
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 module m_initylmg 24 25 use defs_basis 26 use m_abicore 27 use m_errors 28 use m_xmpi 29 30 use defs_abitypes, only : MPI_type 31 use m_paw_sphharm, only : ass_leg_pol, plm_dtheta, plm_dphi, plm_coeff 32 use m_mpinfo, only : proc_distrb_cycle, destroy_mpi_enreg, initmpi_seq 33 34 implicit none 35 36 private