TABLE OF CONTENTS


ABINIT/initylmg [ Functions ]

[ Top ] [ 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}

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

110 subroutine initylmg(gprimd,kg,kptns,mkmem,mpi_enreg,mpsang,mpw,&
111 &  nband,nkpt,npwarr,nsppol,optder,rprimd,ylm,ylm_gr)
112 
113 
114 !This section has been created automatically by the script Abilint (TD).
115 !Do not modify the following lines by hand.
116 #undef ABI_FUNC
117 #define ABI_FUNC 'initylmg'
118 !End of the abilint section
119 
120  implicit none
121 
122 !Arguments ------------------------------------
123 !scalars
124  integer,intent(in) :: mkmem,mpsang,mpw,nkpt,nsppol,optder
125  type(MPI_type),intent(in) :: mpi_enreg
126 !arrays
127  integer,intent(in) :: kg(3,mpw*mkmem),nband(nkpt*nsppol)
128  integer,intent(in) :: npwarr(nkpt)
129  real(dp),intent(in) :: gprimd(3,3),kptns(3,nkpt),rprimd(3,3)
130  real(dp),intent(out) :: ylm(mpw*mkmem,mpsang*mpsang)
131  real(dp),intent(out) :: ylm_gr(mpw*mkmem,3+6*(optder/2),mpsang*mpsang)
132 
133 !Local variables ------------------------------
134 !scalars
135  integer :: dimgr,ia,ib,ii,ikg,ikpt,ilang,ipw
136  integer :: jj,kk,l0,ll
137  integer :: me_distrb,mm,npw_k
138  real(dp),parameter :: tol=1.d-10
139  real(dp) :: cphi,ctheta,fact,onem,rr,sphi,stheta,work1,work2
140  real(dp) :: xx,ylmcst,ylmcst2
141  real(dp) :: yy,zz
142  !character(len=500) :: message
143 !arrays
144  integer,parameter :: alpha(6)=(/1,2,3,3,3,2/)
145  integer,parameter :: beta(6)=(/1,2,3,2,1,1/)
146  integer,allocatable :: kg_k(:,:)
147  real(dp) :: dphi(3),dtheta(3),iphase(mpsang-1),kpg(3)
148  real(dp) :: rphase(mpsang-1)
149  real(dp),allocatable :: blm(:,:)
150  real(dp),allocatable :: ylmgr2_cart(:,:,:),ylmgr2_tmp(:,:)
151  real(dp),allocatable :: ylmgr_cart(:,:)
152  real(dp),allocatable :: ylmgr_red(:,:)
153 
154 !*****************************************************************
155 
156 !Begin executable
157  me_distrb=mpi_enreg%me_kpt
158 !Initialisation of spherical harmonics (and gradients)
159  if (optder>=0) ylm(:,:)  =zero
160  if (optder/=0) ylm_gr(:,:,:)=zero
161  dimgr=3+6*(optder/2)
162 
163 !Allocate some memory
164  if (optder/=0) then
165    ABI_ALLOCATE(ylmgr_cart,(3,2))
166  end if
167  if (optder/=0.and.optder/=2) then
168    ABI_ALLOCATE(ylmgr_red,(3,2))
169  end if
170  if (optder==2) then
171    ABI_ALLOCATE(ylmgr2_cart,(3,3,2))
172    ABI_ALLOCATE(ylmgr2_tmp,(3,3))
173    ABI_ALLOCATE(ylmgr_red,(6,2))
174    ABI_ALLOCATE(blm,(5,mpsang*mpsang))
175  end if
176 
177 !Loop over k-points:
178  ikg=0
179  do ikpt=1,nkpt
180 
181    if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband(ikpt),-1,me_distrb)) cycle
182 
183 
184 !  Get k+G-vectors, for this k-point:
185    npw_k=npwarr(ikpt)
186    ABI_ALLOCATE(kg_k,(3,npw_k))
187    kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
188 
189 !  Special case for l=0
190    if (optder>=0) ylm(1+ikg:npw_k+ikg,1)=1._dp/sqrt(four_pi)
191    if (optder/=0) ylm_gr(1+ikg:npw_k+ikg,1:dimgr,1)=zero
192 
193    if (mpsang>1) then
194 !    Loop over all k+G
195      do ipw=1,npw_k
196 
197 !      Load k+G
198        kpg(1)=kptns(1,ikpt)+real(kg_k(1,ipw),dp)
199        kpg(2)=kptns(2,ikpt)+real(kg_k(2,ipw),dp)
200        kpg(3)=kptns(3,ikpt)+real(kg_k(3,ipw),dp)
201 
202 !      Calculate module of k+G
203        xx=gprimd(1,1)*kpg(1)+gprimd(1,2)*kpg(2)+gprimd(1,3)*kpg(3)
204        yy=gprimd(2,1)*kpg(1)+gprimd(2,2)*kpg(2)+gprimd(2,3)*kpg(3)
205        zz=gprimd(3,1)*kpg(1)+gprimd(3,2)*kpg(2)+gprimd(3,3)*kpg(3)
206        rr=sqrt(xx**2+yy**2+zz**2)
207 
208 !      Continue only for k+G<>0
209        if (rr>tol) then
210 
211 !        Determine theta and phi
212          cphi=one
213          sphi=zero
214          ctheta=zz/rr
215          stheta=sqrt(abs((one-ctheta)*(one+ctheta)))
216          if (stheta>tol) then
217            cphi=xx/(rr*stheta)
218            sphi=yy/(rr*stheta)
219          end if
220          do mm=1,mpsang-1
221            rphase(mm)=dreal(dcmplx(cphi,sphi)**mm)
222            iphase(mm)=aimag(dcmplx(cphi,sphi)**mm)
223          end do
224 
225 !        Determine gradients of theta and phi
226          if (optder/=0) then
227            dtheta(1)=ctheta*cphi
228            dtheta(2)=ctheta*sphi
229            dtheta(3)=-stheta
230            dphi(1)=-sphi
231            dphi(2)=cphi
232            dphi(3)=zero
233          end if
234 
235 !        COMPUTE Ylm(K)
236 !        ============================================
237          if (optder>=0) then
238 !          Loop over angular momentum l
239            do ilang=2,mpsang
240              ll=ilang-1
241              l0=ll**2+ll+1
242              fact=1._dp/real(ll*(ll+1),dp)
243              ylmcst=sqrt(real(2*ll+1,dp)/four_pi)
244 !            Special case m=0
245              ylm(ikg+ipw,l0)=ylmcst*ass_leg_pol(ll,0,ctheta)
246 !            Compute for m>0
247              onem=one
248              do mm=1,ll
249                onem=-onem
250                work1=ylmcst*sqrt(fact)*onem*ass_leg_pol(ll,mm,ctheta)*sqrt(2._dp)
251                ylm(ikg+ipw,l0+mm)=work1*rphase(mm)
252                ylm(ikg+ipw,l0-mm)=work1*iphase(mm)
253                if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp)
254              end do ! End loop over m
255            end do  ! End loop over l
256          end if
257 
258 !        COMPUTE dYlm/dKi
259 !        ============================================
260          if (optder/=0) then
261 !          Loop over angular momentum l
262            do ilang=2,mpsang
263              ll=ilang-1
264              l0=ll**2+ll+1
265              fact=1._dp/real(ll*(ll+1),dp)
266              ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/rr
267 !            === Special case m=0 ===
268 !            1-compute gradients in cartesian coordinates
269              work1=ylmcst*plm_dtheta(ll,0,ctheta)
270              ylmgr_cart(1:3,1)=work1*dtheta(1:3)
271 !            2-Transfer gradients into reduced coordinates
272              do ii=1,3
273                ylmgr_red(ii,1)=(rprimd(1,ii)*ylmgr_cart(1,1)+&
274 &               rprimd(2,ii)*ylmgr_cart(2,1)+&
275 &               rprimd(3,ii)*ylmgr_cart(3,1))
276              end do
277 !            3-Store gradients
278              ylm_gr(ikg+ipw,1:3,l0) =ylmgr_red(1:3,1)
279 !            === Compute for m>0 ===
280              onem=one
281              do mm=1,ll
282                onem=-onem
283 !              1-compute gradients in cartesian coordinates
284                work1=ylmcst*sqrt(fact)*onem*plm_dtheta(ll,mm,ctheta)*sqrt(2._dp)
285                work2=ylmcst*sqrt(fact)*onem*plm_dphi  (ll,mm,ctheta)*sqrt(2._dp)
286                ylmgr_cart(1:3,1)=rphase(mm)*work1*dtheta(1:3)-iphase(mm)*work2*dphi(1:3)
287                ylmgr_cart(1:3,2)=iphase(mm)*work1*dtheta(1:3)+rphase(mm)*work2*dphi(1:3)
288 !              2-Transfer gradients into reduced coordinates
289                do kk=1,2
290                  do ii=1,3
291                    ylmgr_red(ii,kk)=(rprimd(1,ii)*ylmgr_cart(1,kk)+&
292 &                   rprimd(2,ii)*ylmgr_cart(2,kk)+&
293 &                   rprimd(3,ii)*ylmgr_cart(3,kk))
294                  end do
295                end do
296 !              3-Store gradients
297                ylm_gr(ikg+ipw,1:3,l0+mm) =ylmgr_red(1:3,1)
298                ylm_gr(ikg+ipw,1:3,l0-mm) =ylmgr_red(1:3,2)
299                if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp)
300              end do ! End loop over m
301            end do  ! End loop over l
302          end if
303 
304 !        COMPUTE d2Ylm/dKidKj
305 !        ============================================
306          if (optder==2) then
307            call plm_coeff(blm,mpsang,ctheta)
308 !          Loop over angular momentum l
309            do ilang=2,mpsang
310              ll=ilang-1
311              l0=ll**2+ll+1
312              fact=1._dp/real(ll*(ll+1),dp)
313              ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/(rr**2)
314 !            === Special case m=0 ===
315 !            1-compute gradients in cartesian coordinates
316              ylmgr2_cart(1,1,1)=ylmcst*(-blm(3,l0)*sphi*sphi+blm(4,l0)*cphi*cphi)
317              ylmgr2_cart(2,2,1)=ylmcst*(-blm(3,l0)*cphi*cphi+blm(4,l0)*sphi*sphi)
318              ylmgr2_cart(3,3,1)=ylmcst*blm(1,l0)
319              ylmgr2_cart(3,1,1)=ylmcst*blm(2,l0)*cphi
320              ylmgr2_cart(3,2,1)=ylmcst*blm(2,l0)*sphi
321              ylmgr2_cart(2,1,1)=ylmcst*(blm(3,l0)+blm(4,l0))*sphi*cphi
322              ylmgr2_cart(1,3,1)=ylmgr2_cart(3,1,1)
323              ylmgr2_cart(1,2,1)=ylmgr2_cart(2,1,1)
324              ylmgr2_cart(2,3,1)=ylmgr2_cart(3,2,1)
325 !            2-Transfer gradients into reduced coordinates
326              do jj=1,3
327                do ii=1,3
328                  ylmgr2_tmp(ii,jj)=(rprimd(1,jj)*ylmgr2_cart(1,ii,1)+&
329 &                 rprimd(2,jj)*ylmgr2_cart(2,ii,1)+&
330 &                 rprimd(3,jj)*ylmgr2_cart(3,ii,1))
331                end do
332              end do
333              do ii=1,6
334                ia=alpha(ii);ib=beta(ii)
335                ylmgr_red(ii,1)=(rprimd(1,ia)*ylmgr2_tmp(1,ib)+&
336 &               rprimd(2,ia)*ylmgr2_tmp(2,ib)+&
337 &               rprimd(3,ia)*ylmgr2_tmp(3,ib))
338              end do
339              ylm_gr(ikg+ipw,4:9,l0) =ylmgr_red(1:6,1)
340 !            === Compute for m>0 ===
341              onem=one
342              do mm=1,ll
343                onem=-onem;ylmcst2=ylmcst*sqrt(fact)*sqrt(two)
344                ylmgr2_cart(1,1,1)=ylmcst2*((-blm(3,l0+mm)*sphi*sphi+blm(4,l0+mm)*cphi*cphi)*rphase(mm)-&
345 &               blm(5,l0+mm)*2.d0*cphi*sphi*mm*iphase(mm))
346                ylmgr2_cart(1,1,2)=ylmcst2*((-blm(3,l0+mm)*sphi*sphi+blm(4,l0+mm)*cphi*cphi)*iphase(mm)+&
347 &               blm(5,l0+mm)*2.d0*cphi*sphi*mm*rphase(mm))
348                ylmgr2_cart(2,2,1)=ylmcst2*((-blm(3,l0+mm)*cphi*cphi+blm(4,l0+mm)*sphi*sphi)*rphase(mm)+&
349 &               blm(5,l0+mm)*2.d0*cphi*sphi*mm*iphase(mm))
350                ylmgr2_cart(2,2,2)=ylmcst2*((-blm(3,l0+mm)*cphi*cphi+blm(4,l0+mm)*sphi*sphi)*iphase(mm)-&
351 &               blm(5,l0+mm)*2.d0*cphi*sphi*mm*rphase(mm))
352                ylmgr2_cart(3,3,1)=ylmcst2*blm(1,l0+mm)*rphase(mm)
353                ylmgr2_cart(3,3,2)=ylmcst2*blm(1,l0+mm)*iphase(mm)
354                ylmgr2_cart(3,1,1)=ylmcst2*(blm(2,l0+mm)*cphi*rphase(mm)-&
355 &               mm*iphase(mm)*sphi*onem*plm_dtheta(ll,mm,ctheta))
356                ylmgr2_cart(3,1,2)=ylmcst2*(blm(2,l0+mm)*cphi*iphase(mm)+&
357 &               mm*rphase(mm)*sphi*onem*plm_dtheta(ll,mm,ctheta))
358                ylmgr2_cart(3,2,1)=ylmcst2*(blm(2,l0+mm)*sphi*rphase(mm)+&
359 &               mm*iphase(mm)*cphi*onem*plm_dtheta(ll,mm,ctheta))
360                ylmgr2_cart(3,2,2)=ylmcst2*(blm(2,l0+mm)*sphi*iphase(mm)-&
361 &               mm*rphase(mm)*cphi*onem*plm_dtheta(ll,mm,ctheta))
362                ylmgr2_cart(2,1,1)=ylmcst2*((blm(3,l0+mm)+blm(4,l0+mm))*sphi*cphi*rphase(mm)-&
363 &               blm(5,l0+mm)*(sphi*sphi-cphi*cphi)*mm*iphase(mm))
364                ylmgr2_cart(2,1,2)=ylmcst2*((blm(3,l0+mm)+blm(4,l0+mm))*sphi*cphi*iphase(mm)+&
365 &               blm(5,l0+mm)*(sphi*sphi-cphi*cphi)*mm*rphase(mm))
366                ylmgr2_cart(1,3,:)=ylmgr2_cart(3,1,:)
367                ylmgr2_cart(1,2,:)=ylmgr2_cart(2,1,:)
368                ylmgr2_cart(2,3,:)=ylmgr2_cart(3,2,:)
369 !              2-Transfer gradients into reduced coordinates
370                do kk=1,2
371                  do jj=1,3
372                    do ii=1,3
373                      ylmgr2_tmp(ii,jj)=(rprimd(1,jj)*ylmgr2_cart(1,ii,kk)+&
374 &                     rprimd(2,jj)*ylmgr2_cart(2,ii,kk)+&
375 &                     rprimd(3,jj)*ylmgr2_cart(3,ii,kk))
376                    end do
377                  end do
378                  do ii=1,6
379                    ia=alpha(ii);ib=beta(ii)
380                    ylmgr_red(ii,kk)=(rprimd(1,ia)*ylmgr2_tmp(1,ib)+&
381 &                   rprimd(2,ia)*ylmgr2_tmp(2,ib)+&
382 &                   rprimd(3,ia)*ylmgr2_tmp(3,ib))
383                  end do
384                end do
385                ylm_gr(ikg+ipw,4:9,l0+mm) =ylmgr_red(1:6,1)
386                ylm_gr(ikg+ipw,4:9,l0-mm) =ylmgr_red(1:6,2)
387                if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp)
388              end do ! End loop over m
389            end do  ! End loop over l
390          end if
391 
392 !        End condition r<>0
393        end if
394 
395 !      End loop over k+G
396      end do
397 
398 !    End condition l<>0
399    end if
400 
401    ABI_DEALLOCATE(kg_k)
402 
403    ikg=ikg+npw_k
404  end do !  End Loop over k-points
405 
406 !Release the temporary memory
407 !Allocate some memory
408  if (optder/=0) then
409    ABI_DEALLOCATE(ylmgr_cart)
410  end if
411  if (optder/=0.and.optder/=2) then
412    ABI_DEALLOCATE(ylmgr_red)
413  end if
414  if (optder==2) then
415    ABI_DEALLOCATE(ylmgr2_cart)
416    ABI_DEALLOCATE(ylmgr2_tmp)
417    ABI_DEALLOCATE(ylmgr_red)
418    ABI_DEALLOCATE(blm)
419  end if
420 
421 end subroutine initylmg

ABINIT/m_initylmg [ Modules ]

[ Top ] [ 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-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 .

PARENTS

CHILDREN

SOURCE

22 #if defined HAVE_CONFIG_H
23 #include "config.h"
24 #endif
25 
26 #include "abi_common.h"
27 
28 module m_initylmg
29 
30  use defs_basis
31  use defs_abitypes
32  use m_abicore
33  use m_errors
34  use m_xmpi
35 
36  use m_paw_sphharm, only : ass_leg_pol, plm_dtheta, plm_dphi, plm_coeff
37  use m_mpinfo,      only : proc_distrb_cycle
38 
39  implicit none
40 
41  private