TABLE OF CONTENTS


ABINIT/m_paw_sphharm [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_sphharm

FUNCTION

  This module contains a set of routines to compute the complex (resp. real)
  spherical harmonics Ylm (resp. Slm) (and gradients).

COPYRIGHT

 Copyright (C) 2013-2022 ABINIT group (MT, FJ, NH, TRangel)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

NOTES

  FOR DEVELOPPERS: in order to preserve the portability of libPAW library,
  please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt

SOURCE

21 #include "libpaw.h"
22 
23 MODULE m_paw_sphharm
24 
25  USE_DEFS
26  USE_MSG_HANDLING
27  USE_MEMORY_PROFILING
28 
29  implicit none
30 
31  private
32 
33 !Public procedures.
34  public :: ylmc             ! Complex Spherical harmonics for l<=3.
35  public :: ylmcd            ! First derivative of complex Ylm wrt theta and phi up to l<=3
36  public :: ylm_cmplx        ! All (complex) spherical harmonics for lx<=4
37  public :: initylmr         ! Real Spherical Harmonics on a set of vectors
38  public :: ys               ! Matrix element <Yl'm'|Slm>
39  public :: lxyz             ! Matrix element <Yl'm'|L_idir|Ylm>
40  public :: slxyzs           ! Matrix element <Sl'm'|L_idir|Slm>
41  public :: lsylm            ! Compute the LS operator in the real spherical harmonics basis
42  public :: plm_coeff        ! Coefficients depending on Plm used to compute the 2nd der of Ylm
43  public :: ass_leg_pol      ! Associated Legendre Polynomial Plm(x)
44  public :: plm_dphi         ! m*P_lm(x)/sqrt((1-x^2)  (P_lm= associatedLegendre polynomial)
45  public :: plm_dtheta       ! -(1-x^2)^1/2*d/dx{P_lm(x)} (P_lm= associated Legendre polynomial)
46  public :: plm_d2theta      ! d2(Plm (cos(theta)))/d(theta)2 (P_lm= associated Legendre polynomial)
47  public :: pl_deriv         ! d2(Pl (x)))/d(x)2  where P_l is a Legendre polynomial
48  public :: ylm_angular_mesh ! Build (theta, phi) angular mesh
49  public :: mat_mlms2jmj     ! Change a matrix from the Ylm basis to the J,M_J basis
50  public :: mat_slm2ylm      ! Change a matrix from the Slm to the Ylm basis or from Ylm to Slm
51  public :: setsym_ylm       ! Compute rotation matrices expressed in the basis of real spherical harmonics
52  public :: setnabla_ylm     ! Evaluate several integrals involving spherical harmonics and their gradient
53  public :: gaunt            ! Gaunt coeffients for complex Yml
54  public :: realgaunt        ! Compute "real Gaunt coefficients" with "real spherical harmonics"
55  public :: nablarealgaunt   ! Compute the integrals Grad(Slimi).Grad(Sjmj) Slkmk of real spherical harmonics
56 
57 !Private functions
58  private :: create_slm2ylm  ! For a given angular momentum lcor, compute slm2ylm
59  private :: create_mlms2jmj ! For a given angular momentum lcor, give the rotation matrix msml2jmj
60  private :: mkeuler         ! For a given symmetry operation, determines the corresponding Euler angles
61  private :: dbeta           ! Calculate the rotation matrix d^l_{m{\prim}m}(beta)
62  private :: phim            ! Computes Phi_m[t]=Sqrt2.cos[m.t] (m>0), Sqrt2.sin[|m|.t] (m<0), 1 (m=0)
63  private :: gauleg          ! Compute the coefficients for Gauss-Legendre integration
64  private :: perms           ! Returns N!/(N-k)! if N>=0 and N>k
65  private :: rfactorial      ! Calculates N! as a double precision real

m_paw_sphharm/ass_leg_pol [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 ass_leg_pol

FUNCTION

 Compute the associated Legendre Polynomial Plm(x),
 using a stable recursion formula.
 Here m and l are integers satisfying 0<=m<=l,
 while x lies in the range -1<=x<=1

INPUTS

  l,m= l,m numbers
  xarg=argument of the polynom

OUTPUT

SOURCE

1202 function ass_leg_pol(l,m,xarg)
1203 
1204 !Arguments ------------------------------------
1205 !scalars
1206  integer, intent(in) ::  l,m
1207  real(dp), intent(in) :: xarg
1208  real(dp) :: ass_leg_pol
1209 
1210 !Local variables-------------------------------
1211 !scalars
1212  integer :: i,ll
1213  real(dp) :: pll,polmm,tmp1,sqrx,x
1214  character(len=100) :: msg
1215 
1216 ! *************************************************************************
1217 
1218  x=xarg
1219  if (m.lt.0.or.m.gt.l.or.abs(x).gt.1.d0) then
1220    if (m.lt.0.or.m.gt.l.or.abs(x).gt.1.d0+1.d-10) then
1221     msg='Bad choice of l, m or x !'
1222     LIBPAW_BUG(msg)
1223    endif
1224    x=1.d0
1225  endif
1226 
1227  polmm=1.d0
1228  if (m>0) then
1229   sqrx=sqrt(abs((1.d0-x)*(1.d0+x)))
1230   do i=1,m
1231    polmm=polmm*(1.0d0-2.0d0*i)*sqrx
1232   enddo
1233  endif
1234 
1235  if (l==m) then
1236   ass_leg_pol=polmm
1237  else
1238   tmp1=x*(2.0d0*m+1.0d0)*polmm
1239   if (l==(m+1)) then
1240    ass_leg_pol=tmp1
1241   else
1242    do ll=m+2,l
1243     pll=(x*(2.0d0*ll-1.0d0)*tmp1-(ll+m-1.0d0)*polmm)/dble(ll-m)
1244     polmm=tmp1
1245     tmp1=pll
1246    enddo
1247    ass_leg_pol=pll
1248   endif
1249  endif
1250 
1251 end function ass_leg_pol

m_paw_sphharm/create_mlms2jmj [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 create_mlms2jmj

FUNCTION

 For a given angular momentum lcor, give the rotation matrix msml2jmj

INPUTS

  lcor= angular momentum

SIDE EFFECTS

  mlms2jmj= rotation matrix

SOURCE

3179 subroutine create_mlms2jmj(lcor,mlmstwojmj)
3180 
3181 !Arguments ---------------------------------------------
3182 !scalars
3183  integer,intent(in) :: lcor
3184 !arrays
3185  complex(dpc),intent(out) :: mlmstwojmj(2*(2*lcor+1),2*(2*lcor+1))
3186 
3187 !Local variables ---------------------------------------
3188 !scalars
3189  integer :: jc1,jj,jm,ll,ml1,ms1
3190  real(dp) :: invsqrt2lp1,xj,xmj
3191  character(len=500) :: msg
3192 !arrays
3193  integer, allocatable :: ind_msml(:,:)
3194  complex(dpc),allocatable :: mat_mlms2(:,:)
3195 
3196 !*********************************************************************
3197 
3198 !--------------- Built indices + allocations
3199  ll=lcor
3200  mlmstwojmj=czero
3201  LIBPAW_BOUND2_ALLOCATE(ind_msml,BOUNDS(1,2),BOUNDS(-ll,ll))
3202  LIBPAW_ALLOCATE(mat_mlms2,(2*(2*lcor+1),2*(2*lcor+1)))
3203  mlmstwojmj=czero
3204  jc1=0
3205  do ms1=1,2
3206    do ml1=-ll,ll
3207      jc1=jc1+1
3208      ind_msml(ms1,ml1)=jc1
3209    end do
3210  end do
3211 
3212 !--------------- built mlmstwojmj
3213 !do jj=ll,ll+1    ! the physical value of j are ll-0.5,ll+0.5
3214 !xj(jj)=jj-0.5
3215  if(ll==0)then
3216    msg=' ll should not be equal to zero !'
3217    LIBPAW_BUG(msg)
3218  end if
3219  jc1=0
3220  invsqrt2lp1=one/sqrt(float(2*lcor+1))
3221  do jj=ll,ll+1
3222    xj=float(jj)-half
3223    do jm=-jj,jj-1
3224      xmj=float(jm)+half
3225      jc1=jc1+1
3226      if(nint(xj+0.5)==ll+1) then
3227        if(nint(xmj+0.5)==ll+1)  then
3228          mlmstwojmj(ind_msml(2,ll),jc1)=1.0   !  J=L+0.5 and m_J=L+0.5
3229        else if(nint(xmj-0.5)==-ll-1) then
3230          mlmstwojmj(ind_msml(1,-ll),jc1)=1.0   !  J=L+0.5 and m_J=-L-0.5
3231        else
3232          mlmstwojmj(ind_msml(2,nint(xmj-0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5))
3233          mlmstwojmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)-xmj+0.5))
3234        end if
3235      end if
3236      if(nint(xj+0.5)==ll) then
3237        mlmstwojmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5))
3238        mlmstwojmj(ind_msml(2,nint(xmj-0.5)),jc1)=-invsqrt2lp1*(sqrt(float(ll)-xmj+0.5))
3239      end if
3240    end do
3241  end do
3242 
3243  LIBPAW_DEALLOCATE(ind_msml)
3244  LIBPAW_DEALLOCATE(mat_mlms2)
3245 
3246 end subroutine create_mlms2jmj

m_paw_sphharm/create_slm2ylm [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 create_slm2ylm

FUNCTION

 For a given angular momentum lcor, compute slm2ylm.

INPUTS

  lcor= angular momentum, size of the matrix is 2(2*lcor+1)

OUTPUT

  slm2ylm(2lcor+1,2lcor+1) = rotation matrix.

NOTES

  useful only in ndij==4

SOURCE

3124 subroutine create_slm2ylm(lcor,slmtwoylm)
3125 
3126 !Arguments ---------------------------------------------
3127 !scalars
3128  integer,intent(in) :: lcor
3129 !arrays
3130  complex(dpc),intent(out) :: slmtwoylm(2*lcor+1,2*lcor+1)
3131 
3132 !Local variables ---------------------------------------
3133 !scalars
3134  integer :: jm,ll,mm,im
3135  real(dp),parameter :: invsqrt2=one/sqrt2
3136  real(dp) :: onem
3137 !arrays
3138 
3139 ! *********************************************************************
3140 
3141  ll=lcor
3142  slmtwoylm=czero
3143  do im=1,2*ll+1
3144    mm=im-ll-1;jm=-mm+ll+1
3145    onem=dble((-1)**mm)
3146    if (mm> 0) then
3147      slmtwoylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp)
3148      slmtwoylm(jm,im)= cmplx(invsqrt2,     zero,kind=dp)
3149    end if
3150    if (mm==0) then
3151      slmtwoylm(im,im)=cone
3152    end if
3153    if (mm< 0) then
3154      slmtwoylm(im,im)= cmplx(zero,     invsqrt2,kind=dp)
3155      slmtwoylm(jm,im)=-cmplx(zero,onem*invsqrt2,kind=dp)
3156    end if
3157  end do
3158 
3159 end subroutine create_slm2ylm

m_paw_sphharm/dbeta [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 dbeta

FUNCTION

  Calculate the rotation matrix d^l_{m{\prim}m}(beta) using Eq. 4.14 of
  M.E. Rose, Elementary Theory of Angular Momentum,
             John Wiley & Sons, New-York, 1957

INPUTS

  cosbeta= cosinus of beta (=Euler angle)
  sinbeta= sinus of beta (=Euler angle)
  ll= index l
  mm= index m
  mp= index m_prime

OUTPUT

  dbeta= rotation matrix

NOTES

  - This file comes from the file crystal_symmetry.f
    by N.A.W. Holzwarth and A. Tackett for the code pwpaw
  - Assume l relatively small so that factorials do not cause
    roundoff error
  - XG20200718 This routine was inaccurate when cosbeta was close to one or minus one. 
    This has been fixed by adding sinbeta argument obtained from mkeuler. 
    Tolerances have been adjusted as well.

SOURCE

3451 function dbeta(cosbeta,sinbeta,ll,mp,mm)
3452 
3453 !Arguments ---------------------------------------------
3454 !scalars
3455  integer,intent(in) :: ll,mm,mp
3456  real(dp) :: dbeta
3457  real(dp),intent(in) :: cosbeta,sinbeta
3458 
3459 !Local variables ------------------------------
3460 !scalars
3461  integer,parameter :: mxterms=200
3462  integer :: ii,ina,inb,inc,ml,ms
3463  real(dp) :: arg,cosbetab2,pref,sinbetab2,sum,tt
3464 
3465 !************************************************************************
3466  dbeta=zero
3467 
3468 !Special cases
3469  if (abs(cosbeta-1._dp).lt.tol10) then
3470    if (mp.eq.mm) dbeta=1
3471  else if (abs(cosbeta+1._dp).lt.tol10) then
3472    if (mp.eq.-mm) dbeta=(-1)**(ll+mm)
3473  else
3474 !  General case
3475 
3476 !!!!! Old coding
3477 !!  This is inaccurate when cosbeta is close to -1
3478 !   cosbetab2=sqrt((1+cosbeta)*0.5_dp)
3479 !!  This is inaccurate when cosbeta is close to +1
3480 !   sinbetab2=sqrt((1-cosbeta)*0.5_dp)
3481 !!!!! End old coding, begin new coding
3482   if(cosbeta>-tol8)then
3483     !If cosbeta is positive, cosbeta2 is positive with value >0.7, so one can divide by cosbetab2
3484     cosbetab2=sqrt((1+cosbeta)*half)
3485     sinbetab2=sinbeta*half/cosbetab2
3486   else
3487     !If cosbeta is negative, sinbeta2 is positive with value >0.7, so one can divide by sinbetab2
3488     sinbetab2=sqrt((1-cosbeta)*half)
3489     cosbetab2=sinbeta*half/sinbetab2
3490   endif
3491 !!!!! End of new coding
3492 
3493    ml=max(mp,mm)
3494    ms=min(mp,mm)
3495    if (ml.ne.mp) sinbetab2=-sinbetab2
3496    tt=-(sinbetab2/cosbetab2)**2
3497    pref=sqrt((rfactorial(ll-ms)*rfactorial(ll+ml))&
3498 &   /(rfactorial(ll+ms)*rfactorial(ll-ml)))&
3499 &   /rfactorial(ml-ms)*(cosbetab2**(2*ll+ms-ml))&
3500 &   *((-sinbetab2)**(ml-ms))
3501    sum=1._dp
3502    arg=1._dp
3503    ina=ml-ll
3504    inb=-ms-ll
3505    inc=ml-ms+1
3506    do ii=1,mxterms
3507      if (ina.eq.0.or.inb.eq.0) exit
3508      arg=(arg*ina*inb*tt)/(ii*inc)
3509      sum=sum+arg
3510      ina=ina+1
3511      inb=inb+1
3512      inc=inc+1
3513    end do
3514    dbeta=pref*sum
3515  end if
3516 
3517 end function dbeta

m_paw_sphharm/gauleg [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 gauleg

FUNCTION

 Private function
 Compute the coefficients (supports and weights) for Gauss-Legendre integration

INPUTS

  xmin=lower bound of integration
  xmax=upper bound of integration
  nn=order of integration

OUTPUT

  x(nn)=array of support points
  weights(n)=array of integration weights

SOURCE

3541  subroutine gauleg(xmin,xmax,x,weights,nn)
3542 
3543 !Arguments ---------------------------------------------
3544 !scalars
3545  integer,intent(in) :: nn
3546  real(dp),intent(in) :: xmax,xmin
3547 !arrays
3548  real(dp),intent(out) :: weights(nn),x(nn)
3549 
3550 !Local variables ------------------------------
3551 !scalars
3552  integer :: ii,jj
3553  real(dp),parameter :: tol=1.d-13
3554  real(dp) :: p1,p2,p3,pi,xl,pp,xmean,z,z1
3555 
3556 !************************************************************************
3557 
3558  pi=4._dp*atan(1._dp)
3559  xl=(xmax-xmin)*0.5_dp
3560  xmean=(xmax+xmin)*0.5_dp
3561 
3562  do ii=1,(nn+1)/2
3563    z=cos(pi*(ii-0.25_dp)/(nn+0.5_dp))
3564    do
3565      p1=1._dp
3566      p2=0._dp
3567      do jj=1,nn
3568        p3=p2
3569        p2=p1
3570        p1=((2._dp*jj-1._dp)*z*p2-(jj-1._dp)*p3)/jj
3571      end do
3572      pp=nn*(p2-z*p1)/(1._dp-z**2)
3573      z1=z
3574      z=z1-p1/pp
3575      if(abs(z-z1) < tol) exit
3576    end do
3577    x(ii)=xmean-xl*z
3578    x(nn+1-ii)=xmean+xl*z
3579    weights(ii)=2._dp*xl/((1._dp-z**2)*pp**2)
3580    weights(nn+1-ii)=weights(ii)
3581  end do
3582 
3583  end subroutine gauleg

m_paw_sphharm/gaunt [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 gaunt

FUNCTION

 Returns gaunt coefficient, i.e.
   the integral of Sqrt[4 \pi] Y*(l_i,m_i) Y*(ll,mm) Y(l_j,m_j)
   See the 3-j and 6-j symbols by Rotenberg, etc., (Technology Press, 1959), pg.5.

INPUTS

   ll,mm,l1,l2,m1,m2= six quantum numbers defining the Gaunt coef.

OUTPUT

   gaunt(ll,mm,l1,l2,m1,m2)=the value of the integral

SOURCE

2648 function gaunt(ll,mm,l1,m1,l2,m2)
2649 
2650 !Arguments ---------------------------------------------
2651 !scalars
2652  integer,intent(in) :: l1,l2,ll,m1,m2,mm
2653  real(dp) :: gaunt
2654 
2655 !Local variables ------------------------------
2656 !scalars
2657  integer :: i1,i2,j1,j1half,j2,j2half,j3,j3half,j_half,jj,k1,k2,n1,n2
2658  real(dp) :: argument,sign,sum,xx,yy
2659  logical :: ok
2660 
2661 !************************************************************************
2662 
2663  gaunt=zero;sum=zero;ok =.true.
2664 
2665  if((-m1-mm+m2) /= 0) ok = .false.
2666  if(abs(m1) > l1) ok = .false.
2667  if(abs(mm) > ll) ok = .false.
2668  if(abs(m2) > l2) ok = .false.
2669 
2670  jj = l1 + ll + l2
2671  if (mod(jj,2)/=0) ok = .false.
2672  j1 = jj-2*l2
2673  j2 = jj-2*ll
2674  j3 = jj-2*l1
2675 
2676  if (j1<0 .or. j2<0 .or. j3<0) ok = .false.
2677 
2678  if (ok) then
2679 
2680    xx = (2 * l1 + 1) * (2 * ll + 1) * (2 * l2 + 1)
2681 
2682    j1half = j1/2
2683    j2half = j2/2
2684    j3half = j3/2
2685    j_half = jj/2
2686 
2687    gaunt = (-1)**j1half * sqrt(xx)
2688    gaunt = gaunt * rfactorial(j2)*rfactorial(j3)/rfactorial(jj+1)
2689    gaunt = gaunt * rfactorial(j_half)/(rfactorial(j1half)&
2690 &                * rfactorial(j2half)*rfactorial(j3half))
2691 
2692    yy = rfactorial(l2 + m2) * rfactorial(l2 - m2)
2693 
2694    if (mm>=0) then
2695      yy = yy * perms(ll+mm,2*mm)
2696    else
2697      yy = yy / perms(ll-mm,-2*mm)
2698    end if
2699 
2700    if (m1>=0) then
2701      yy = yy / perms(l1+m1,2*m1)
2702    else
2703      yy = yy * perms(l1-m1,-2*m1)
2704    end if
2705 
2706    gaunt = gaunt * sqrt(yy)
2707 
2708    i1 = l2 - ll - m1
2709    i2 = l2 - l1 + mm
2710    k1 = -min(0, i1, i2)
2711    n1 = l1 + m1
2712    n2 = ll - mm
2713    k2 = min(j1, n1, n2)
2714 
2715    sign = 1._dp
2716    if(k1>0) sign = (-1._dp)**k1
2717 
2718    argument = sign     * perms(n1,k1)/rfactorial(k1)
2719    argument = argument * perms(n2,k1)/rfactorial(i1 + k1)
2720    argument = argument * perms(j1,k1)/rfactorial(i2 + k1)
2721    sum = sum + argument
2722 
2723    sign = -sign
2724    k1 = k1 + 1
2725    do while(k1 <= k2)
2726      argument = sign     * perms(n1, k1)/rfactorial(k1)
2727      argument = argument * perms(n2, k1)/rfactorial(i1 + k1)
2728      argument = argument * perms(j1, k1)/rfactorial(i2 + k1)
2729      sum = sum + argument
2730      sign = -sign
2731      k1 = k1 + 1
2732    end do
2733 
2734  end if
2735 
2736  gaunt = gaunt * sum
2737 
2738  end function gaunt

m_paw_sphharm/initylmr [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 initylmr

FUNCTION

 Calculate the real spherical harmonics Ylm (and gradients)
 over a set of (r) vectors given in Cartesian coordinates.

INPUTS

  mpsang= 1+maximum angular momentum for nonlocal pseudopotential
  normchoice=0  the input rr vectors are normalized
            =1  the norm of the input vector is in nrm() array
  nrm(npts) = Depending of normchoice, this array contains
              either the weight of the point or the norm of rr.
  npts = number of rr vectors
  option= 1=compute Ylm(R), 2=compute Ylm(R) and dYlm/dRi (cartesian derivatives),
          3=compute Ylm(R), dYlm/dRi and d2Ylm/dRidRj (cartesian derivatives)
  rr(3,npts)=  vectors for which ylmr have to be calculated
               For each point of the spherical mesh, gives the
               Cartesian coordinates of the corresponding point.

OUTPUT

  if (option=1, 2 or 3)
    ylm(mpsang*mpsang,npts)     = real spherical harmonics for each r point
  if (option=2 or 3)
    ylmr_gr(1:3,mpsang*mpsang,npts)= gradients of real spherical harmonics
  if (option=3)
    ylmr_gr(4:9,mpsang*mpsang,npts)= first and second gradients of real spherical harmonics

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 imaginary 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

536 subroutine initylmr(mpsang,normchoice,npts,nrm,option,rr,ylmr,ylmr_gr)
537 
538 !Arguments ------------------------------------
539 !scalars
540  integer,intent(in) :: mpsang,normchoice,npts,option
541 !arrays
542  real(dp),intent(in) :: nrm(npts),rr(3,npts)
543  real(dp),intent(out) :: ylmr(mpsang*mpsang,npts)
544  real(dp),optional,intent(out) :: ylmr_gr(3*(option/2)+6*(option/3),mpsang*mpsang,npts)
545 
546 !Local variables ------------------------------
547 !scalars
548  integer :: dimgr,ilang,inpt,l0,ll,mm
549  real(dp) :: cphi,ctheta,fact,onem,rnorm,sphi,stheta,work1,work2,ylmcst,ylmcst2
550  logical :: compute_ylm,compute_ylm2gr,compute_ylmgr
551 !arrays
552  real(dp) :: dphi(3),dtheta(3),iphase(mpsang-1),rphase(mpsang-1)
553  real(dp),allocatable :: blm(:,:)
554 
555 !************************************************************************
556 
557 !What has to be computed ?
558  compute_ylm   = (option==1.or.option==2.or.option==3)
559  compute_ylmgr =((             option==2.or.option==3).and.present(ylmr_gr))
560  compute_ylm2gr=((                          option==3).and.present(ylmr_gr))
561  dimgr=3*(option/2)+6*(option/3)
562 
563 !Initialisation of spherical harmonics
564  if (compute_ylm  ) ylmr   (:  ,1:npts)=zero
565  if (compute_ylmgr) ylmr_gr(:,:,1:npts)=zero
566 
567 !Special case for l=0
568  if (compute_ylm  ) ylmr(1,1:npts)=1._dp/sqrt(four_pi)
569  if (compute_ylmgr) ylmr_gr(1:dimgr,1,1:npts)=zero
570  if (mpsang>1) then
571 
572 !  Loop over all rr
573    do inpt=1,npts
574 
575 !    Load module of rr
576      rnorm=one
577      if (normchoice==1) rnorm=nrm(inpt)
578 
579 !    Continue only for r<>0
580 
581      if (rnorm>tol10) then
582 
583 !      Determine theta and phi
584        cphi=one
585        sphi=zero
586        ctheta=rr(3,inpt)/rnorm
587 !      MM030519 : abs is needed to prevent very small negative arg
588        stheta=sqrt(abs((one-ctheta)*(one+ctheta)))
589        if (stheta>tol10) then
590          cphi=rr(1,inpt)/(rnorm*stheta)
591          sphi=rr(2,inpt)/(rnorm*stheta)
592        end if
593        do mm=1,mpsang-1
594          rphase(mm)=dreal(dcmplx(cphi,sphi)**mm)
595          iphase(mm)=aimag(dcmplx(cphi,sphi)**mm)
596        end do
597 
598 !      Determine gradients of theta and phi
599        if (compute_ylmgr) then
600          dtheta(1)=ctheta*cphi
601          dtheta(2)=ctheta*sphi
602          dtheta(3)=-stheta
603          dphi(1)=-sphi
604          dphi(2)=cphi
605          dphi(3)=zero
606        end if
607 
608 !      COMPUTE Ylm(R)
609        if (compute_ylm) then
610 !        Loop over angular momentum l
611          do ilang=2,mpsang
612            ll=ilang-1
613            l0=ll**2+ll+1
614            fact=1._dp/real(ll*(ll+1),dp)
615            ylmcst=sqrt(real(2*ll+1,dp)/four_pi)
616 !          Special case m=0
617            ylmr(l0,inpt)=ylmcst*ass_leg_pol(ll,0,ctheta)
618 !          Compute for m>0
619            onem=one
620            do mm=1,ll
621              onem=-onem
622              work1=ylmcst*sqrt(fact)*onem*ass_leg_pol(ll,mm,ctheta)*sqrt(2._dp)
623              ylmr(l0+mm,inpt)=work1*rphase(mm)
624              ylmr(l0-mm,inpt)=work1*iphase(mm)
625              if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp)
626            end do ! End loop over m
627          end do  ! End loop over l
628        end if
629 
630 !      COMPUTE dYlm/dRi
631        if (compute_ylmgr) then
632 !        Loop over angular momentum l
633          do ilang=2,mpsang
634            ll=ilang-1
635            l0=ll**2+ll+1
636            fact=1._dp/real(ll*(ll+1),dp)
637            ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/rnorm
638 !          Special case m=0
639            work1=ylmcst*plm_dtheta(ll,0,ctheta)
640            ylmr_gr(1:3,l0,inpt)=work1*dtheta(1:3)
641 !          Compute for m>0
642            onem=one
643            do mm=1,ll
644              onem=-onem
645              work1=ylmcst*sqrt(fact)*onem*plm_dtheta(ll,mm,ctheta)*sqrt(2._dp)
646              work2=ylmcst*sqrt(fact)*onem*plm_dphi  (ll,mm,ctheta)*sqrt(2._dp)
647              ylmr_gr(1:3,l0+mm,inpt)=rphase(mm)*work1*dtheta(1:3)-iphase(mm)*work2*dphi(1:3)
648              ylmr_gr(1:3,l0-mm,inpt)=iphase(mm)*work1*dtheta(1:3)+rphase(mm)*work2*dphi(1:3)
649              if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp)
650            end do ! End loop over m
651          end do  ! End loop over l
652        end if
653 
654 !      COMPUTE d2Ylm/dRidRj
655        if (compute_ylm2gr) then
656          LIBPAW_ALLOCATE(blm,(5,mpsang*mpsang))
657          call plm_coeff(blm,mpsang,ctheta)
658 
659 !        Loop over angular momentum l
660          do ilang=2,mpsang
661            ll=ilang-1
662            l0=ll**2+ll+1
663            fact=1._dp/real(ll*(ll+1),dp)
664            ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/(rnorm**2)
665 !          Special case m=0
666            ylmr_gr(4,l0,inpt)=ylmcst*(-blm(3,l0)*sphi*sphi+blm(4,l0)*cphi*cphi)
667            ylmr_gr(5,l0,inpt)=ylmcst*(-blm(3,l0)*cphi*cphi+blm(4,l0)*sphi*sphi)
668            ylmr_gr(6,l0,inpt)=ylmcst*blm(1,l0)
669            ylmr_gr(7,l0,inpt)=ylmcst*blm(2,l0)*sphi
670            ylmr_gr(8,l0,inpt)=ylmcst*blm(2,l0)*cphi
671            ylmr_gr(9,l0,inpt)=ylmcst*(blm(3,l0)+blm(4,l0))*sphi*cphi
672 !          Compute for m>0
673            onem=one
674            do mm=1,ll
675              onem=-onem;ylmcst2=ylmcst*sqrt(fact)*sqrt(two)
676              ylmr_gr(4,l0+mm,inpt)=ylmcst2*((-blm(3,l0+mm)*sphi*sphi+blm(4,l0+mm)*cphi*cphi)*rphase(mm)-&
677 &             blm(5,l0+mm)*2.d0*cphi*sphi*mm*iphase(mm))
678              ylmr_gr(4,l0-mm,inpt)=ylmcst2*((-blm(3,l0+mm)*sphi*sphi+blm(4,l0+mm)*cphi*cphi)*iphase(mm)+&
679 &             blm(5,l0+mm)*2.d0*cphi*sphi*mm*rphase(mm))
680              ylmr_gr(5,l0+mm,inpt)=ylmcst2*((-blm(3,l0+mm)*cphi*cphi+blm(4,l0+mm)*sphi*sphi)*rphase(mm)+&
681 &             blm(5,l0+mm)*2.d0*cphi*sphi*mm*iphase(mm))
682              ylmr_gr(5,l0-mm,inpt)=ylmcst2*((-blm(3,l0+mm)*cphi*cphi+blm(4,l0+mm)*sphi*sphi)*iphase(mm)-&
683 &             blm(5,l0+mm)*2.d0*cphi*sphi*mm*rphase(mm))
684              ylmr_gr(6,l0+mm,inpt)=ylmcst2*blm(1,l0+mm)*rphase(mm)
685              ylmr_gr(6,l0-mm,inpt)=ylmcst2*blm(1,l0+mm)*iphase(mm)
686              ylmr_gr(7,l0+mm,inpt)=ylmcst2*(blm(2,l0+mm)*sphi*rphase(mm)+&
687 &             mm*iphase(mm)*cphi*onem*plm_dtheta(ll,mm,ctheta))
688              ylmr_gr(7,l0-mm,inpt)=ylmcst2*(blm(2,l0+mm)*sphi*iphase(mm)-&
689 &             mm*rphase(mm)*cphi*onem*plm_dtheta(ll,mm,ctheta))
690              ylmr_gr(8,l0+mm,inpt)=ylmcst2*(blm(2,l0+mm)*cphi*rphase(mm)-&
691 &             mm*iphase(mm)*sphi*onem*plm_dtheta(ll,mm,ctheta))
692              ylmr_gr(8,l0-mm,inpt)=ylmcst2*(blm(2,l0+mm)*cphi*iphase(mm)+&
693 &             mm*rphase(mm)*sphi*onem*plm_dtheta(ll,mm,ctheta))
694              ylmr_gr(9,l0+mm,inpt)=ylmcst2*((blm(3,l0+mm)+blm(4,l0+mm))*sphi*cphi*rphase(mm)-&
695 &             blm(5,l0+mm)*(sphi*sphi-cphi*cphi)*mm*iphase(mm))
696              ylmr_gr(9,l0-mm,inpt)=ylmcst2*((blm(3,l0+mm)+blm(4,l0+mm))*sphi*cphi*iphase(mm)+&
697 &             blm(5,l0+mm)*(sphi*sphi-cphi*cphi)*mm*rphase(mm))
698              if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp)
699            end do ! End loop over m
700          end do  ! End loop over l
701          LIBPAW_DEALLOCATE(blm)
702        end if
703 
704 !      End condition r<>0
705      end if
706 
707 !    End loop over rr
708    end do
709 
710 !  End condition l<>0
711  end if
712 
713 end subroutine initylmr

m_paw_sphharm/lsylm [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 lsylm

FUNCTION

 Compute the LS operator in the real spherical harmonics basis
 ls_ylm(ilm1,ilm2,ispin)= <sigma, S_lm1| L.S |S_lm2, sigma_prime>
   ilm,1m2=(l,m1,m2) with -l<=m1<=l, -l<=m2<=l and 0<l<=lmax
   ispin=(sigma,sigma_prime) 1=(up,up), 2=(up,dn), 3=(dn,up), 4=(dn,dn)

INPUTS

  lmax= max. value of angular momentum l

OUTPUT

  ls_ylm(2,l_max**2*(l_max**2+1)/2,2)=LS operator in the real spherical harmonics basis
        ls_ylm(:,:,1)=<up, S_lm1| L.S |S_lm2, up>
        ls_ylm(:,:,2)=<up, S_lm1| L.S |S_lm2, down>
        One can deduce:
        <down, S_lm1| L.S |S_lm2, down>=-<up, S_lm1| L.S |S_lm2, up>
        <down, S_lm1| L.S |S_lm2, up>  =-Conjg[<up, S_lm1| L.S |S_lm2, down>]
        Also, only ilm1<=ilm2 terms are stored, because:
         <sigma, S_lm1| L.S |S_lm2, sigma_prime>=-<sigma_prime, S_lm1| L.S |S_lm2, sigma>

SOURCE

 915 subroutine lsylm(ls_ylm,lmax)
 916 
 917 !Arguments ---------------------------------------------
 918 !scalars
 919  integer,intent(in) :: lmax
 920 !arrays
 921  real(dp),allocatable :: ls_ylm(:,:,:)
 922 
 923 !Local variables ---------------------------------------
 924 !scalars
 925  integer :: ii,ilm,im,j0lm,jj,jlm,jm,klm,ll,lm0,mm,ispden
 926  real(dp),parameter :: invsqrt2=one/sqrt2
 927  real(dp) :: onem
 928  character(len=500) :: msg
 929  logical,parameter :: tso=.false. ! use true to Test Spin Orbit and
 930 !                                   write the matrix of L.S in different basis
 931 !arrays
 932  complex(dpc) :: tmp(2)
 933  complex(dpc),allocatable :: ls_cplx(:,:,:),slm2ylm(:,:)
 934  complex(dpc),allocatable :: mat_inp_c(:,:,:),mat_out_c(:,:,:)
 935  complex(dpc),allocatable :: mat_ls_ylm(:,:,:),mat_jmj(:,:)
 936  character(len=9),parameter :: dspin2(2)=(/"up-up    ","up-dn    "/)
 937  character(len=9),parameter :: dspin6(6)=(/"dn       ","up       ","dn-dn    ","up-up    ","dn-up    ","up-dn    "/)
 938  character(len=9),parameter :: dspinm(6)=(/"dn       ","up       ","n        ","mx       ","my       ","mz       "/)
 939 
 940 ! *************************************************************************
 941 
 942  if (.not.allocated(ls_ylm)) then
 943    msg='ls_ylm is not allocated!'
 944    LIBPAW_BUG(msg)
 945  end if
 946  if ( size(ls_ylm) < 2*(lmax+1)**2 * ((lmax+1)**2+1) ) then
 947    msg='wrong size for ls_ylm!'
 948    LIBPAW_BUG(msg)
 949  end if
 950 
 951 !Initialization
 952  ls_ylm=zero
 953 
 954 !Nothing to do if lmax=0
 955  if (lmax<=0) return
 956 
 957 !Loop on l quantum number
 958  do ll=1,lmax
 959 
 960 !  Transformation matrixes: real->complex spherical harmonics
 961    LIBPAW_ALLOCATE(slm2ylm,(2*ll+1,2*ll+1))
 962    slm2ylm=czero
 963    do im=1,2*ll+1
 964      mm=im-ll-1;jm=-mm+ll+1
 965      onem=dble((-1)**mm)
 966      if (mm> 0) then
 967        slm2ylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp)
 968        slm2ylm(jm,im)= cmplx(invsqrt2,     zero,kind=dp)
 969      end if
 970      if (mm==0) then
 971        slm2ylm(im,im)=cone
 972      end if
 973      if (mm< 0) then
 974        slm2ylm(im,im)= cmplx(zero,     invsqrt2,kind=dp)
 975        slm2ylm(jm,im)=-cmplx(zero,onem*invsqrt2,kind=dp)
 976      end if
 977    end do
 978 
 979 !  Compute <sigma, Y_lm1|L.S|Y_lm2, sigma_prime> (Y_lm=complex spherical harmonics)
 980 !  1= <up|L.S|up>  ;  2= <up|L.S|dn>
 981    LIBPAW_ALLOCATE(ls_cplx,(2*ll+1,2*ll+1,2))
 982    ls_cplx=czero
 983    if(tso)  then
 984      LIBPAW_ALLOCATE(mat_ls_ylm,(2*ll+1,2*ll+1,4))
 985      if(tso) mat_ls_ylm=czero
 986    end if
 987    if(tso)  then
 988      LIBPAW_ALLOCATE(mat_jmj,(2*(2*ll+1),2*(2*ll+1)))
 989      if(tso) mat_jmj=czero
 990    end if
 991    do im=1,2*ll+1
 992      mm=im-ll-1
 993      ls_cplx(im,im,1)=half*mm
 994      if(tso) mat_ls_ylm(im,im,1)=-half*mm ! dn dn
 995      if(tso) mat_ls_ylm(im,im,2)=half*mm  ! up up
 996      if ((mm+1)<= ll) then
 997        ls_cplx(im,im+1,2)=half*sqrt(real((ll-mm)*(ll+mm+1),kind=dp))
 998        if(tso) mat_ls_ylm(im,im+1,4)=half*sqrt(real((ll-mm)*(ll+mm+1),kind=dp))  ! up dn
 999        if(tso) mat_ls_ylm(im+1,im,3)=half*sqrt(real((ll-mm)*(ll+mm+1),kind=dp))  ! dn up
1000      end if
1001      if ((mm-1)>=-ll) then
1002        ls_cplx(im-1,im,2)=half*sqrt(real((ll+mm)*(ll-mm+1),kind=dp))
1003        if(tso) mat_ls_ylm(im-1,im,4)=half*sqrt(real((ll+mm)*(ll-mm+1),kind=dp))  ! up dn
1004        if(tso) mat_ls_ylm(im,im-1,3)=half*sqrt(real((ll+mm)*(ll-mm+1),kind=dp))  ! dn up
1005      end if
1006    end do
1007 
1008 !  test : print LS in J,M_J basis
1009    if(tso) then
1010      do ispden=1,4
1011        write(msg,'(3a)') ch10,"value of LS in the Ylm basis for " ,trim(dspin6(ispden+2*(4/4)))
1012        call wrtout(std_out,msg,'COLL')
1013        do im=1,ll*2+1
1014          write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') (mat_ls_ylm(im,jm,ispden),jm=1,ll*2+1)
1015          call wrtout(std_out,msg,'COLL')
1016        end do
1017      end do
1018      call mat_mlms2jmj(ll,mat_ls_ylm,mat_jmj,4,1,2,3,std_out,'COLL')  ! optspin=2 : dn spin are first
1019    end if
1020 
1021 !  Compute <sigma, S_lm1|L.S|S_lm2, sigma_prime> (S_lm=real spherical harmonics)
1022 !  1= <up|L.S|up>  ;  2= <up|L.S|dn>
1023    if(tso) then
1024      LIBPAW_ALLOCATE(mat_inp_c,(2*ll+1,2*ll+1,4))
1025      LIBPAW_ALLOCATE(mat_out_c,(2*ll+1,2*ll+1,4))
1026    end if
1027    lm0=ll**2
1028    do jm=1,2*ll+1
1029      jlm=lm0+jm;j0lm=jlm*(jlm-1)/2
1030      do im=1,jm
1031        ilm=lm0+im;klm=j0lm+ilm
1032        tmp(:)=czero
1033        do ii=1,2*ll+1
1034          do jj=1,2*ll+1
1035            tmp(:)=tmp(:)+ls_cplx(ii,jj,:)*CONJG(slm2ylm(ii,im))*slm2ylm(jj,jm)
1036          end do
1037        end do
1038        ls_ylm(1,klm,:)=REAL(tmp(:),kind=dp)
1039        ls_ylm(2,klm,:)=AIMAG(tmp(:))
1040      end do
1041    end do
1042 
1043 !  Test: print LS in Slm basis
1044    if(tso) then
1045      call mat_slm2ylm(ll,mat_ls_ylm,mat_inp_c,4,2,2,3,std_out,'COLL') ! from Ylm to Slm, and dn spin are first
1046      do ispden=1,4
1047        write(msg,'(3a)') ch10,"value of LS in the Slm basis for " ,trim(dspin6(ispden+2*(4/4)))
1048        call wrtout(std_out,msg,'COLL')
1049        do im=1,ll*2+1
1050          write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') (mat_inp_c(im,jm,ispden),jm=1,ll*2+1)
1051          call wrtout(std_out,msg,'COLL')
1052        end do
1053      end do
1054 !    change into n,m basis
1055      mat_ls_ylm(:,:,1)=(mat_inp_c(:,:,1)+mat_inp_c(:,:,2))
1056      mat_ls_ylm(:,:,2)=(mat_inp_c(:,:,3)+mat_inp_c(:,:,4))
1057      mat_ls_ylm(:,:,3)=-cmplx(0.d0,1.d0)*(mat_inp_c(:,:,4)-mat_inp_c(:,:,3))
1058      mat_ls_ylm(:,:,4)=(mat_inp_c(:,:,1)-mat_inp_c(:,:,2))
1059      do ispden=1,4
1060        write(msg,'(3a)') ch10,"value of LS in the Slm basis for " ,trim(dspinm(ispden+2*(4/4)))
1061        call wrtout(std_out,msg,'COLL')
1062        do im=1,ll*2+1
1063          write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') (mat_ls_ylm(im,jm,ispden),jm=1,ll*2+1)
1064          call wrtout(std_out,msg,'COLL')
1065        end do
1066      end do
1067      LIBPAW_DEALLOCATE(mat_inp_c)
1068      LIBPAW_DEALLOCATE(mat_ls_ylm)
1069      LIBPAW_DEALLOCATE(mat_jmj)
1070      LIBPAW_DEALLOCATE(mat_out_c)
1071    end if ! tso
1072 
1073    LIBPAW_DEALLOCATE(ls_cplx)
1074    LIBPAW_DEALLOCATE(slm2ylm)
1075 
1076 !  End loop on l
1077  end do
1078 
1079  end subroutine lsylm

m_paw_sphharm/lxyz.F90 [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 lxyz

FUNCTION

 Computes the matrix element <Yl'm'|L_idir|Ylm>

INPUTS

   integer :: lp,mp,idir,ll,mm

OUTPUT

   complex(dpc) :: lidir

NOTES

  Ylm is the standard complex-valued spherical harmonic,
  idir is the direction in space of L

SOURCE

799 subroutine lxyz(lp,mp,idir,ll,mm,lidir)
800 
801 !Arguments ---------------------------------------------
802 !scalars
803  integer,intent(in) :: idir,ll,lp,mm,mp
804  complex(dpc),intent(out) :: lidir
805 
806 !Local variables ---------------------------------------
807 !scalars
808  complex(dpc) :: jme, jmme, jpme
809 
810 ! *********************************************************************
811 
812  lidir = czero
813  if ( lp /= ll ) return
814 
815  jpme=czero; jmme=czero; jme=czero
816  if (mp==mm) then 
817    jme=cone*mm
818  else if (mp==mm+1) then
819    jpme=-cone*sqrt(half*((ll*(ll+1))-mm*(mm+1)))
820  else if (mp==mm-1) then
821    jmme= cone*sqrt(half*((ll*(ll+1))-mm*(mm-1)))
822  end if
823 
824  select case (idir)
825    case (1) ! Lx
826      lidir = -sqrthalf*(jpme - jmme)
827    case (2) ! Ly
828      lidir = j_dpc*sqrthalf*(jpme + jmme)
829    case (3) ! Lz
830      lidir = jme
831  end select
832 
833 end subroutine lxyz

m_paw_sphharm/mat_mlms2jmj [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 mat_mlms2jmj

FUNCTION

 For a given angular momentum lcor, change a matrix of dimension 2(2*lcor+1)
 from the Ylm basis to the J,M_J basis if option==1

INPUTS

  lcor= angular momentum
  ndij= ndij = 4
  option=  1 matrix in |l,s,m_l,m_s> basis is changed into |l,s,j,m_j> basis
           2 matrix in |l,s,j,m_j> basis is changed into |l,s,m_l,m_s> basis
  optspin=  1  Spin up are first
            2  Spin dn are first
  prtvol=printing volume
  unitfi=printing file unit ; -1 for no printing
  wrt_mode=printing mode in parallel ('COLL' or 'PERS')

SIDE EFFECTS

  mat_mlms= Input/Output matrix in the Ylm basis, size of the matrix is (2*lcor+1,2*lcor+1,ndij)
  mat_jmj= Input/Output matrix in the J,M_J basis, size is 2*(2*lcor+1),2*(2*lcor+1)

NOTES

  usefull only in ndij==4

SOURCE

1668 subroutine mat_mlms2jmj(lcor,mat_mlms,mat_jmj,ndij,option,optspin,prtvol,unitfi,wrt_mode)
1669 
1670 !Arguments ---------------------------------------------
1671 !scalars
1672  integer,intent(in) :: ndij,lcor,option,optspin,prtvol,unitfi
1673  character(len=4),intent(in) :: wrt_mode
1674 !arrays
1675  complex(dpc),intent(inout) :: mat_mlms(2*lcor+1,2*lcor+1,ndij)
1676  complex(dpc),intent(inout) :: mat_jmj(2*(2*lcor+1),2*(2*lcor+1))
1677 
1678 !Local variables ---------------------------------------
1679 !scalars
1680  integer :: ii,im,im1,im2,ispden,jc1,jc2,jj,jm,ll,ml1,ml2,ms1,ms2
1681  real(dp),parameter :: invsqrt2=one/sqrt2
1682  real(dp) :: invsqrt2lp1,xj,xmj
1683  complex(dpc) :: mat_tmp,tmp2
1684  character(len=9),parameter :: dspinold(6)=(/"up       ","down     ","up-up    ","down-down","up-dn    ","dn-up    "/)
1685  character(len=9),parameter :: dspin(6)=(/"dn       ","up       ","dn-dn    ","up-up    ","dn-up    ","up-dn    "/)
1686  character(len=500) :: msg
1687 !arrays
1688  integer, allocatable :: ind_msml(:,:)
1689  complex(dpc),allocatable :: mat_mlms2(:,:),mlms2jmj(:,:)
1690 
1691 !*********************************************************************
1692 
1693  if(ndij/=4) then
1694    msg=" ndij/=4 !"
1695    LIBPAW_BUG(msg)
1696  end if
1697  if (option/=1.and.option/=2) then
1698    msg=' option=/1 and =/2 !'
1699    LIBPAW_BUG(msg)
1700  end if
1701  if (optspin/=1.and.optspin/=2) then
1702    msg=' optspin=/1 and =/2 !'
1703    LIBPAW_BUG(msg)
1704  end if
1705 
1706  if (unitfi/=-1) then
1707    if(option==1) then
1708      write(msg,'(3a)') ch10,&
1709 &     "matrix in |l,s,m_l,m_s> basis is changed into |l,s,j,m_j> basis"
1710      call wrtout(unitfi,msg,wrt_mode)
1711    else if(option==2) then
1712      write(msg,'(3a)') ch10,&
1713 &     "matrix in |l,s,j,m_j> basis is changed into |l,s,m_l,m_s> basis"
1714      call wrtout(unitfi,msg,wrt_mode)
1715    end if
1716  end if
1717 
1718  if(option==1) then
1719    if(optspin==2) then
1720      if(abs(prtvol)>2.and.unitfi/=-1)&
1721 &     write(msg,'(3a)') ch10,"assume spin dn is the first in the array"
1722    else if (optspin==1) then
1723      if(abs(prtvol)>2.and.unitfi/=-1)&
1724 &     write(msg,'(3a)') ch10,"change array in order that spin dn is the first in the array"
1725      do ii=1,2*lcor+1
1726        do jj=1,2*lcor+1
1727          mat_tmp=mat_mlms(ii,jj,2)
1728          mat_mlms(ii,jj,2)=mat_mlms(ii,jj,1)
1729          mat_mlms(ii,jj,1)=mat_tmp
1730          mat_tmp=mat_mlms(ii,jj,4)
1731          mat_mlms(ii,jj,4)=mat_mlms(ii,jj,3)
1732          mat_mlms(ii,jj,3)=mat_tmp
1733        end do
1734      end do
1735 !    mat_tmp(:,:,1)=mat_mlms(:,:,2);mat_tmp(:,:,2)=mat_mlms(:,:,1)
1736 !    mat_tmp(:,:,3)=mat_mlms(:,:,4);mat_tmp(:,:,4)=mat_mlms(:,:,3)
1737 !    mat_mlms(:,:,:)=mat_tmp(:,:,:)
1738    end if
1739    if(abs(prtvol)>2.and.unitfi/=-1) then
1740      call wrtout(unitfi,msg,wrt_mode)
1741    end if
1742  end if
1743 
1744  if(option==1.and.abs(prtvol)>2.and.unitfi/=-1) then
1745    do ispden=1,ndij
1746      write(msg,'(3a)') ch10,&
1747 &     "Input matrix in the Ylm basis for component ",trim(dspin(ispden+2*(ndij/4)))
1748      call wrtout(unitfi,msg,wrt_mode)
1749      do im1=1,lcor*2+1
1750        write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')&
1751 &       (mat_mlms(im1,im2,ispden),im2=1,lcor*2+1)
1752        call wrtout(unitfi,msg,wrt_mode)
1753      end do
1754    end do
1755  end if ! option==1
1756 
1757 !--------------- Built indices + allocations
1758  ll=lcor
1759  LIBPAW_ALLOCATE(mlms2jmj,(2*(2*ll+1),2*(2*ll+1)))
1760  mlms2jmj=czero
1761  LIBPAW_BOUND2_ALLOCATE(ind_msml,BOUNDS(1,2),BOUNDS(-ll,ll))
1762  LIBPAW_ALLOCATE(mat_mlms2,(2*(2*lcor+1),2*(2*lcor+1)))
1763  mlms2jmj=czero
1764  jc1=0
1765  do ms1=1,2
1766    do ml1=-ll,ll
1767      jc1=jc1+1
1768      ind_msml(ms1,ml1)=jc1
1769    end do
1770  end do
1771 !--------------- Change representation of input matrix for ndij==4
1772  if(option==1) then
1773    jc1=0
1774    do ms1=1,2
1775      do ml1=1,2*ll+1
1776        jc1=jc1+1
1777        jc2=0
1778        do ms2=1,2
1779          do ml2=1,2*ll+1
1780            jc2=jc2+1
1781            if(ms1==ms2) mat_mlms2(jc1,jc2)=mat_mlms(ml1,ml2,ms1)
1782            if(ms1<ms2) mat_mlms2(jc1,jc2)=mat_mlms(ml1,ml2,3)
1783            if(ms1>ms2) mat_mlms2(jc1,jc2)=mat_mlms(ml1,ml2,4)
1784          end do
1785        end do
1786      end do
1787    end do
1788    if(abs(prtvol)>1.and.unitfi/=-1) then
1789      write(msg,'(3a)') ch10,"Input matrix in the lms basis for all component"
1790      call wrtout(unitfi,msg,wrt_mode)
1791      do im1=1,2*(lcor*2+1)
1792        write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')&
1793 &       (mat_mlms2(im1,im2),im2=1,2*(lcor*2+1))
1794        call wrtout(unitfi,msg,wrt_mode)
1795      end do
1796    end if
1797  end if  ! option==1
1798 
1799 !--------------- built mlms2jmj
1800 !do jj=ll,ll+1    ! the physical value of j are ll-0.5,ll+0.5
1801 !xj(jj)=jj-0.5
1802  if(ll==0)then
1803    msg=' ll should not be equal to zero !'
1804    LIBPAW_BUG(msg)
1805  end if
1806  jc1=0
1807  invsqrt2lp1=one/sqrt(float(2*lcor+1))
1808  do jj=ll,ll+1
1809    xj=float(jj)-half
1810    do jm=-jj,jj-1
1811      xmj=float(jm)+half
1812      jc1=jc1+1
1813      if(nint(xj+0.5)==ll+1) then
1814        if(nint(xmj+0.5)==ll+1)  then
1815          mlms2jmj(ind_msml(2,ll),jc1)=1.0   !  J=L+0.5 and m_J=L+0.5
1816        else if(nint(xmj-0.5)==-ll-1) then
1817          mlms2jmj(ind_msml(1,-ll),jc1)=1.0   !  J=L+0.5 and m_J=-L-0.5
1818        else
1819          mlms2jmj(ind_msml(2,nint(xmj-0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5))
1820          mlms2jmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)-xmj+0.5))
1821        end if
1822      end if
1823      if(nint(xj+0.5)==ll) then
1824        mlms2jmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5))
1825        mlms2jmj(ind_msml(2,nint(xmj-0.5)),jc1)=-invsqrt2lp1*(sqrt(float(ll)-xmj+0.5))
1826      end if
1827    end do
1828  end do
1829  if(abs(prtvol)>2.and.unitfi/=-1) then
1830    write(msg,'(3a)') ch10,"Matrix to go from |M_L,M_S> to |J,M_J>"
1831    call wrtout(unitfi,msg,wrt_mode)
1832    do im1=1,2*(lcor*2+1)
1833      write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') (mlms2jmj(im1,im2),im2=1,2*(lcor*2+1))
1834      call wrtout(unitfi,msg,wrt_mode)
1835    end do
1836  end if
1837 
1838  do jm=1,2*(2*ll+1)
1839    do im=1,2*(2*ll+1)
1840      tmp2=czero
1841      do ii=1,2*(2*ll+1)
1842        do jj=1,2*(2*ll+1)
1843          if(option==1) then
1844            tmp2=tmp2+mat_mlms2(ii,jj)*CONJG(mlms2jmj(ii,im))*(mlms2jmj(jj,jm))
1845          else if(option==2) then
1846            tmp2=tmp2+mat_jmj(ii,jj)*(mlms2jmj(im,ii))*CONJG(mlms2jmj(jm,jj)) ! inv=t*
1847          end if
1848        end do
1849      end do
1850      if(option==1) then
1851        mat_jmj(im,jm)=tmp2
1852      else if(option==2) then
1853        mat_mlms2(im,jm)=tmp2
1854      end if
1855    end do
1856  end do
1857  if(option==1) then
1858    if (abs(prtvol)>=1.and.unitfi/=-1) then
1859      write(msg,'(3a)') ch10," Matrix in the J,M_J basis"
1860      call wrtout(unitfi,msg,wrt_mode)
1861      do im1=1,2*(lcor*2+1)
1862        write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') (mat_jmj(im1,im2),im2=1,2*(lcor*2+1))
1863        call wrtout(unitfi,msg,wrt_mode)
1864      end do
1865    end if
1866  else if(option==2) then
1867    if (abs(prtvol)>=1.and.unitfi/=-1) then
1868      write(msg,'(3a)') ch10," Matrix in the m_s m_l basis"
1869      call wrtout(unitfi,msg,wrt_mode)
1870      do im1=1,2*(lcor*2+1)
1871        write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') (mat_mlms2(im1,im2),im2=1,2*(lcor*2+1))
1872        call wrtout(unitfi,msg,wrt_mode)
1873      end do
1874    end if
1875    jc1=0
1876    do ms1=1,2
1877      do ml1=1,2*ll+1
1878        jc1=jc1+1
1879        jc2=0
1880        do ms2=1,2
1881          do ml2=1,2*ll+1
1882            jc2=jc2+1
1883            if(ms1==ms2) mat_mlms(ml1,ml2,ms1)=mat_mlms2(jc1,jc2)
1884            if(ms1<ms2) mat_mlms(ml1,ml2,3)=mat_mlms2(jc1,jc2)
1885            if(ms1>ms2) mat_mlms(ml1,ml2,4)=mat_mlms2(jc1,jc2)
1886          end do
1887        end do
1888      end do
1889    end do
1890  end if
1891  LIBPAW_DEALLOCATE(mlms2jmj)
1892  LIBPAW_DEALLOCATE(mat_mlms2)
1893  LIBPAW_DEALLOCATE(ind_msml)
1894 
1895  end subroutine mat_mlms2jmj

m_paw_sphharm/mat_slm2ylm [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 mat_slm2ylm

FUNCTION

 For a given angular momentum lcor, change a matrix  of dimension (2*lcor+1)
 from the Slm to the Ylm basis if option==1 or from Ylm to Slm if !option==2

INPUTS

  lcor= angular momentum, size of the matrix is 2(2*lcor+1)
  mat_inp_c= Input matrix
  ndij= ndij = 4
  option= -1  Change matrix from Slm to Ylm basis
           1  Change matrix from Ylm to Slm basis
  optspin=  1  Spin up are first
            2  Spin dn are first
  prtvol=printing volume
  unitfi=printing file unit ; -1 for no printing
  wrt_mode=printing mode in parallel ('COLL' or 'PERS')

OUTPUT

  mat_inp_c= Output matrix in Ylm or Slm basis according to option

NOTES

  usefull only in ndij==4

SOURCE

1928 subroutine mat_slm2ylm(lcor,mat_inp_c,mat_out_c,ndij,option,optspin,prtvol,unitfi,wrt_mode)
1929 
1930 !Arguments ---------------------------------------------
1931 !scalars
1932  integer,intent(in) :: ndij,lcor,option,optspin,prtvol,unitfi
1933  character(len=4),intent(in) :: wrt_mode
1934 !arrays
1935  complex(dpc) :: mat_inp_c(2*lcor+1,2*lcor+1,ndij),mat_out(2*lcor+1,2*lcor+1,ndij)
1936  complex(dpc) :: mat_out_c(2*lcor+1,2*lcor+1,ndij)
1937 
1938 !Local variables ---------------------------------------
1939 !scalars
1940  integer :: jm,ii,jj,ll,mm,ispden,im,im1,im2
1941  real(dp),parameter :: invsqrt2=one/sqrt2
1942  real(dp) :: onem
1943  complex(dpc) :: tmp2
1944  character(len=9),parameter :: dspinc(6)=(/"up       ","down     ","up-up    ","down-down","up-dn    ","dn-up    "/)! optspin 1
1945  character(len=9),parameter :: dspinc2(6)=(/"up       ","down     ","dn-dn    ","up-up    ","dn-up    ","up-dn    "/)! optspin 2
1946  character(len=500) :: msg
1947 !arrays
1948  complex(dpc),allocatable :: slm2ylm(:,:)
1949 
1950 ! *********************************************************************
1951 
1952  if(ndij/=4) then
1953    msg=' ndij:=4 !'
1954    LIBPAW_BUG(msg)
1955  end if
1956  if (option/=1.and.option/=2.and.option/=3.and.option/=4) then
1957    msg=' option=/1 or 2 or 3 or 4 !'
1958    LIBPAW_BUG(msg)
1959  end if
1960 
1961  if(abs(prtvol)>2.and.unitfi/=-1) then
1962    write(msg,'(3a)') ch10, "   mat_slm2ylm"
1963    call wrtout(unitfi,msg,wrt_mode)
1964  end if
1965 
1966  if(abs(prtvol)>2.and.unitfi/=-1) then
1967    if(option==1.or.option==3) then
1968      write(msg,'(3a)') ch10,"matrix in Slm basis is changed into Ylm basis"
1969      call wrtout(unitfi,msg,wrt_mode)
1970    else if(option==2.or.option==4) then
1971      write(msg,'(3a)') ch10,"matrix in Ylm basis is changed into Slm basis"
1972      call wrtout(unitfi,msg,wrt_mode)
1973    end if
1974  end if
1975 
1976  ll=lcor
1977  LIBPAW_ALLOCATE(slm2ylm,(2*ll+1,2*ll+1))
1978  slm2ylm=czero
1979  mat_out=zero
1980  mat_out_c=czero
1981  do im=1,2*ll+1
1982    mm=im-ll-1;jm=-mm+ll+1
1983    onem=dble((-1)**mm)
1984    if (mm> 0) then
1985      slm2ylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp)
1986      slm2ylm(jm,im)= cmplx(invsqrt2,     zero,kind=dp)
1987    end if
1988    if (mm==0) then
1989      slm2ylm(im,im)=cone
1990    end if
1991    if (mm< 0) then
1992      slm2ylm(im,im)= cmplx(zero,     invsqrt2,kind=dp)
1993      slm2ylm(jm,im)=-cmplx(zero,onem*invsqrt2,kind=dp)
1994    end if
1995  end do
1996  if(abs(prtvol)>2.and.unitfi/=-1) then
1997    do ispden=1,ndij
1998      if(optspin==1) then
1999        if(option==1.or.option==3)&
2000 &       write(msg,'(3a)') ch10,&
2001 &       "Input matrix in the Slm basis for component ",trim(dspinc(ispden+2*(ndij/4)))
2002        if(option==2.or.option==3)&
2003 &       write(msg,'(3a)') ch10,&
2004 &       "Input matrix in the Ylm basis for component ",trim(dspinc(ispden+2*(ndij/4)))
2005      else
2006        if(option==1.or.option==3)&
2007 &       write(msg,'(3a)') ch10,&
2008 &       "Input matrix in the Slm basis for component ",trim(dspinc2(ispden+2*(ndij/4)))
2009        if(option==2.or.option==3)&
2010 &       write(msg,'(3a)') ch10,&
2011 &       "Input matrix in the Ylm basis for component ",trim(dspinc2(ispden+2*(ndij/4)))
2012      end if
2013      call wrtout(unitfi,msg,wrt_mode)
2014      do im1=1,lcor*2+1
2015        write(msg,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
2016 &       (mat_inp_c(im1,im2,ispden),im2=1,lcor*2+1)
2017        call wrtout(unitfi,msg,wrt_mode)
2018      end do
2019    end do
2020  end if
2021  do ispden=1,ndij
2022    do jm=1,2*ll+1
2023      do im=1,2*ll+1
2024        tmp2=czero
2025        do ii=1,2*ll+1
2026          do jj=1,2*ll+1
2027            if(option==1) then
2028              tmp2=tmp2+mat_inp_c(ii,jj,ispden)*(slm2ylm(im,ii))*CONJG(slm2ylm(jm,jj))
2029            else if(option==2) then
2030              tmp2=tmp2+mat_inp_c(ii,jj,ispden)*CONJG(slm2ylm(ii,im))*(slm2ylm(jj,jm))
2031            end if
2032          end do
2033        end do
2034        mat_out_c(im,jm,ispden)=tmp2
2035      end do
2036    end do
2037  end do ! ispden
2038  do ii=1,2*ll+1
2039    do jj=1,2*ll+1
2040      mat_out(ii,jj,1)=real(mat_out_c(ii,jj,1))
2041      mat_out(ii,jj,2)=real(mat_out_c(ii,jj,2))
2042      mat_out(ii,jj,3)=real(mat_out_c(ii,jj,3))
2043      mat_out(ii,jj,4)=aimag(mat_out_c(ii,jj,3))
2044 !    check that n_{m,m'}^{alpha,beta}=conjg(n_{m',m"}^{beta,alpha}).
2045      if((abs(aimag(mat_out_c(ii,jj,3))+aimag(mat_out_c(jj,ii,4))).ge.0.0001).or. &
2046 &     (abs(real(mat_out_c(ii,jj,3))-real(mat_out_c(jj,ii,4))).ge.0.0001)) then
2047        write(msg,'(a,4f10.4)') &
2048 &       ' prb with mat_out_c ',mat_out_c(ii,jj,3),mat_out_c(ii,jj,4)
2049        LIBPAW_BUG(msg)
2050      end if
2051    end do
2052  end do
2053 
2054  LIBPAW_DEALLOCATE(slm2ylm)
2055 
2056 end subroutine mat_slm2ylm

m_paw_sphharm/mkeuler [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 mkeuler

FUNCTION

 Private function
 For a given symmetry operation, determines the corresponding Euler angles

INPUTS

  rot(3,3)= symmetry matrix

OUTPUT

  cosalp=  cos(alpha) with alpha=Euler angle 1
  cosbeta= cos(beta)  with beta =Euler angle 2
  cosgam=  cos(gamma) with gamma=Euler angle 3
  isn= error code (0 if the routine exit normally)
  sinalp= sin(alpha) with alpha=Euler angle 1
  sinbeta= sin(beta)  with beta=Euler angle 2
  singam= sin(gamma) with gamma=Euler angle 3

NOTES

  This file comes from the file crystal_symmetry.f
  by N.A.W. Holzwarth and A. Tackett for the code pwpaw
  XG20200718 However, this routine was not accurate in the determination
  of beta when cosbeta was close to one (indeed this is a special case). 
  This has been corrected. Moreover, sinbeta has been made an output in order
  to allow accurate calculations in dbeta. Also, tolerances have been made consistent.

SOURCE

3281 subroutine mkeuler(rot,cosbeta,sinbeta,cosalp,sinalp,cosgam,singam,isn)
3282 
3283 !Arguments ---------------------------------------------
3284 !scalars
3285  integer,intent(out) :: isn
3286  real(dp),intent(out) :: cosalp,cosbeta,cosgam,sinalp,sinbeta,singam
3287 !arrays
3288  real(dp),intent(in) :: rot(3,3)
3289 
3290 !Local variables ---------------------------------------
3291 !scalars
3292  integer :: ier
3293  real(dp) :: check,sinbeta2
3294  character(len=500) :: msg
3295 
3296 ! *********************************************************************
3297 
3298  do isn= -1,1,2
3299 
3300 !Old coding, inaccurate
3301 !  cosbeta=real(isn)*rot(3,3)
3302 !  if(abs(1._dp-cosbeta*cosbeta)<tol10) then
3303 !    sinbeta=zero
3304 !  else
3305 !    sinbeta=sqrt(1._dp-cosbeta*cosbeta)
3306 !  end if
3307 !  if (abs(sinbeta).gt.tol10)  then
3308 !    cosalp=isn*rot(3,1)/sinbeta
3309 !    sinalp=isn*rot(3,2)/sinbeta
3310 !    cosgam=-isn*rot(1,3)/sinbeta
3311 !    singam=isn*rot(2,3)/sinbeta
3312 !  else
3313 !    cosalp=isn*rot(1,1)/cosbeta
3314 !    sinalp=isn*rot(1,2)/cosbeta
3315 !    cosgam=one
3316 !    singam=zero
3317 !  end if
3318 
3319 !New coding, more accurate
3320    cosbeta=real(isn)*rot(3,3)
3321    sinbeta2=rot(1,3)**2+rot(2,3)**2
3322    if(sinbeta2<tol8**2)then
3323      sinbeta=zero
3324      cosalp=isn*rot(1,1)/cosbeta
3325      sinalp=isn*rot(1,2)/cosbeta
3326      cosgam=one
3327      singam=zero
3328    else
3329      sinbeta=sqrt(sinbeta2)
3330      cosalp=isn*rot(3,1)/sinbeta
3331      sinalp=isn*rot(3,2)/sinbeta
3332      cosgam=-isn*rot(1,3)/sinbeta
3333      singam=isn*rot(2,3)/sinbeta
3334    end if
3335 !
3336 
3337 !  Check matrix:
3338    ier=0
3339    check=cosalp*cosbeta*cosgam-sinalp*singam
3340    if (abs(check-isn*rot(1,1))>tol8) ier=ier+1
3341    check=sinalp*cosbeta*cosgam+cosalp*singam
3342    if (abs(check-isn*rot(1,2))>tol8) ier=ier+1
3343    check=-sinbeta*cosgam
3344    if (abs(check-isn*rot(1,3))>tol8) ier=ier+1
3345    check=-cosalp*cosbeta*singam-sinalp*cosgam
3346    if (abs(check-isn*rot(2,1))>tol8) ier=ier+1
3347    check=-sinalp*cosbeta*singam+cosalp*cosgam
3348    if (abs(check-isn*rot(2,2))>tol8) ier=ier+1
3349    check=sinbeta*singam
3350    if (abs(check-isn*rot(2,3))>tol8) ier=ier+1
3351    check=cosalp*sinbeta
3352    if (abs(check-isn*rot(3,1))>tol8) ier=ier+1
3353    check=sinalp*sinbeta
3354    if (abs(check-isn*rot(3,2))>tol8) ier=ier+1
3355    if (ier.eq.0) return
3356  end do
3357 
3358  isn=0
3359  write(msg, '(7a)' )&
3360 & 'Error during determination of symetries!',ch10,&
3361 & 'Action: check your input file:',ch10,&
3362 & 'unit cell vectors and/or atoms positions',ch10,&
3363 & 'have to be given with a better precision.'
3364  LIBPAW_ERROR(msg)
3365 
3366 end subroutine mkeuler

m_paw_sphharm/nablarealgaunt [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 nablarealgaunt

FUNCTION

 Evaluate integrals involving spherical harmonics and their gradient.
 These integrals are angular part for <nablaphi|nablaphj> and <tnablaphi|tnablaphj>
 Nabla_RealGaunt(ilm,ilm_i,ilm_j) = Int[ Slm Grad(Slm_i).Grad(Slm_j)]

INPUTS

  l_max = 1 + max. l value for Slm (see description above)
  l_max_ij = 1 + max. l value for Slm_i and Slm_j (see description above)

OUTPUT

 nnablagnt= number of non-zero integrals
 nabgauntselect(l_max**2,l_max_ij**2,l_max_ij**2)= indexes of the non-zero integrals
 nablagaunt(l_max**2*l_max_ij**4)= values of the integrals

SOURCE

2937 subroutine nablarealgaunt(l_max,l_max_ij,nnablagnt,nabgauntselect,nablagaunt) 
2938 
2939 !Arguments ---------------------------------------------
2940 !scalars
2941  integer, intent(in) :: l_max,l_max_ij
2942  integer, intent(out) :: nnablagnt
2943 !array
2944  integer,intent(out) :: nabgauntselect(:,:,:)
2945  real(dp),intent(out) :: nablagaunt(:)
2946 
2947 !Local variables ---------------------------------------
2948  logical,parameter :: debug=.false.
2949  integer :: angl_size,ii,ilm,ilm_i,ilm_j,ipt,mpsang,ntheta,nphi,ylm_size
2950  real(dp) :: nabla_rg, yylmgr
2951  character(len=500) :: msg
2952  real(dp),allocatable :: ang_wgth(:),cart_coord(:,:),ylmr(:,:),ylmrgr(:,:,:)
2953 
2954 !************************************************************************
2955 
2956  if ( size(nabgauntselect)< (l_max**2)*(l_max_ij**4) .or. &
2957 &     size(nablagaunt)    < (l_max**2)*(l_max_ij**4) ) then
2958    msg='Too small sizes for nabgauntselect/nablagaunt!'
2959    LIBPAW_BUG(msg)
2960  end if
2961 
2962  nabgauntselect(:,:,:)=-1
2963  nablagaunt(:)=zero
2964 
2965  ii=0
2966  if (l_max>1) then
2967    if (l_max_ij>=1) then
2968      ii=ii+1 ; nabgauntselect(1,2,2)=ii ; nablagaunt(ii)=0.5641895835477563_dp !(1/sqrt(pi)) 
2969      ii=ii+1 ; nabgauntselect(1,3,3)=ii ; nablagaunt(ii)=0.5641895835477563_dp !(1/sqrt(pi))
2970      ii=ii+1 ; nabgauntselect(1,4,4)=ii ; nablagaunt(ii)=0.5641895835477563_dp !(1/sqrt(pi))
2971    end if
2972    if (l_max_ij>2) then
2973      ii=ii+1 ; nabgauntselect(1,5,5)=ii ; nablagaunt(ii)=1.692568750643269_dp !\frac{3}{\sqrt{\pi}}
2974      ii=ii+1 ; nabgauntselect(1,6,6)=ii ; nablagaunt(ii)=1.692568750643269_dp !\dfrac{3}{\sqrt{\pi}}
2975      ii=ii+1 ; nabgauntselect(1,7,7)=ii ; nablagaunt(ii)=1.692568750643269_dp !\frac{3}{\sqrt{\pi}}
2976      ii=ii+1 ; nabgauntselect(1,8,8)=ii ; nablagaunt(ii)=1.692568750643269_dp !\frac{3}{\sqrt{\pi}}
2977      ii=ii+1 ; nabgauntselect(1,9,9)=ii ; nablagaunt(ii)=1.692568750643269_dp !\frac{3}{\sqrt{\pi}}
2978      ii=ii+1 ; nabgauntselect(2,2,7)=ii ; nablagaunt(ii)=-0.37846987830302403_dp !\frac{-3}{2\sqrt{5\pi}}
2979      ii=ii+1 ; nabgauntselect(2,2,9)=ii ; nablagaunt(ii)=-0.6555290583552474_dp !\frac{-1.5\sqrt{3}}{\sqrt{5\pi}}
2980      ii=ii+1 ; nabgauntselect(2,3,6)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
2981      ii=ii+1 ; nabgauntselect(2,4,5)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
2982      ii=ii+1 ; nabgauntselect(2,5,4)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
2983      ii=ii+1 ; nabgauntselect(2,6,3)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
2984      ii=ii+1 ; nabgauntselect(2,7,2)=ii ; nablagaunt(ii)=-0.37846987830302403_dp!\frac{-3}{2\sqrt{5\pi}}
2985      ii=ii+1 ; nabgauntselect(2,9,2)=ii ; nablagaunt(ii)=-0.6555290583552474_dp !\frac{3}{2}\sqrt{\frac{3}{5\pi}}
2986      ii=ii+1 ; nabgauntselect(3,2,6)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
2987      ii=ii+1 ; nabgauntselect(3,3,7)=ii ; nablagaunt(ii)=0.75693974607408354_dp !\frac{3}{\sqrt{5\pi}
2988      ii=ii+1 ; nabgauntselect(3,4,8)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
2989      ii=ii+1 ; nabgauntselect(3,6,2)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
2990      ii=ii+1 ; nabgauntselect(3,7,3)=ii ; nablagaunt(ii)=0.7569397566060481_dp !\frac{3}{\sqrt{5\pi}}
2991      ii=ii+1 ; nabgauntselect(3,8,4)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{3}{2}\sqrt{\frac{3}{5\pi}}
2992      ii=ii+1 ; nabgauntselect(4,2,5)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
2993      ii=ii+1 ; nabgauntselect(4,3,8)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
2994      ii=ii+1 ; nabgauntselect(4,4,7)=ii ; nablagaunt(ii)=-0.37846987830302403_dp !\frac{-3}{2\sqrt{5\pi}}
2995      ii=ii+1 ; nabgauntselect(4,4,9)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
2996      ii=ii+1 ; nabgauntselect(4,5,2)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{1.5\sqrt{3}}{\sqrt{5\pi}}
2997      ii=ii+1 ; nabgauntselect(4,7,4)=ii ; nablagaunt(ii)=-0.37846987830302403_dp!\frac{-3}{2\sqrt{5\pi}}
2998      ii=ii+1 ; nabgauntselect(4,8,3)=ii ; nablagaunt(ii)=0.6555290583552474_dp  !\frac{3}{2}\sqrt{\frac{3}{5\pi}}
2999      ii=ii+1 ; nabgauntselect(4,9,4)=ii ; nablagaunt(ii)=0.6555290583552474_dp !\frac{3}{2}\sqrt{\frac{3}{5\pi}}
3000    end if
3001  end if
3002 
3003  if (l_max>2) then
3004    if (l_max_ij>1) then
3005      ii=ii+1 ; nabgauntselect(5,2,4)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3006      ii=ii+1 ; nabgauntselect(5,4,2)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3007      ii=ii+1 ; nabgauntselect(6,2,3)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3008      ii=ii+1 ; nabgauntselect(6,3,2)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3009      ii=ii+1 ; nabgauntselect(7,2,2)=ii ; nablagaunt(ii)=0.126156626101008_dp !\frac{1}{2\sqrt{5\pi}}
3010      ii=ii+1 ; nabgauntselect(7,3,3)=ii ; nablagaunt(ii)=-0.252313252202016_dp !-\frac{1}{\sqrt{5\pi}}
3011      ii=ii+1 ; nabgauntselect(7,4,4)=ii ; nablagaunt(ii)=0.126156626101008_dp !\frac{1}{2\sqrt{5\pi}}
3012      ii=ii+1 ; nabgauntselect(8,3,4)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3013      ii=ii+1 ; nabgauntselect(8,4,3)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3014      ii=ii+1 ; nabgauntselect(9,2,2)=ii ; nablagaunt(ii)=0.2185096861184158_dp !\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3015      ii=ii+1 ; nabgauntselect(9,4,4)=ii ; nablagaunt(ii)=-0.2185096861184158_dp !-\frac{0.5\sqrt{3}}{\sqrt{5\pi}}
3016    end if
3017    if (l_max_ij>2) then
3018      ii=ii+1 ; nabgauntselect(5,5,7)=ii ; nablagaunt(ii)=-0.5406712547186058_dp !-\frac{3}{7}\sqrt{\frac{5}{\pi}}
3019      ii=ii+1 ; nabgauntselect(5,6,8)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3020      ii=ii+1 ; nabgauntselect(5,7,5)=ii ; nablagaunt(ii)=-0.5406712547186058_dp !-\frac{3}{7}\sqrt{\frac{5}{\pi}}
3021      ii=ii+1 ; nabgauntselect(5,8,6)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3022      ii=ii+1 ; nabgauntselect(6,5,8)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3023      ii=ii+1 ; nabgauntselect(6,6,7)=ii ; nablagaunt(ii)=0.2703356273593029_dp !\frac{3}{14}\sqrt{\frac{5}{\pi}}
3024      ii=ii+1 ; nabgauntselect(6,6,9)=ii ; nablagaunt(ii)=-0.4682350416823196_dp !-\frac{3}{14}\sqrt{\frac{15}{\pi}}
3025      ii=ii+1 ; nabgauntselect(6,7,6)=ii ; nablagaunt(ii)=0.2703356273593029_dp !\frac{3}{14}\sqrt{\frac{5}{\pi}}
3026      ii=ii+1 ; nabgauntselect(6,8,5)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3027      ii=ii+1 ; nabgauntselect(6,9,6)=ii ; nablagaunt(ii)=-0.4682350416823196_dp !-\frac{3}{14}\sqrt{\frac{15}{\pi}}
3028      ii=ii+1 ; nabgauntselect(7,5,5)=ii ; nablagaunt(ii)=-0.5406712547186058_dp !-\frac{3}{7}\sqrt{\frac{5}{\pi}}
3029      ii=ii+1 ; nabgauntselect(7,6,6)=ii ; nablagaunt(ii)=0.2703356273593029_dp !\frac{3}{14}\sqrt{\frac{5}{\pi}}
3030      ii=ii+1 ; nabgauntselect(7,7,7)=ii ; nablagaunt(ii)=0.5406712547186058_dp !\frac{3}{7}\sqrt{\frac{5}{\pi}}
3031      ii=ii+1 ; nabgauntselect(7,8,8)=ii ; nablagaunt(ii)=0.2703356273593029_dp !\frac{3}{14}\sqrt{\frac{5}{\pi}}
3032      ii=ii+1 ; nabgauntselect(7,9,9)=ii ; nablagaunt(ii)=-0.5406712547186058_dp !-\frac{3}{7}\sqrt{\frac{5}{\pi}}
3033      ii=ii+1 ; nabgauntselect(8,5,6)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3034      ii=ii+1 ; nabgauntselect(8,6,5)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3035      ii=ii+1 ; nabgauntselect(8,7,8)=ii ; nablagaunt(ii)=0.2703356273593029_dp !\frac{3}{14}\sqrt{\frac{5}{\pi}}
3036      ii=ii+1 ; nabgauntselect(8,8,7)=ii ; nablagaunt(ii)=0.2703356273593029_dp !\frac{3}{14}\sqrt{\frac{5}{\pi}}
3037      ii=ii+1 ; nabgauntselect(8,8,9)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3038      ii=ii+1 ; nabgauntselect(8,9,8)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3039      ii=ii+1 ; nabgauntselect(9,6,6)=ii ; nablagaunt(ii)=-0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3040      ii=ii+1 ; nabgauntselect(9,7,9)=ii ; nablagaunt(ii)=-0.5406712547186058_dp !-\frac{3}{7}\sqrt{\frac{5}{\pi}}
3041      ii=ii+1 ; nabgauntselect(9,8,8)=ii ; nablagaunt(ii)=0.4682350416823196_dp !\frac{3}{14}\sqrt{\frac{15}{\pi}}
3042      ii=ii+1 ; nabgauntselect(9,9,7)=ii ; nablagaunt(ii)=-0.5406712547186058_dp !-\frac{3}{7}\sqrt{\frac{5}{\pi}}
3043    end if
3044  end if
3045 
3046  nnablagnt=ii
3047 
3048 !If not tabulated, compute the integrals
3049  if (l_max>3.or.l_max_ij>3) then
3050  
3051    ntheta=25 ; nphi=25
3052    call ylm_angular_mesh(ntheta,nphi,angl_size,cart_coord,ang_wgth)
3053 
3054    mpsang=1+max(l_max,l_max_ij)
3055    ylm_size=mpsang**2
3056    LIBPAW_ALLOCATE(ylmr,(ylm_size,angl_size))
3057    LIBPAW_ALLOCATE(ylmrgr,(3,ylm_size,angl_size))
3058    call initylmr(mpsang,0,angl_size,ang_wgth,2,cart_coord,ylmr,ylmr_gr=ylmrgr)
3059  
3060    if (debug) open(unit=111,file='nablarealgaunt.dat',form='formatted')
3061 
3062    do ilm=1,l_max**2
3063      do ilm_i=1,l_max_ij**2
3064        do ilm_j=1,l_max_ij**2
3065 
3066          if (ilm<10.and.ilm_i<10.and.ilm_j<10) cycle  ! Already stored (tabulated)
3067 
3068          ! Compute integral
3069          nabla_rg=zero
3070          do ipt=1,angl_size
3071            yylmgr=ylmrgr(1,ilm_i,ipt)*ylmrgr(1,ilm_j,ipt) &
3072 &                +ylmrgr(2,ilm_i,ipt)*ylmrgr(2,ilm_j,ipt) &
3073 &                +ylmrgr(3,ilm_i,ipt)*ylmrgr(3,ilm_j,ipt)
3074            nabla_rg=nabla_rg+ang_wgth(ipt)*ylmr(ilm,ipt)*yylmgr
3075          end do
3076          nabla_rg=four_pi*nabla_rg
3077 
3078          ! Store it if non-zero
3079          if (abs(nabla_rg)>tol12) then
3080            if (debug) then
3081              write(111,'(5x,a,i2,a,i2,a,i2,a,f19.15,a)') &
3082 &              "ii=ii+1 ; nabgauntselect(",ilm,",",ilm_i,",",ilm_j, &
3083 &              ")=ii ; nablagaunt(ii)=",nabla_rg,"_dp"
3084            end if
3085            nnablagnt=nnablagnt+1
3086            nabgauntselect(ilm,ilm_i,ilm_j)=nnablagnt
3087            nablagaunt(nnablagnt)=nabla_rg
3088          end if
3089 
3090        end do ! ilm_j
3091      end do ! ilm_i
3092    end do !ilm
3093 
3094    if (debug) close(111)
3095    LIBPAW_DEALLOCATE(ylmr)
3096    LIBPAW_DEALLOCATE(ylmrgr)
3097    LIBPAW_DEALLOCATE(cart_coord)
3098    LIBPAW_DEALLOCATE(ang_wgth)
3099  end if
3100 
3101 end subroutine nablarealgaunt

m_paw_sphharm/perms [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 perms

FUNCTION

 Private function
 Returns N!/(N-k)!  if N>=0 and N>k ; otherwise 0 is returned

INPUTS

   kk=number k to use
   nn=number N to use

OUTPUT

   perms= n!/(n-k)!

SOURCE

3605 function perms(nn,kk)
3606 
3607 !Arguments ---------------------------------------------
3608 !scalars
3609  integer,intent(in) :: kk,nn
3610  real(dp) :: perms
3611 
3612 !Local variables ---------------------------------------
3613 !scalars
3614  integer :: ii
3615  real(dp) :: pp
3616 
3617 ! *********************************************************************
3618 
3619  if (nn>=0.and.nn>=kk) then
3620    pp=1._dp
3621    do ii=nn-kk+1,nn
3622      pp=pp*ii
3623    end do
3624  else
3625    pp=0._dp
3626  end if
3627 
3628  perms=pp
3629 
3630 end function perms

m_paw_sphharm/phim [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 phim

FUNCTION

 Computes Phi_m[theta]=Sqrt[2] cos[m theta],      if m>0
                       Sqrt[2] sin[Abs(m) theta], if m<0
                       1                        , if m=0

INPUTS

  costeta= cos(theta)  (theta= input angle)
  mm = index m
  sinteta= sin(theta)  (theta= input angle)

OUTPUT

  phim= Phi_m(theta) (see above)

NOTES

  - This file comes from the file crystal_symmetry.f
    by N.A.W. Holzwarth and A. Tackett for the code pwpaw

SOURCE

3394 pure function phim(costheta,sintheta,mm)
3395 
3396 !Arguments ---------------------------------------------
3397 !scalars
3398  integer,intent(in) :: mm
3399  real(dp) :: phim
3400  real(dp),intent(in) :: costheta,sintheta
3401 
3402 ! *********************************************************************
3403 
3404  if (mm==0)  phim=one
3405  if (mm==1)  phim=sqrt2*costheta
3406  if (mm==-1) phim=sqrt2*sintheta
3407  if (mm==2)  phim=sqrt2*(costheta*costheta-sintheta*sintheta)
3408  if (mm==-2) phim=sqrt2*two*sintheta*costheta
3409  if (mm==3)  phim=sqrt2*&
3410 & (costheta*(costheta*costheta-sintheta*sintheta)&
3411 & -sintheta*two*sintheta*costheta)
3412  if (mm==-3) phim=sqrt2*&
3413 & (sintheta*(costheta*costheta-sintheta*sintheta)&
3414 & +costheta*two*sintheta*costheta)
3415 
3416  end function phim

m_paw_sphharm/pl_deriv [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 pl_deriv

FUNCTION

 Compute d2(Pl (x)))/d(x)2  where P_l is a legendre polynomial

INPUTS

  mpsang=1+ maximum l quantum number
  xx= input value

OUTPUT

  pl_d2(mpsang*mpsang)

SOURCE

1516 subroutine pl_deriv(mpsang,pl_d2,xx)
1517 
1518 !Arguments ---------------------------------------------
1519 !scalars
1520  integer,intent(in) :: mpsang
1521  real(dp),intent(in) :: xx
1522 !arrays
1523  real(dp),intent(out) :: pl_d2(mpsang)
1524 
1525 !Local variables ---------------------------------------
1526 !scalars
1527  integer :: il,ilm
1528  real(dp) :: il_,il_m1,il_2m1
1529  character(len=500) :: msg
1530 !arrays
1531  real(dp) :: pl(mpsang),pl_d1(mpsang)
1532 
1533 ! *********************************************************************
1534 
1535  if (abs(xx).gt.1.d0) then
1536    msg = 'pl_deriv : xx > 1 !'
1537    LIBPAW_ERROR(msg)
1538  end if
1539 
1540  pl_d2=zero; pl_d1=zero; pl=zero
1541  pl(1)=one; pl(2)=xx
1542  pl_d1(1)=zero; pl_d1(2)=one
1543  pl_d2(1)=zero; pl_d2(2)=zero
1544  if (mpsang>2) then
1545    do il=2,mpsang-1
1546      il_=dble(il);il_m1=dble(il-1);il_2m1=dble(2*il-1)
1547      ilm=il+1
1548      pl(ilm)=(il_2m1*xx*pl(ilm-1)-il_m1*pl(ilm-2))/il_
1549      pl_d1(ilm)=(il_2m1*(xx*pl_d1(ilm-1)+pl(ilm-1))-il_m1*pl_d1(ilm-2))/il_
1550      pl_d2(ilm)=(il_2m1*(xx*pl_d2(ilm-1)+two*pl_d1(ilm-1))-il_m1*pl_d2(ilm-2))/il_
1551    end do
1552  end if
1553 
1554 end subroutine pl_deriv

m_paw_sphharm/plm_coeff [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 plm_coeff

FUNCTION

 Compute coefficients depending on Plm and its derivatives where P_lm is a legendre polynomial.
 They are used to compute the second derivatives of spherical harmonics

INPUTS

  mpsang=1+ maximum l quantum number
  xx= input value

OUTPUT

  blm(5,mpsang*mpsang)=coefficients depending on Plm and its derivatives where P_lm is a legendre polynome

SOURCE

1101 subroutine plm_coeff(blm,mpsang,xx)
1102 
1103 !Arguments ---------------------------------------------
1104 !scalars
1105  integer,intent(in) :: mpsang
1106  real(dp),intent(in) :: xx
1107 !arrays
1108  real(dp),intent(out) :: blm(5,mpsang*mpsang)
1109 
1110 !Local variables ---------------------------------------
1111 !scalars
1112  integer :: il,ilm,ilm0,ilm1,im
1113  real(dp) :: dplm_dt,d2plm_dt2,llp1,onemx2,plm,sqrx,xsqrx,xx2,yy
1114  logical :: is_one
1115  character(len=500) :: msg
1116 !arrays
1117  real(dp) :: pl_d2(mpsang),plm_d2t(mpsang*mpsang)
1118 
1119 !************************************************************************
1120 
1121  if (abs(xx).gt.1.d0) then
1122    msg = ' plm_coeff :  xx > 1 !'
1123    LIBPAW_ERROR(msg)
1124  end if
1125 
1126  blm=zero
1127  is_one=(abs(abs(xx)-one)<=tol12)
1128  xx2=xx**2
1129  onemx2=abs(one-xx2)
1130  sqrx=sqrt(onemx2)
1131  xsqrx=xx*sqrt(onemx2)
1132 
1133  call plm_d2theta(mpsang,plm_d2t,xx)
1134  if (is_one) then
1135    yy=sign(one,xx)
1136    call pl_deriv(mpsang,pl_d2,yy)
1137  end if
1138 
1139  do il=0,mpsang-1
1140    llp1=dble(il*(il+1))
1141    ilm0=il*il+il+1
1142    do im=0,il
1143      ilm=ilm0+im;ilm1=ilm0-im
1144 
1145      plm      =(-1)**im*ass_leg_pol(il,im,xx)
1146      dplm_dt  =(-1)**im*plm_dtheta(il,im,xx)
1147      d2plm_dt2=         plm_d2t(ilm)
1148 
1149      blm(1,ilm)=         two*xsqrx    *dplm_dt+onemx2*d2plm_dt2
1150      blm(2,ilm)=         (one-two*xx2)*dplm_dt-xsqrx *d2plm_dt2
1151      blm(3,ilm)=llp1*plm+                             d2plm_dt2
1152      blm(4,ilm)=        -two*xsqrx    *dplm_dt+xx2   *d2plm_dt2
1153 
1154 
1155      if (is_one) then
1156        if (im==1) then
1157          blm(5,ilm)=llp1*plm+d2plm_dt2
1158        end if
1159        if (im==2) then
1160          blm(5,ilm)=d2plm_dt2-three*pl_d2(il+1)
1161        end if
1162      else
1163        if(im>0) then
1164          blm(5,ilm)=plm/onemx2-dplm_dt*xx/sqrx
1165        end if
1166      end if
1167 
1168      if (im>0) then
1169        blm(1,ilm1)=blm(1,ilm)
1170        blm(2,ilm1)=blm(2,ilm)
1171        blm(3,ilm1)=blm(3,ilm)
1172        blm(4,ilm1)=blm(4,ilm)
1173        blm(5,ilm1)=blm(5,ilm)
1174      end if
1175 
1176    end do
1177  end do
1178 
1179 end subroutine plm_coeff

m_paw_sphharm/plm_d2theta [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 plm_d2theta

FUNCTION

 Compute d2(Plm (cos(theta)))/d(theta)2  where P_lm is a legendre polynome

INPUTS

  mpsang=1+ maximum l quantum number
  xx= input value

OUTPUT

  plm_d2t(mpsang*mpsang)

SOURCE

1442 subroutine plm_d2theta(mpsang,plm_d2t,xx)
1443 
1444 !Arguments ---------------------------------------------
1445 !scalars
1446  integer,intent(in) :: mpsang
1447  real(dp),intent(in) :: xx
1448 !arrays
1449  real(dp),intent(out) :: plm_d2t(mpsang*mpsang)
1450 
1451 !Local variables ---------------------------------------
1452 !scalars
1453  integer :: il,ilm,ilmm1,ilmm2,im
1454  real(dp) :: sqrx
1455  character(len=500) :: msg
1456 
1457 !************************************************************************
1458  if (abs(xx).gt.1.d0) then
1459    msg = 'plm_d2theta : xx > 1 !'
1460    LIBPAW_ERROR(msg)
1461  end if
1462 
1463  plm_d2t=zero
1464  if (mpsang>1) then
1465    sqrx=sqrt(abs((1.d0-xx)*(1.d0+xx)))
1466 
1467    do il=1,mpsang-1
1468      ilm=il*il+2*il+1
1469      ilmm1=(il-1)*(il-1)+2*(il-1)+1
1470 !    terme d2(Pll)/dtet2
1471      plm_d2t(ilm)=(2*il-1)*(sqrx*(plm_d2t(ilmm1)-(-1)**(il-1)*ass_leg_pol(il-1,il-1,xx))+&
1472 &     2.d0*xx*(-1)**(il-1)*plm_dtheta(il-1,il-1,xx))
1473      plm_d2t(ilm-2*il)=plm_d2t(ilm)
1474 !    terme d2(Pl(l-1))/dtet2
1475      plm_d2t(ilm-1)=(2*il-1)*(xx*(plm_d2t(ilmm1)-(-1)**(il-1)*ass_leg_pol(il-1,il-1,xx))-&
1476 &     2.d0*sqrx*(-1)**(il-1)*plm_dtheta(il-1,il-1,xx))
1477      if(il>1) plm_d2t(il*il+2)=plm_d2t(ilm-1)
1478    end do
1479 !  terme d2(Plm)/dtet2
1480    if(mpsang>2) then
1481      do il=2,mpsang-1
1482        do im=0,il-2
1483          ilm=il*il+il+1+im
1484          ilmm1=(il-1)*(il-1)+il+im
1485          ilmm2=(il-2)*(il-2)+il-1+im
1486          plm_d2t(ilm)=dble(2*il-1)/dble(il-im)*(xx*(plm_d2t(ilmm1)-(-1)**im*ass_leg_pol(il-1,im,xx))-&
1487 &         2.d0*sqrx*(-1)**im*plm_dtheta(il-1,im,xx))-&
1488 &         dble(il+im-1)/dble(il-im)*plm_d2t(ilmm2)
1489          plm_d2t(ilm-2*im)=plm_d2t(ilm)
1490        end do
1491      end do
1492    end if
1493  end if
1494 
1495 end subroutine plm_d2theta

m_paw_sphharm/plm_dphi [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 plm_dphi

FUNCTION

 Compute  m*P_lm(x)/sqrt((1-x^2)where P_lm is a legendre polynome

INPUTS

  ll= l quantum number
  mm= m quantum number
  xx= input value

OUTPUT

  plm_dphi(xx)

NOTES

  This routine comes from Function Der_Phi_P(L,m,x)
  (pwpaw code from N. Holzwarth, implemented by Y. Abraham))

SOURCE

1277 function plm_dphi(ll,mm,xx)
1278 
1279 !Arguments ---------------------------------------------
1280 !scalars
1281  integer,intent(in) :: ll,mm
1282  real(dp) :: plm_dphi
1283  real(dp),intent(in) :: xx
1284 
1285 !Local variables ---------------------------------------
1286 !scalars
1287  integer :: il,im
1288  real(dp) :: dosomx2,fact,pll,pmm,pmmp1,somx2
1289  character(len=500) :: msg
1290 
1291 ! *********************************************************************
1292 
1293  if (mm.lt.0.or.mm.gt.ll.or.abs(xx).gt.1.d0) then
1294    msg = 'plm_dphi : mm < 0 or mm > ll or xx > 1 !'
1295    LIBPAW_ERROR(msg)
1296  end if
1297 
1298  plm_dphi=zero
1299  if (mm==0) return
1300 
1301  pmm=one
1302  dosomx2=one
1303  if (mm > 0) then
1304    somx2=sqrt((1-xx)*(1+xx))
1305    fact=one
1306    do im=1,mm
1307      pmm=-pmm*fact
1308      fact=fact+2
1309    end do
1310    if (mm > 1) then
1311      do im=2,mm
1312        dosomx2=somx2*dosomx2
1313      end do
1314    end if
1315    pmm=pmm*dosomx2 !due to one more term (-1^M)
1316  end if
1317  if(ll==mm) then
1318    plm_dphi=pmm*mm
1319  else
1320    pmmp1=xx*(2*mm+1)*pmm
1321    if(ll==mm+1) then
1322      plm_dphi=pmmp1*mm
1323    else if(ll>=mm+2) then
1324      do il=mm+2,ll
1325        pll=(xx*(2*il-1)*pmmp1-(il+mm-1)*pmm)/(il-mm)
1326        pmm=pmmp1
1327        pmmp1=pll
1328      end do
1329      plm_dphi=pll*mm
1330    end if
1331  end if
1332 
1333 end function plm_dphi

m_paw_sphharm/plm_dtheta [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 plm_dtheta

FUNCTION

 Compute -(1-x^2)^1/2*d/dx{P_lm(x)} where P_lm is a legendre polynome

INPUTS

  ll= l quantum number
  mm= m quantum number
  xx= input value

OUTPUT

  plm_dtheta(xx)

NOTES

  This routine comes from Function Der_Theta_P(L,m,x)
  (pwpaw code from N. Holzwarth, implemented by Y. Abraham))

SOURCE

1359 function plm_dtheta(ll,mm,xx)
1360 
1361 !Arguments ---------------------------------------------
1362 !scalars
1363  integer,intent(in) :: ll,mm
1364  real(dp) :: plm_dtheta
1365  real(dp),intent(in) :: xx
1366 
1367 !Local variables ---------------------------------------
1368 !scalars
1369  integer :: il,im
1370  real(dp) :: dosomx2,dpll,dpmm,dpmmp1,fact,pll,pmm,pmmp1,somx2
1371  character(len=500) :: msg
1372 
1373 ! *********************************************************************
1374 
1375  if (mm.lt.0.or.mm.gt.ll.or.abs(xx).gt.1.d0) then
1376    msg = 'plm_dtheta : mm < 0 or mm > ll or xx > 1 !'
1377    LIBPAW_ERROR(msg)
1378  end if
1379 
1380  plm_dtheta=zero
1381  pmm=one
1382  dpmm=one
1383  dosomx2=one
1384  somx2=sqrt((1-xx)*(1+xx))
1385  if(mm==0)then
1386    dpmm=zero
1387  elseif (mm > 0) then
1388    fact=one
1389    do im=1,mm
1390      pmm=-pmm*fact*somx2
1391      dpmm=-dpmm*fact
1392      fact=fact+2
1393    end do
1394    if(mm>1)then
1395      do im=2,mm
1396        dosomx2=dosomx2*somx2
1397      end do
1398    end if
1399    dpmm= dpmm*mm*xx*dosomx2
1400  end if
1401  if(ll==mm)then
1402    plm_dtheta=dpmm
1403  else
1404    pmmp1=xx*(2*mm+1)*pmm
1405    dpmmp1=-(2*mm+1)*somx2*pmm+xx*(2*mm+1)*dpmm
1406    if(ll==mm+1) then
1407      plm_dtheta=dpmmp1
1408    else if(ll>=mm+2)then
1409      do il=mm+2,ll
1410        pll=(xx*(2*il-1)*pmmp1-(il+mm-1)*pmm)/(il-mm)
1411        dpll=(-somx2*(2*il-1)*pmmp1+(xx*(2*il-1)*dpmmp1-(il+mm-1)*dpmm))/(il-mm)
1412        pmm=pmmp1
1413        pmmp1=pll
1414        dpmm=dpmmp1
1415        dpmmp1=dpll
1416      end do
1417      plm_dtheta=dpll
1418    end if
1419  end if
1420 
1421 end function plm_dtheta

m_paw_sphharm/rfactorial [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 rfactorial

FUNCTION

 Private function
 Calculates N! as a double precision real.

INPUTS

   nn=number to use

OUTPUT

   factorial= n! (real)

SOURCE

3651 elemental function rfactorial(nn)
3652 
3653 !Arguments ---------------------------------------------
3654 !scalars
3655  integer,intent(in) :: nn
3656  real(dp) :: rfactorial
3657 
3658 !Local variables ---------------------------------------
3659 !scalars
3660  integer :: ii
3661 
3662 ! *********************************************************************
3663 
3664  rfactorial=one
3665  do ii=2,nn
3666    rfactorial=rfactorial*ii
3667  end do
3668 
3669 end function rfactorial

m_paw_sphharm/setnabla_ylm [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 setnabla_ylm

FUNCTION

 Evaluate several inegrals involving spherical harmonics and their gradient.
 These integrals are angular part for <phi_i|nabla|phi_j> and <tphi_i|nabla|tphi_j>.

INPUTS

  mpsang=1+ max. angular momentum

OUTPUT

  ang_phipphj :: angular part for <phi_i|nabla|phi_j> and <tphi_i|nabla|tphi_j>
  ang_phipphj(i,j,1)=\int sin\theta cos\phi Si Sj d\omega
  ang_phipphj(i,j,2)=\int cos\theta cos\phi Si \frac{d}{d\theta}Sj d\Omega
  ang_phipphj(i,j,3)=\int -sin\phi  Si \frac{d}{d\phi}Sj d\Omega
  ang_phipphj(i,j,4)=\int sin\theta sin\phi Si Sj d\Omega
  ang_phipphj(i,j,5)=\int cos\theta sin\phi Si \frac{d}{d\theta}Sj d\Omega
  ang_phipphj(i,j,6)=\int cos\phi Si \frac{d}{d\phi}Sj d\Omega
  ang_phipphj(i,j,7)=\int cos\theta  Si Sj d\Omega
  ang_phipphj(i,j,8)=\int -sin\theta Si \frac{d}{d\theta}Sj d\Omega

NOTES

   See : Mazevet, S., Torrent, M., Recoules, V. and Jollet, F., High Energy Density Physics, 6, 84-88 (2010)
         Calculations of the Transport Properties within the PAW Formalism

SOURCE

2250  subroutine setnabla_ylm(ang_phipphj,mpsang)
2251 
2252 !Arguments ------------------------------------
2253 !scalars
2254  integer,intent(in) :: mpsang
2255 !arrays
2256  real(dp),intent(out) :: ang_phipphj(mpsang**2,mpsang**2,8)
2257 
2258 !Local variables-------------------------------
2259  character(len=500) :: msg
2260  real(dp) :: ang_phipphj_tmp(16,16,8)
2261 
2262 ! ************************************************************************
2263 
2264  if (mpsang>4) then
2265    msg='  Not designed for angular momentum greater than 3!'
2266    LIBPAW_ERROR(msg)
2267  end if
2268 
2269 !8 angular integrals for l=0..3, m=-l..+l
2270 !ang_phipphj(1,4,1)=\frac{1}{\sqrt{3}}
2271 !ang_phipphj(2,5,1)=\frac{1}{\sqrt{5}}
2272 !ang_phipphj(3,8,1)=\frac{1}{\sqrt{5}}
2273 !ang_phipphj(4,1,1)=\frac{1}{\sqrt{3}}
2274 !ang_phipphj(4,7,1)=-\frac{1}{\sqrt{15}}
2275 !ang_phipphj(4,9,1)=\frac{1}{\sqrt{5}}
2276 !ang_phipphj(5,2,1)=\frac{1}{\sqrt{5}}
2277 !ang_phipphj(5,10,1)=\sqrt{\frac{3}{14}}
2278 !ang_phipphj(5,12,1)=-\frac{1}{\sqrt{70}}
2279 !ang_phipphj(6,11,1)=\frac{1}{\sqrt{7}}
2280 !ang_phipphj(7,4,1)=-\frac{1}{\sqrt{15}}
2281 !ang_phipphj(7,14,1)=\sqrt{\frac{6}{35}}
2282 !ang_phipphj(8,3,1)=\frac{1}{\sqrt{5}}
2283 !ang_phipphj(8,13,1)=-\sqrt{\frac{3}{35}}
2284 !ang_phipphj(8,15,1)=\frac{1}{\sqrt{7}}
2285 !ang_phipphj(9,4,1)=\frac{1}{\sqrt{5}}
2286 !ang_phipphj(9,14,1)=-\frac{1}{\sqrt{70}}
2287 !ang_phipphj(9,16,1)=\sqrt{\frac{3}{14}}
2288 !ang_phipphj(10,5,1)=\sqrt{\frac{3}{14}}
2289 !ang_phipphj(11,6,1)=\frac{1}{\sqrt{7}}
2290 !ang_phipphj(12,5,1)=-\frac{1}{\sqrt{70}}
2291 !ang_phipphj(13,8,1)=-\sqrt{\frac{3}{35}}
2292 !ang_phipphj(14,7,1)=\sqrt{\frac{6}{35}}
2293 !ang_phipphj(14,9,1)=-\frac{1}{\sqrt{70}}
2294 !ang_phipphj(15,8,1)=\frac{1}{\sqrt{7}}
2295 !ang_phipphj(16,9,1)=\sqrt{\frac{3}{14}}
2296 !ang_phipphj(1,4,2)=\frac{1}{2 \sqrt{3}}
2297 !ang_phipphj(1,14,2)=-\frac{\sqrt{\frac{7}{6}}}{2}
2298 !ang_phipphj(2,5,2)=\frac{1}{2 \sqrt{5}}
2299 !ang_phipphj(3,8,2)=\frac{1}{2 \sqrt{5}}
2300 !ang_phipphj(4,7,2)=-\sqrt{\frac{3}{5}}
2301 !ang_phipphj(4,9,2)=\frac{1}{2 \sqrt{5}}
2302 !ang_phipphj(5,2,2)=\frac{1}{4 \sqrt{5}}
2303 !ang_phipphj(5,10,2)=\frac{\sqrt{\frac{3}{14}}}{2}
2304 !ang_phipphj(5,12,2)=-2 \sqrt{\frac{2}{35}}
2305 !ang_phipphj(6,11,2)=\frac{1}{2 \sqrt{7}}
2306 !ang_phipphj(7,4,2)=\frac{1}{\sqrt{15}}
2307 !ang_phipphj(7,14,2)=\frac{13}{2 \sqrt{210}}
2308 !ang_phipphj(8,3,2)=-\frac{1}{\sqrt{5}}
2309 !ang_phipphj(8,13,2)=-4 \sqrt{\frac{3}{35}}
2310 !ang_phipphj(8,15,2)=\frac{1}{2 \sqrt{7}}
2311 !ang_phipphj(9,4,2)=\frac{1}{4 \sqrt{5}}
2312 !ang_phipphj(9,14,2)=-2 \sqrt{\frac{2}{35}}
2313 !ang_phipphj(9,16,2)=\frac{\sqrt{\frac{3}{14}}}{2}
2314 !ang_phipphj(10,5,2)=\frac{1}{\sqrt{42}}
2315 !ang_phipphj(11,6,2)=-\frac{1}{4 \sqrt{7}}
2316 !ang_phipphj(12,5,2)=\sqrt{\frac{2}{35}}
2317 !ang_phipphj(13,8,2)=2 \sqrt{\frac{3}{35}}
2318 !ang_phipphj(14,7,2)=-2 \sqrt{\frac{6}{35}}
2319 !ang_phipphj(14,9,2)=\sqrt{\frac{2}{35}}
2320 !ang_phipphj(15,8,2)=-\frac{1}{4 \sqrt{7}}
2321 !ang_phipphj(16,9,2)=\frac{1}{\sqrt{42}}
2322 !ang_phipphj(1,4,3)=\frac{\sqrt{3}}{2}
2323 !ang_phipphj(1,14,3)=\frac{\sqrt{\frac{7}{6}}}{2}
2324 !ang_phipphj(2,5,3)=\frac{\sqrt{5}}{2}
2325 !ang_phipphj(3,8,3)=\frac{\sqrt{5}}{2}
2326 !ang_phipphj(4,9,3)=\frac{\sqrt{5}}{2}
2327 !ang_phipphj(5,2,3)=-\frac{\sqrt{5}}{4}
2328 !ang_phipphj(5,10,3)=\frac{\sqrt{\frac{21}{2}}}{2}
2329 !ang_phipphj(6,11,3)=\frac{\sqrt{7}}{2}
2330 !ang_phipphj(7,14,3)=\frac{\sqrt{\frac{35}{6}}}{2}
2331 !ang_phipphj(8,15,3)=\frac{\sqrt{7}}{2}
2332 !ang_phipphj(9,4,3)=-\frac{\sqrt{5}}{4}
2333 !ang_phipphj(9,16,3)=\frac{\sqrt{\frac{21}{2}}}{2}
2334 !ang_phipphj(10,5,3)=-\sqrt{\frac{7}{6}}
2335 !ang_phipphj(11,6,3)=-\frac{\sqrt{7}}{4}
2336 !ang_phipphj(15,8,3)=-\frac{\sqrt{7}}{4}
2337 !ang_phipphj(16,9,3)=-\sqrt{\frac{7}{6}}
2338 !ang_phipphj(1,2,4)=\frac{1}{\sqrt{3}}
2339 !ang_phipphj(2,1,4)=\frac{1}{\sqrt{3}}
2340 !ang_phipphj(2,7,4)=-\frac{1}{\sqrt{15}}
2341 !ang_phipphj(2,9,4)=-\frac{1}{\sqrt{5}}
2342 !ang_phipphj(3,6,4)=\frac{1}{\sqrt{5}}
2343 !ang_phipphj(4,5,4)=\frac{1}{\sqrt{5}}
2344 !ang_phipphj(5,4,4)=\frac{1}{\sqrt{5}}
2345 !ang_phipphj(5,14,4)=-\frac{1}{\sqrt{70}}
2346 !ang_phipphj(5,16,4)=-\sqrt{\frac{3}{14}}
2347 !ang_phipphj(6,3,4)=\frac{1}{\sqrt{5}}
2348 !ang_phipphj(6,13,4)=-\sqrt{\frac{3}{35}}
2349 !ang_phipphj(6,15,4)=-\frac{1}{\sqrt{7}}
2350 !ang_phipphj(7,2,4)=-\frac{1}{\sqrt{15}}
2351 !ang_phipphj(7,12,4)=\sqrt{\frac{6}{35}}
2352 !ang_phipphj(8,11,4)=\frac{1}{\sqrt{7}}
2353 !ang_phipphj(9,2,4)=-\frac{1}{\sqrt{5}}
2354 !ang_phipphj(9,10,4)=\sqrt{\frac{3}{14}}
2355 !ang_phipphj(9,12,4)=\frac{1}{\sqrt{70}}
2356 !ang_phipphj(10,9,4)=\sqrt{\frac{3}{14}}
2357 !ang_phipphj(11,8,4)=\frac{1}{\sqrt{7}}
2358 !ang_phipphj(12,7,4)=\sqrt{\frac{6}{35}}
2359 !ang_phipphj(12,9,4)=\frac{1}{\sqrt{70}}
2360 !ang_phipphj(13,6,4)=-\sqrt{\frac{3}{35}}
2361 !ang_phipphj(14,5,4)=-\frac{1}{\sqrt{70}}
2362 !ang_phipphj(15,6,4)=-\frac{1}{\sqrt{7}}
2363 !ang_phipphj(16,5,4)=-\sqrt{\frac{3}{14}}
2364 !ang_phipphj(1,2,5)=\frac{1}{2 \sqrt{3}}
2365 !ang_phipphj(1,12,5)=-\frac{\sqrt{\frac{7}{6}}}{2}
2366 !ang_phipphj(2,7,5)=-\sqrt{\frac{3}{5}}
2367 !ang_phipphj(2,9,5)=-\frac{1}{2 \sqrt{5}}
2368 !ang_phipphj(3,6,5)=\frac{1}{2 \sqrt{5}}
2369 !ang_phipphj(4,5,5)=\frac{1}{2 \sqrt{5}}
2370 !ang_phipphj(5,4,5)=\frac{1}{4 \sqrt{5}}
2371 !ang_phipphj(5,14,5)=-2 \sqrt{\frac{2}{35}}
2372 !ang_phipphj(5,16,5)=-\frac{\sqrt{\frac{3}{14}}}{2}
2373 !ang_phipphj(6,3,5)=-\frac{1}{\sqrt{5}}
2374 !ang_phipphj(6,13,5)=-4 \sqrt{\frac{3}{35}}
2375 !ang_phipphj(6,15,5)=-\frac{1}{2 \sqrt{7}}
2376 !ang_phipphj(7,2,5)=\frac{1}{\sqrt{15}}
2377 !ang_phipphj(7,12,5)=\frac{13}{2 \sqrt{210}}
2378 !ang_phipphj(8,11,5)=\frac{1}{2 \sqrt{7}}
2379 !ang_phipphj(9,2,5)=-\frac{1}{4 \sqrt{5}}
2380 !ang_phipphj(9,10,5)=\frac{\sqrt{\frac{3}{14}}}{2}
2381 !ang_phipphj(9,12,5)=2 \sqrt{\frac{2}{35}}
2382 !ang_phipphj(10,9,5)=\frac{1}{\sqrt{42}}
2383 !ang_phipphj(11,8,5)=-\frac{1}{4 \sqrt{7}}
2384 !ang_phipphj(12,7,5)=-2 \sqrt{\frac{6}{35}}
2385 !ang_phipphj(12,9,5)=-\sqrt{\frac{2}{35}}
2386 !ang_phipphj(13,6,5)=2 \sqrt{\frac{3}{35}}
2387 !ang_phipphj(14,5,5)=\sqrt{\frac{2}{35}}
2388 !ang_phipphj(15,6,5)=\frac{1}{4 \sqrt{7}}
2389 !ang_phipphj(16,5,5)=-\frac{1}{\sqrt{42}}
2390 !ang_phipphj(1,2,6)=\frac{\sqrt{3}}{2}
2391 !ang_phipphj(1,12,6)=\frac{\sqrt{\frac{7}{6}}}{2}
2392 !ang_phipphj(2,9,6)=-\frac{\sqrt{5}}{2}
2393 !ang_phipphj(3,6,6)=\frac{\sqrt{5}}{2}
2394 !ang_phipphj(4,5,6)=\frac{\sqrt{5}}{2}
2395 !ang_phipphj(5,4,6)=-\frac{\sqrt{5}}{4}
2396 !ang_phipphj(5,16,6)=-\frac{\sqrt{\frac{21}{2}}}{2}
2397 !ang_phipphj(6,15,6)=-\frac{\sqrt{7}}{2}
2398 !ang_phipphj(7,12,6)=\frac{\sqrt{\frac{35}{6}}}{2}
2399 !ang_phipphj(8,11,6)=\frac{\sqrt{7}}{2}
2400 !ang_phipphj(9,2,6)=\frac{\sqrt{5}}{4}
2401 !ang_phipphj(9,10,6)=\frac{\sqrt{\frac{21}{2}}}{2}
2402 !ang_phipphj(10,9,6)=-\sqrt{\frac{7}{6}}
2403 !ang_phipphj(11,8,6)=-\frac{\sqrt{7}}{4}
2404 !ang_phipphj(15,6,6)=\frac{\sqrt{7}}{4}
2405 !ang_phipphj(16,5,6)=\sqrt{\frac{7}{6}}
2406 !ang_phipphj(1,3,7)=\frac{1}{\sqrt{3}}
2407 !ang_phipphj(2,6,7)=\frac{1}{\sqrt{5}}
2408 !ang_phipphj(3,1,7)=\frac{1}{\sqrt{3}}
2409 !ang_phipphj(3,7,7)=\frac{2}{\sqrt{15}}
2410 !ang_phipphj(4,8,7)=\frac{1}{\sqrt{5}}
2411 !ang_phipphj(5,11,7)=\frac{1}{\sqrt{7}}
2412 !ang_phipphj(6,2,7)=\frac{1}{\sqrt{5}}
2413 !ang_phipphj(6,12,7)=2 \sqrt{\frac{2}{35}}
2414 !ang_phipphj(7,3,7)=\frac{2}{\sqrt{15}}
2415 !ang_phipphj(7,13,7)=\frac{3}{\sqrt{35}}
2416 !ang_phipphj(8,4,7)=\frac{1}{\sqrt{5}}
2417 !ang_phipphj(8,14,7)=2 \sqrt{\frac{2}{35}}
2418 !ang_phipphj(9,15,7)=\frac{1}{\sqrt{7}}
2419 !ang_phipphj(11,5,7)=\frac{1}{\sqrt{7}}
2420 !ang_phipphj(12,6,7)=2 \sqrt{\frac{2}{35}}
2421 !ang_phipphj(13,7,7)=\frac{3}{\sqrt{35}}
2422 !ang_phipphj(14,8,7)=2 \sqrt{\frac{2}{35}}
2423 !ang_phipphj(15,9,7)=\frac{1}{\sqrt{7}}
2424 !ang_phipphj(1,3,8)=\frac{2}{\sqrt{3}}
2425 !ang_phipphj(2,6,8)=\frac{3}{\sqrt{5}}
2426 !ang_phipphj(3,7,8)=2 \sqrt{\frac{3}{5}}
2427 !ang_phipphj(4,8,8)=\frac{3}{\sqrt{5}}
2428 !ang_phipphj(5,11,8)=\frac{4}{\sqrt{7}}
2429 !ang_phipphj(6,2,8)=-\frac{1}{\sqrt{5}}
2430 !ang_phipphj(6,12,8)=8 \sqrt{\frac{2}{35}}
2431 !ang_phipphj(7,3,8)=-\frac{2}{\sqrt{15}}
2432 !ang_phipphj(7,13,8)=\frac{12}{\sqrt{35}}
2433 !ang_phipphj(8,4,8)=-\frac{1}{\sqrt{5}}
2434 !ang_phipphj(8,14,8)=8 \sqrt{\frac{2}{35}}
2435 !ang_phipphj(9,15,8)=\frac{4}{\sqrt{7}}
2436 !ang_phipphj(11,5,8)=-\frac{2}{\sqrt{7}}
2437 !ang_phipphj(12,6,8)=-4 \sqrt{\frac{2}{35}}
2438 !ang_phipphj(13,7,8)=-\frac{6}{\sqrt{35}}
2439 !ang_phipphj(14,8,8)=-4 \sqrt{\frac{2}{35}}
2440 !ang_phipphj(15,9,8)=-\frac{2}{\sqrt{7}}
2441 
2442 
2443  ang_phipphj_tmp=zero
2444 !
2445  ang_phipphj_tmp(1,4,1)=0.57735026918962576451_dp
2446  ang_phipphj_tmp(2,5,1)=0.44721359549995793928_dp
2447  ang_phipphj_tmp(3,8,1)=0.44721359549995793928_dp
2448  ang_phipphj_tmp(4,1,1)=0.57735026918962576451_dp
2449  ang_phipphj_tmp(4,7,1)=-0.25819888974716112568_dp
2450  ang_phipphj_tmp(4,9,1)=0.44721359549995793928_dp
2451  ang_phipphj_tmp(5,2,1)=0.44721359549995793928_dp
2452  ang_phipphj_tmp(5,10,1)=0.46291004988627573078_dp
2453  ang_phipphj_tmp(5,12,1)=-0.11952286093343936400_dp
2454  ang_phipphj_tmp(6,11,1)=0.37796447300922722721_dp
2455  ang_phipphj_tmp(7,4,1)=-0.25819888974716112568_dp
2456  ang_phipphj_tmp(7,14,1)=0.41403933560541253068_dp
2457  ang_phipphj_tmp(8,3,1)=0.44721359549995793928_dp
2458  ang_phipphj_tmp(8,13,1)=-0.29277002188455995381_dp
2459  ang_phipphj_tmp(8,15,1)=0.37796447300922722721_dp
2460  ang_phipphj_tmp(9,4,1)=0.44721359549995793928_dp
2461  ang_phipphj_tmp(9,14,1)=-0.11952286093343936400_dp
2462  ang_phipphj_tmp(9,16,1)=0.46291004988627573078_dp
2463  ang_phipphj_tmp(10,5,1)=0.46291004988627573078_dp
2464  ang_phipphj_tmp(11,6,1)=0.37796447300922722721_dp
2465  ang_phipphj_tmp(12,5,1)=-0.11952286093343936400_dp
2466  ang_phipphj_tmp(13,8,1)=-0.29277002188455995381_dp
2467  ang_phipphj_tmp(14,7,1)=0.41403933560541253068_dp
2468  ang_phipphj_tmp(14,9,1)=-0.11952286093343936400_dp
2469  ang_phipphj_tmp(15,8,1)=0.37796447300922722721_dp
2470  ang_phipphj_tmp(16,9,1)=0.46291004988627573078_dp
2471 !
2472  ang_phipphj_tmp(1,4,2)=0.28867513459481288225_dp
2473  ang_phipphj_tmp(1,14,2)=-0.54006172486732168591_dp
2474  ang_phipphj_tmp(2,5,2)=0.22360679774997896964_dp
2475  ang_phipphj_tmp(3,8,2)=0.22360679774997896964_dp
2476  ang_phipphj_tmp(4,7,2)=-0.77459666924148337704_dp
2477  ang_phipphj_tmp(4,9,2)=0.22360679774997896964_dp
2478  ang_phipphj_tmp(5,2,2)=0.11180339887498948482_dp
2479  ang_phipphj_tmp(5,10,2)=0.23145502494313786539_dp
2480  ang_phipphj_tmp(5,12,2)=-0.47809144373375745599_dp
2481  ang_phipphj_tmp(6,11,2)=0.18898223650461361361_dp
2482  ang_phipphj_tmp(7,4,2)=0.25819888974716112568_dp
2483  ang_phipphj_tmp(7,14,2)=0.44854261357253024157_dp
2484  ang_phipphj_tmp(8,3,2)=-0.44721359549995793928_dp
2485  ang_phipphj_tmp(8,13,2)=-1.1710800875382398152_dp
2486  ang_phipphj_tmp(8,15,2)=0.18898223650461361361_dp
2487  ang_phipphj_tmp(9,4,2)=0.11180339887498948482_dp
2488  ang_phipphj_tmp(9,14,2)=-0.47809144373375745599_dp
2489  ang_phipphj_tmp(9,16,2)=0.23145502494313786539_dp
2490  ang_phipphj_tmp(10,5,2)=0.15430334996209191026_dp
2491  ang_phipphj_tmp(11,6,2)=-0.094491118252306806804_dp
2492  ang_phipphj_tmp(12,5,2)=0.23904572186687872799_dp
2493  ang_phipphj_tmp(13,8,2)=0.58554004376911990761_dp
2494  ang_phipphj_tmp(14,7,2)=-0.82807867121082506136_dp
2495  ang_phipphj_tmp(14,9,2)=0.23904572186687872799_dp
2496  ang_phipphj_tmp(15,8,2)=-0.094491118252306806804_dp
2497  ang_phipphj_tmp(16,9,2)=0.15430334996209191026_dp
2498 !
2499  ang_phipphj_tmp(1,4,3)=0.86602540378443864676_dp
2500  ang_phipphj_tmp(1,14,3)=0.54006172486732168591_dp
2501  ang_phipphj_tmp(2,5,3)=1.1180339887498948482_dp
2502  ang_phipphj_tmp(3,8,3)=1.1180339887498948482_dp
2503  ang_phipphj_tmp(4,9,3)=1.1180339887498948482_dp
2504  ang_phipphj_tmp(5,2,3)=-0.55901699437494742410_dp
2505  ang_phipphj_tmp(5,10,3)=1.6201851746019650577_dp
2506  ang_phipphj_tmp(6,11,3)=1.3228756555322952953_dp
2507  ang_phipphj_tmp(7,14,3)=1.2076147288491198811_dp
2508  ang_phipphj_tmp(8,15,3)=1.3228756555322952953_dp
2509  ang_phipphj_tmp(9,4,3)=-0.55901699437494742410_dp
2510  ang_phipphj_tmp(9,16,3)=1.6201851746019650577_dp
2511  ang_phipphj_tmp(10,5,3)=-1.0801234497346433718_dp
2512  ang_phipphj_tmp(11,6,3)=-0.66143782776614764763_dp
2513  ang_phipphj_tmp(15,8,3)=-0.66143782776614764763_dp
2514  ang_phipphj_tmp(16,9,3)=-1.0801234497346433718_dp
2515 !
2516  ang_phipphj_tmp(1,2,4)=0.57735026918962576451_dp
2517  ang_phipphj_tmp(2,1,4)=0.57735026918962576451_dp
2518  ang_phipphj_tmp(2,7,4)=-0.25819888974716112568_dp
2519  ang_phipphj_tmp(2,9,4)=-0.44721359549995793928_dp
2520  ang_phipphj_tmp(3,6,4)=0.44721359549995793928_dp
2521  ang_phipphj_tmp(4,5,4)=0.44721359549995793928_dp
2522  ang_phipphj_tmp(5,4,4)=0.44721359549995793928_dp
2523  ang_phipphj_tmp(5,14,4)=-0.11952286093343936400_dp
2524  ang_phipphj_tmp(5,16,4)=-0.46291004988627573078_dp
2525  ang_phipphj_tmp(6,3,4)=0.44721359549995793928_dp
2526  ang_phipphj_tmp(6,13,4)=-0.29277002188455995381_dp
2527  ang_phipphj_tmp(6,15,4)=-0.37796447300922722721_dp
2528  ang_phipphj_tmp(7,2,4)=-0.25819888974716112568_dp
2529  ang_phipphj_tmp(7,12,4)=0.41403933560541253068_dp
2530  ang_phipphj_tmp(8,11,4)=0.37796447300922722721_dp
2531  ang_phipphj_tmp(9,2,4)=-0.44721359549995793928_dp
2532  ang_phipphj_tmp(9,10,4)=0.46291004988627573078_dp
2533  ang_phipphj_tmp(9,12,4)=0.11952286093343936400_dp
2534  ang_phipphj_tmp(10,9,4)=0.46291004988627573078_dp
2535  ang_phipphj_tmp(11,8,4)=0.37796447300922722721_dp
2536  ang_phipphj_tmp(12,7,4)=0.41403933560541253068_dp
2537  ang_phipphj_tmp(12,9,4)=0.11952286093343936400_dp
2538  ang_phipphj_tmp(13,6,4)=-0.29277002188455995381_dp
2539  ang_phipphj_tmp(14,5,4)=-0.11952286093343936400_dp
2540  ang_phipphj_tmp(15,6,4)=-0.37796447300922722721_dp
2541  ang_phipphj_tmp(16,5,4)=-0.46291004988627573078_dp
2542 !
2543  ang_phipphj_tmp(1,2,5)=0.28867513459481288225_dp
2544  ang_phipphj_tmp(1,12,5)=-0.54006172486732168591_dp
2545  ang_phipphj_tmp(2,7,5)=-0.77459666924148337704_dp
2546  ang_phipphj_tmp(2,9,5)=-0.22360679774997896964_dp
2547  ang_phipphj_tmp(3,6,5)=0.22360679774997896964_dp
2548  ang_phipphj_tmp(4,5,5)=0.22360679774997896964_dp
2549  ang_phipphj_tmp(5,4,5)=0.11180339887498948482_dp
2550  ang_phipphj_tmp(5,14,5)=-0.47809144373375745599_dp
2551  ang_phipphj_tmp(5,16,5)=-0.23145502494313786539_dp
2552  ang_phipphj_tmp(6,3,5)=-0.44721359549995793928_dp
2553  ang_phipphj_tmp(6,13,5)=-1.1710800875382398152_dp
2554  ang_phipphj_tmp(6,15,5)=-0.18898223650461361361_dp
2555  ang_phipphj_tmp(7,2,5)=0.25819888974716112568_dp
2556  ang_phipphj_tmp(7,12,5)=0.44854261357253024157_dp
2557  ang_phipphj_tmp(8,11,5)=0.18898223650461361361_dp
2558  ang_phipphj_tmp(9,2,5)=-0.11180339887498948482_dp
2559  ang_phipphj_tmp(9,10,5)=0.23145502494313786539_dp
2560  ang_phipphj_tmp(9,12,5)=0.47809144373375745599_dp
2561  ang_phipphj_tmp(10,9,5)=0.15430334996209191026_dp
2562  ang_phipphj_tmp(11,8,5)=-0.094491118252306806804_dp
2563  ang_phipphj_tmp(12,7,5)=-0.82807867121082506136_dp
2564  ang_phipphj_tmp(12,9,5)=-0.23904572186687872799_dp
2565  ang_phipphj_tmp(13,6,5)=0.58554004376911990761_dp
2566  ang_phipphj_tmp(14,5,5)=0.23904572186687872799_dp
2567  ang_phipphj_tmp(15,6,5)=0.094491118252306806804_dp
2568  ang_phipphj_tmp(16,5,5)=-0.15430334996209191026_dp
2569 !
2570  ang_phipphj_tmp(1,2,6)=0.86602540378443864676_dp
2571  ang_phipphj_tmp(1,12,6)=0.54006172486732168591_dp
2572  ang_phipphj_tmp(2,9,6)=-1.1180339887498948482_dp
2573  ang_phipphj_tmp(3,6,6)=1.1180339887498948482_dp
2574  ang_phipphj_tmp(4,5,6)=1.1180339887498948482_dp
2575  ang_phipphj_tmp(5,4,6)=-0.55901699437494742410_dp
2576  ang_phipphj_tmp(5,16,6)=-1.6201851746019650577_dp
2577  ang_phipphj_tmp(6,15,6)=-1.3228756555322952953_dp
2578  ang_phipphj_tmp(7,12,6)=1.2076147288491198811_dp
2579  ang_phipphj_tmp(8,11,6)=1.3228756555322952953_dp
2580  ang_phipphj_tmp(9,2,6)=0.55901699437494742410_dp
2581  ang_phipphj_tmp(9,10,6)=1.6201851746019650577_dp
2582  ang_phipphj_tmp(10,9,6)=-1.0801234497346433718_dp
2583  ang_phipphj_tmp(11,8,6)=-0.66143782776614764763_dp
2584  ang_phipphj_tmp(15,6,6)=0.66143782776614764763_dp
2585  ang_phipphj_tmp(16,5,6)=1.0801234497346433718_dp
2586 !
2587  ang_phipphj_tmp(1,3,7)=0.57735026918962576451_dp
2588  ang_phipphj_tmp(2,6,7)=0.44721359549995793928_dp
2589  ang_phipphj_tmp(3,1,7)=0.57735026918962576451_dp
2590  ang_phipphj_tmp(3,7,7)=0.51639777949432225136_dp
2591  ang_phipphj_tmp(4,8,7)=0.44721359549995793928_dp
2592  ang_phipphj_tmp(5,11,7)=0.37796447300922722721_dp
2593  ang_phipphj_tmp(6,2,7)=0.44721359549995793928_dp
2594  ang_phipphj_tmp(6,12,7)=0.47809144373375745599_dp
2595  ang_phipphj_tmp(7,3,7)=0.51639777949432225136_dp
2596  ang_phipphj_tmp(7,13,7)=0.50709255283710994651_dp
2597  ang_phipphj_tmp(8,4,7)=0.44721359549995793928_dp
2598  ang_phipphj_tmp(8,14,7)=0.47809144373375745599_dp
2599  ang_phipphj_tmp(9,15,7)=0.37796447300922722721_dp
2600  ang_phipphj_tmp(11,5,7)=0.37796447300922722721_dp
2601  ang_phipphj_tmp(12,6,7)=0.47809144373375745599_dp
2602  ang_phipphj_tmp(13,7,7)=0.50709255283710994651_dp
2603  ang_phipphj_tmp(14,8,7)=0.47809144373375745599_dp
2604  ang_phipphj_tmp(15,9,7)=0.37796447300922722721_dp
2605 !
2606  ang_phipphj_tmp(1,3,8)=1.1547005383792515290_dp
2607  ang_phipphj_tmp(2,6,8)=1.3416407864998738178_dp
2608  ang_phipphj_tmp(3,7,8)=1.5491933384829667541_dp
2609  ang_phipphj_tmp(4,8,8)=1.3416407864998738178_dp
2610  ang_phipphj_tmp(5,11,8)=1.5118578920369089089_dp
2611  ang_phipphj_tmp(6,2,8)=-0.44721359549995793928_dp
2612  ang_phipphj_tmp(6,12,8)=1.9123657749350298240_dp
2613  ang_phipphj_tmp(7,3,8)=-0.51639777949432225136_dp
2614  ang_phipphj_tmp(7,13,8)=2.0283702113484397860_dp
2615  ang_phipphj_tmp(8,4,8)=-0.44721359549995793928_dp
2616  ang_phipphj_tmp(8,14,8)=1.9123657749350298240_dp
2617  ang_phipphj_tmp(9,15,8)=1.5118578920369089089_dp
2618  ang_phipphj_tmp(11,5,8)=-0.75592894601845445443_dp
2619  ang_phipphj_tmp(12,6,8)=-0.95618288746751491198_dp
2620  ang_phipphj_tmp(13,7,8)=-1.0141851056742198930_dp
2621  ang_phipphj_tmp(14,8,8)=-0.95618288746751491198_dp
2622  ang_phipphj_tmp(15,9,8)=-0.75592894601845445443_dp
2623 
2624  ang_phipphj(:,:,:)=ang_phipphj_tmp(1:mpsang**2,1:mpsang**2,:)
2625 
2626  end subroutine setnabla_ylm

m_paw_sphharm/setsym_ylm [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 setsym_ylm

FUNCTION

 Compute rotation matrices expressed in the basis of real spherical harmonics
 This coefficients are used later to symmetrize PAW on-site quantities (rhoij, dij, ...).

INPUTS

  gprimd(3,3)==dimensional primitive translations for reciprocal space (bohr^-1)
  lmax=value of lmax mentioned at the second line of the psp file
  nsym=number of symmetry elements in space group
  pawprtvol=control print volume and debugging output
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  sym(3,3,nsym)=symmetries of group in terms of operations on primitive translations

OUTPUT

  zarot(2*lmax+1,2*lmax+1,lmax+1,nsym)=coefficients of the
      transformation of real spherical harmonics
      under the symmetry operations

NOTES

  Typical use: sym(:,:,:) is symrec(:,:,:) (rotations in reciprocal space)
               because we need symrel^-1 (=transpose[symrec])
               to symmetrize quantities.

  - This file comes from the file crystal_symmetry.f
    by N.A.W. Holzwarth and A. Tackett for the code pwpaw
  - Uses sign & phase convension of  M. E. Rose, Elementary Theory of Angular
    Momentum, John Wiley & Sons,. inc. 1957)
    zalpha = exp(-i*alpha)   zgamma = exp (-i*gamma)
  - Assumes each transformation  can be expressed in terms of 3 Euler
    angles with or without inversion

  Reference for evaluation of rotation matrices in the basis of real SH:
  Blanco M.A., Florez M. and Bermejo M.
  Journal of Molecular Structure: THEOCHEM, Volume 419, Number 1, 8 December 1997 , pp. 19-27(9)
  http://www.unioviedo.es/qcg/art/Theochem419-19-ov-BF97-rotation-matrices.pdf

SOURCE

2102 subroutine setsym_ylm(gprimd,lmax,nsym,pawprtvol,rprimd,sym,zarot)
2103 
2104 !Arguments ---------------------------------------------
2105 !scalars
2106  integer,intent(in) :: lmax,nsym,pawprtvol
2107 !arrays
2108  integer,intent(in) :: sym(3,3,nsym)
2109  real(dp),intent(in) :: gprimd(3,3),rprimd(3,3)
2110  real(dp),intent(out) :: zarot(2*lmax+1,2*lmax+1,lmax+1,nsym)
2111 
2112 !Local variables ------------------------------
2113 !scalars
2114  integer :: i1,ii,il,irot,isn,j1,jj,k1,ll,mm,mp
2115  real(dp) :: cosalp,cosbeta,cosgam,sinalp,sinbeta,singam
2116  character(len=1000) :: msg
2117 !arrays
2118  real(dp) :: prod(3,3),rot(3,3)
2119 
2120 !************************************************************************
2121 
2122  if (abs(pawprtvol)>=3) then
2123    write(msg,'(8a,i4)') ch10,&
2124 &   ' PAW TEST:',ch10,&
2125 &   ' ==== setsym_ylm: rotation matrices in the basis ============',ch10,&
2126 &   ' ====              of real spherical harmonics    ============',ch10,&
2127 &   '  > Number of symmetries (nsym)=',nsym
2128    call wrtout(std_out,msg,'COLL')
2129  end if
2130 
2131  zarot=zero
2132 
2133  do irot=1,nsym
2134 
2135    if (abs(pawprtvol)>=3) then
2136      write(msg,'(a,i2,a,9i2,a)') '   >For symmetry ',irot,' (',sym(:,:,irot),')'
2137      call wrtout(std_out,msg,'COLL')
2138    end if
2139 
2140 !  === l=0 case ===
2141    zarot(1,1,1,irot)=one
2142 
2143 !  === l>0 case ===
2144    if (lmax>0) then
2145 !    Calculate the rotations in the cartesian basis
2146      rot=zero;prod=zero
2147      do k1=1,3
2148        do j1=1,3
2149          do i1=1,3
2150            prod(i1,j1)=prod(i1,j1)+sym(i1,k1,irot)*rprimd(j1,k1)
2151          end do
2152        end do
2153      end do
2154      do j1=1,3
2155        do i1=1,3
2156          do k1=1,3
2157            rot(i1,j1)=rot(i1,j1)+gprimd(i1,k1)*prod(k1,j1)
2158          end do
2159          if(abs(rot(i1,j1))<tol10) rot(i1,j1)=zero
2160        end do
2161      end do
2162      call mkeuler(rot,cosbeta,sinbeta,cosalp,sinalp,cosgam,singam,isn)
2163      do ll=1,lmax
2164        il=(isn)**ll
2165        do mp=-ll,ll
2166          jj=mp+ll+1
2167          do mm=-ll,ll
2168            ii=mm+ll+1
2169 
2170 !          Formula (47) from the paper of Blanco et al
2171            zarot(ii,jj,ll+1,irot)=il&
2172 &           *(phim(cosalp,sinalp,mm)*phim(cosgam,singam,mp)*sign(1,mp)&
2173            *(dbeta(cosbeta,sinbeta,ll,abs(mp),abs(mm))&
2174 &           +(-1._dp)**mm*dbeta(cosbeta,sinbeta,ll,abs(mm),-abs(mp)))*half&
2175 &           -phim(cosalp,sinalp,-mm)*phim(cosgam,singam,-mp)*sign(1,mm)&
2176            *(dbeta(cosbeta,sinbeta,ll,abs(mp),abs(mm))&
2177 &           -(-1._dp)**mm*dbeta(cosbeta,sinbeta,ll,abs(mm),-abs(mp)))*half)
2178          end do
2179        end do
2180      end do
2181    end if   ! lmax case
2182 
2183    if (abs(pawprtvol)>=3) then
2184      if(lmax>0) then
2185        write(msg,'(2a,3(3(2x,f7.3),a))') &
2186 &       '    Rotation matrice for l=1:',ch10,&
2187 &       (zarot(1,jj,2,irot),jj=1,3),ch10,&
2188 &       (zarot(2,jj,2,irot),jj=1,3),ch10,&
2189 &       (zarot(3,jj,2,irot),jj=1,3)
2190        call wrtout(std_out,msg,'COLL')
2191      end if
2192      if(lmax>1) then
2193        write(msg,'(2a,5(5(2x,f7.3),a))') &
2194 &       '    Rotation matrice for l=2:',ch10,&
2195 &       (zarot(1,jj,3,irot),jj=1,5),ch10,&
2196 &       (zarot(2,jj,3,irot),jj=1,5),ch10,&
2197 &       (zarot(3,jj,3,irot),jj=1,5),ch10,&
2198 &       (zarot(4,jj,3,irot),jj=1,5),ch10,&
2199 &       (zarot(5,jj,3,irot),jj=1,5)
2200        call wrtout(std_out,msg,'COLL')
2201      end if
2202      if(lmax>2) then
2203        write(msg,'(2a,7(7(2x,f7.3),a))') &
2204 &       '    Rotation matrice for l=3:',ch10,&
2205 &       (zarot(1,jj,4,irot),jj=1,7),ch10,&
2206 &       (zarot(2,jj,4,irot),jj=1,7),ch10,&
2207 &       (zarot(3,jj,4,irot),jj=1,7),ch10,&
2208 &       (zarot(4,jj,4,irot),jj=1,7),ch10,&
2209 &       (zarot(5,jj,4,irot),jj=1,7),ch10,&
2210 &       (zarot(6,jj,4,irot),jj=1,7),ch10,&
2211 &       (zarot(7,jj,4,irot),jj=1,7)
2212        call wrtout(std_out,msg,'COLL')
2213      end if
2214    end if
2215 
2216  end do  ! isym loop
2217 
2218 end subroutine setsym_ylm

m_paw_sphharm/slxyzs [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 slxyzs

FUNCTION

 computes the matrix element <Sl'm'|L_idir|Slm>

INPUTS

   integer :: lp,mp,idir,ll,mm

OUTPUT

   complex(dpc) :: sls_val

NOTES

 Slm is the real spherical harmonic used througout abinit,
 L_idir is a component of the angular momentum operator.
 The subroutine computes <S_l'm'|L_idir|S_lm>

SOURCE

858 subroutine slxyzs(lp,mp,idir,ll,mm,sls_val)
859 
860 !Arguments ---------------------------------------------
861 !scalars
862  integer,intent(in) :: idir,ll,lp,mm,mp
863  complex(dpc),intent(out) :: sls_val
864 
865 !Local variables ---------------------------------------
866 !scalars
867  integer :: mpp,mppp
868  complex(dpc) :: lidir,sy_val,ys_val
869 
870 ! *********************************************************************
871 
872  sls_val = czero
873 
874  if ( lp /= ll ) return
875 
876  do mpp = -ll, ll
877    call ys(ll,mpp,ll,mp,sy_val)
878    do mppp = -ll, ll
879      call lxyz(ll,mpp,idir,ll,mppp,lidir)
880      call ys(ll,mppp,ll,mm,ys_val)
881      sls_val = sls_val + conjg(sy_val)*lidir*ys_val
882    end do
883  end do
884 
885 end subroutine slxyzs

m_paw_sphharm/ylm_angular_mesh [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 ylm_angular_mesh

FUNCTION

  Build (theta, phi) angular mesh from (ntheta, nphi)

INPUTS

   ntheta= number of sample points in the theta dir
   nphi= number of sample points in the phi dir

OUTPUT

   angl_size= total number of sample points in the angular mesh, i.e. (ntheta * nphi)
   cart_coord(3, angl_size)= for each point of the angular mesh, gives the Cartesian coordinates
     of the corresponding point on an unitary sphere.
   ang_wgth(angl_size)= for each point of the angular mesh, gives the weight
       of the corresponding point on an unitary sphere.

 NOTE
   Summing over f * angwgth gives the spherical average 1/(4pi) \int domega f(omega)

SOURCE

1582 subroutine ylm_angular_mesh(ntheta, nphi, angl_size, cart_coord, ang_wgth)
1583 
1584 !Arguments ------------------------------------
1585  integer,intent(in) :: ntheta, nphi
1586  integer,intent(out) :: angl_size
1587  real(dp),allocatable,intent(out) :: cart_coord(:,:)
1588  real(dp),allocatable,intent(out) :: ang_wgth(:)
1589 
1590 !Local variables ------------------------------
1591 !scalars
1592  integer :: it, ip, npoints
1593  real(dp) :: ang, con, cos_phi, cos_theta, sin_phi, sin_theta
1594  character(len=500) :: msg
1595 !arrays
1596  real(dp),allocatable :: th(:),wth(:)
1597 
1598 ! *************************************************************************
1599 
1600  LIBPAW_ALLOCATE(th, (ntheta))
1601  LIBPAW_ALLOCATE(wth, (ntheta))
1602 
1603  con = two_pi / nphi
1604  call gauleg(-one, one, th, wth, ntheta)
1605 
1606  angl_size = ntheta * nphi
1607  LIBPAW_ALLOCATE(cart_coord, (3, angl_size))
1608  LIBPAW_ALLOCATE(ang_wgth, (angl_size))
1609  npoints = 0
1610  do it = 1, ntheta
1611    cos_theta = th(it)
1612    sin_theta = sqrt(one - cos_theta*cos_theta)
1613    do ip = 1, nphi
1614      ang = con * (ip-1)
1615      cos_phi = cos(ang); sin_phi = sin(ang)
1616      npoints = npoints + 1
1617      cart_coord(1, npoints) = sin_theta * cos_phi
1618      cart_coord(2, npoints) = sin_theta * sin_phi
1619      cart_coord(3, npoints) = cos_theta
1620      ! Normalization required
1621      ang_wgth(npoints) = wth(it) / (two * nphi)
1622    end do
1623  end do
1624 
1625  LIBPAW_DEALLOCATE(th)
1626  LIBPAW_DEALLOCATE(wth)
1627 
1628 !Error if npoints exceeds angl_size
1629  if (npoints > angl_size) then
1630    write(msg, '(a,i4,a,a,i4)' ) 'npoints =',npoints,ch10,&
1631 &                               'angl_size =',angl_size
1632    LIBPAW_BUG(msg)
1633  end if
1634 
1635 end subroutine ylm_angular_mesh

m_paw_sphharm/ylm_cmplx [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 ylm_cmplx

FUNCTION

  Calculate all (complex) spherical harmonics for lx<=4

INPUTS

  lx= quantum numbers.
  xx= cartesian coordinate in the x direction
  yy= cartesian coordinate in the y direction
  zz= cartesian coordinate in the z direction

 cartesian coordinates

OUTPUT

  ylm((lx+1)*(lx+1)) complex spherical harmonics for all l<=lx and all
                     possible values of m.

NOTES

  We are supressing the so-called Condon-Shortley phase

SOURCE

393 subroutine ylm_cmplx(lx,ylm,xx,yy,zz)
394 
395 !Arguments ------------------------------------
396 !scalars
397  integer,intent(in) :: lx
398  real(dp),intent(in) :: xx,yy,zz
399 !arrays
400  complex(dpc),intent(out) :: ylm((lx+1)*(lx+1))
401 
402 !Local variables-------------------------------
403 !scalars
404  integer :: ii,l1,m1,nc,nn
405  real(dp) :: dc,dl,dm,ds,rr,rrs,rs,sq2,w,x,xs,ya,yi,yr
406 !arrays
407  real(dp) :: cosa(lx+1),fact(2*(lx+1)),plm(lx+2,lx+2),qlm(lx+2,lx+2),sgn(lx+1)
408  real(dp) :: sina(lx+1)
409 
410 ! *************************************************************************
411 
412 !normalization coefficients
413  sq2=sqrt(2.0d0)
414  fact(1)=1.0d0
415  do ii=2,2*lx+1
416    fact(ii)=(ii-1)*fact(ii-1)
417  end do
418  do l1=1,lx+1
419    sgn(l1)=(-1.d0)**(l1-1)
420    do m1=1,l1
421      qlm(l1,m1)=sqrt((2*l1-1)*fact(l1-m1+1)/&
422 &     (four_pi*fact(l1+m1-1)))
423    end do
424  end do
425 
426 !legendre polynomials
427  rs=xx**2 + yy**2 + zz**2
428  if(rs > tol8) then
429    xs=zz**2/rs
430    x=zz/sqrt(rs)
431    w=sqrt(abs(1.0d0 - xs))
432  else
433    x=0.0d0
434 
435    w=1.0d0
436  end if
437  plm(1,1)=1.0d0
438  plm(2,1)=x
439  plm(2,2)=w
440  plm(3,2)=3.0d0*x*w
441  do m1=1,lx
442    dm=m1-1
443    if(m1 > 1) then
444      plm(m1+1,m1)=x*plm(m1,m1) + 2*dm*w*plm(m1,m1-1)
445    end if
446    if(m1 < lx) then
447      do l1=m1+2,lx+1
448        dl=l1-1
449        plm(l1,m1)=((2*dl-1)*x*plm(l1-1,m1)&
450 &       - (dl+dm-1)*plm(l1-2,m1))/(dl-dm)
451      end do
452    end if
453    plm(m1+1,m1+1)=(2*dm+1)*w*plm(m1,m1)
454  end do
455 
456 !azimuthal angle phase factors
457  rrs=xx**2 + yy**2
458  if(rrs > tol8) then
459    rr=sqrt(rrs)
460    dc=xx/rr
461    ds=yy/rr
462  else
463    dc=1.0d0
464    ds=0.0d0
465  end if
466  cosa(1)=1.0d0
467  sina(1)=0.0d0
468  do m1=2,lx+1
469    cosa(m1)=dc*cosa(m1-1) - ds*sina(m1-1)
470    sina(m1)=ds*cosa(m1-1) + dc*sina(m1-1)
471  end do
472 
473 !combine factors
474  do l1=1,lx+1
475    do m1=2,l1
476      nn=(l1-1)**2 + (l1-1) + (m1-1) + 1
477      nc=(l1-1)**2 + (l1-1) - (m1-1) + 1
478 !    note that we are supressing the so-called Condon-Shortley phase
479 !    ya=sgn(m1)*qlm(l1,m1)*plm(l1,m1)
480      ya=qlm(l1,m1)*plm(l1,m1)
481      yr=ya*cosa(m1)
482      yi=ya*sina(m1)
483      ylm(nc)=sgn(m1)*cmplx(yr,-yi)
484      ylm(nn)=cmplx(yr,yi)
485    end do
486  end do
487  do l1=1,lx+1
488    nn=(l1-1)**2 + (l1-1) + 1
489    ya=qlm(l1,1)*plm(l1,1)
490    ylm(nn)=cmplx(ya,0.d0)
491  end do
492 
493 end subroutine ylm_cmplx

m_paw_sphharm/ylmc [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 ylmc

FUNCTION

  Return a complex spherical harmonic with l <= 3

INPUTS

  il=angular quantum number
  im=magnetic quantum number
  kcart=vector in cartesian coordinates defining the value of \theta and \psi
   where calculate the spherical harmonic

OUTPUT

  ylm= spherical harmonic

NOTES

  Note the use of double precision complex.
  Case l>3 not implemented.

SOURCE

 95 function ylmc(il,im,kcart)
 96 
 97 !Arguments ------------------------------------
 98 !scalars
 99  integer,intent(in) :: il,im
100  complex(dpc) :: ylmc
101 !arrays
102  real(dp),intent(in) :: kcart(3)
103 
104 !Local variables-------------------------------
105 !scalars
106  integer,parameter :: LMAX=3
107  real(dp),parameter :: PPAD=tol8
108  real(dp) :: cosphi,costh,costhreephi,costwophi,r,rxy,sinphi,sinth,sinthreephi,sintwophi
109  !complex(dpc) :: new_ylmc
110  character(len=500) :: msg
111  complex(dpc) :: ctmp
112 
113 ! *************************************************************************
114 
115  if (ABS(im)>ABS(il)) then
116    write(msg,'(3(a,i0))') 'm is,',im,' however it should be between ',-il,' and ',il
117    LIBPAW_ERROR(msg)
118  end if
119 
120  ylmc = czero
121 
122  r=SQRT(kcart(1)**2+kcart(2)**2+kcart(3)**2)
123  if (r<PPAD) r=r+PPAD
124  !$if (r<tol10) RETURN
125 
126  rxy=SQRT(kcart(1)**2+kcart(2)**2)
127  if (rxy<PPAD)rxy=r+PPAD
128 !
129 ! Determine theta and phi
130  costh= kcart(3)/r
131 
132 #if 1
133  ! old buggy coding
134  sinth= rxy/r
135  cosphi= kcart(1)/rxy
136  sinphi= kcart(2)/rxy
137 #else
138  sinth=sqrt(abs((one-costh)*(one+costh))) ! abs is needed to prevent very small negative arg
139  cosphi=one
140  sinphi=zero
141  if (sinth>tol10) then
142    cosphi=kcart(1)/(r*sinth)
143    sinphi=kcart(2)/(r*sinth)
144  end if
145 #endif
146 
147  costwophi= two*cosphi**2 - one
148  sintwophi= two*sinphi*cosphi
149  costhreephi=cosphi*costwophi-sinphi*sintwophi
150  sinthreephi=cosphi*sintwophi+sinphi*costwophi
151 
152  select case (il)
153 
154  case (0)
155   ylmc= one/SQRT(four_pi)
156 
157  case (1)
158   if (ABS(im)==0) then
159    ylmc = SQRT(three/(four_pi))*costh
160   else if (ABS(im)==1) then
161    ylmc = -SQRT(three/(8._dp*pi))*sinth*CMPLX(cosphi,sinphi)
162   else
163    msg='wrong im'
164    LIBPAW_ERROR(msg)
165   end if
166 
167  case (2)
168   if (ABS(im)==0) then
169    ylmc = SQRT(5.d0/(16.d0*pi))*(three*costh**2-one)
170   else if (ABS(im)==1) then
171    ylmc = -SQRT(15.d0/(8.d0*pi))*sinth*costh*cmplx(cosphi,sinphi)
172   else if (ABS(im)==2) then
173    ylmc = SQRT(15.d0/(32.d0*pi))*(sinth)**2*CMPLX(costwophi,sintwophi)
174   else
175    msg='wrong im'
176    LIBPAW_ERROR(msg)
177   end if
178 
179  case (3)
180   if (ABS(im)==0) then
181    ylmc= SQRT(7.d0/(16.d0*pi))*(5.d0*costh**3 -3.d0*costh)
182   else if (ABS(im)==1) then
183    ylmc= -SQRT(21.d0/(64.d0*pi))*sinth*(5.d0*costh**2-one)*CMPLX(cosphi,sinphi)
184   else if (ABS(im)==2) then
185    ylmc= SQRT(105.d0/(32.d0*pi))*sinth**2*costh*CMPLX(costwophi,sintwophi)
186   else if (ABS(im)==3) then
187    ylmc=-SQRT(35.d0/(64.d0*pi))*sinth**3*CMPLX(costhreephi,sinthreephi)
188   else
189    msg='wrong im'
190    LIBPAW_ERROR(msg)
191   end if
192 
193  case default
194   !write(msg,'(a,i6,a,i6)')' The maximum allowed value for l is,',LMAX,' however l=',il
195   !LIBPAW_ERROR(msg)
196  end select
197 !
198 !=== Treat the case im < 0 ===
199  if (im < 0) then
200    ctmp = (-one)**(im)*CONJG(ylmc)
201    ylmc = ctmp
202  end if
203 
204  ! FIXME: Use the piece of code below as it works for arbitrary (l,m)
205  ! the implementation above is buggy when the vector is along z!
206  !
207 #if 0
208 ! Remember the expression of complex spherical harmonics:
209 ! $Y_{lm}(\theta,\phi)=sqrt{{(2l+1) over (4\pi)} {fact(l-m)/fact(l+m)} } P_l^m(cos(\theta)) e^{i m\phi}$
210   new_ylmc = SQRT((2*il+1)*rfactorial(il-ABS(im))/(rfactorial(il+ABS(im))*four_pi)) * &
211 &   ass_leg_pol(il,ABS(im),costh) * CMPLX(cosphi,sinphi)**ABS(im)
212   if (im<0) new_ylmc=(-one)**(im)*CONJG(new_ylmc)
213 
214   if (ABS(new_ylmc-ylmc)>tol6) then
215     !LIBPAW_WARNING("Check new_ylmc")
216     !write(std_out,*)"il,im,new_ylmc, ylmc",il,im,new_ylmc,ylmc
217     !write(std_out,*)"fact",SQRT((2*il+1)*rfactorial(il-ABS(im))/(rfactorial(il+ABS(im))*four_pi))
218     !write(std_out,*)"costh,sinth,ass_leg_pol",costh,sinth,ass_leg_pol(il,ABS(im),costh)
219     !write(std_out,*)"cosphi,sinphi,e^{imphi}",cosphi,sinphi,CMPLX(cosphi,sinphi)**ABS(im)
220   end if
221   ylmc = new_ylmc
222 #endif
223 
224 end function ylmc

m_paw_sphharm/ylmcd [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 ylmcd

FUNCTION

  Computes dth and dphi, the first derivatives of complex Ylm as a function of
  th and phi (the angles of the spherical coordinates)
  It works for all spherical harmonics with l <= 3

INPUTS

  il=angular quantum number
  im=magnetic quantum number
  kcart=cartesian coordinates of the vector where the first derivatives of Ylm are evaluated

OUTPUT

  dth =derivative of Y_lm with respect to \theta
  dphi=derivative of Y_lm with respect to \phi

NOTES

  Note the use of double precision complex.
  Case l>3 not implemented.

SOURCE

253 subroutine ylmcd(il,im,kcart,dth,dphi)
254 
255 !Arguments ------------------------------------
256 !scalars
257  integer,intent(in) :: il,im
258  complex(dpc),intent(out) :: dphi,dth
259 !arrays
260  real(dp),intent(in) :: kcart(3)
261 
262 !Local variables-------------------------------
263 !scalars
264  integer,parameter :: LMAX=3
265  real(dp),parameter :: PPAD=tol8
266  real(dp) :: cosphi,costh,costhreephi,costwophi,r,rxy,sinphi,sinth,sinthreephi,sintwophi,c
267  character(len=500) :: msg
268  complex(dpc) :: ctmp 
269 
270 ! *************************************************************************
271 
272  if (ABS(im)>ABS(il))then
273    write(msg,'(3(a,i0))')' m is,',im,' however it should be between ',-il,' and ',il
274    LIBPAW_ERROR(msg)
275  end if
276 
277  dphi=czero; dth=czero
278 
279  r=SQRT(kcart(1)**2+kcart(2)**2+kcart(3)**2)
280  if (r<PPAD) r=r+PPAD
281  !$if (r<tol10) RETURN
282 
283  rxy=SQRT(kcart(1)**2+kcart(2)**2)
284  if (rxy<PPAD) rxy=r+PPAD
285 
286 ! Determine theta and phi
287  costh= kcart(3)/r
288 #if 1
289  ! old buggy coding
290  sinth= rxy/r
291  cosphi= kcart(1)/rxy
292  sinphi= kcart(2)/rxy
293 #else
294  sinth=sqrt(abs((one-costh)*(one+costh))) ! abs is needed to prevent very small negative arg
295  cosphi=one
296  sinphi=zero
297  if (sinth>tol10) then
298    cosphi=kcart(1)/(r*sinth)
299    sinphi=kcart(2)/(r*sinth)
300  end if
301 #endif
302 
303  costwophi= two*cosphi**2 - one
304  sintwophi= two*sinphi*cosphi
305  costhreephi=cosphi*costwophi-sinphi*sintwophi
306  sinthreephi=cosphi*sintwophi+sinphi*costwophi
307 
308  select case (il)
309 
310  case (0)
311    dth  = czero
312    dphi = czero
313 
314  case (1)
315    if (ABS(im)==0) then
316      dth= -SQRT(three/(four_pi))*sinth
317      dphi= czero
318    else if (abs(im)==1) then
319      dth= -SQRT(3.d0/(8.d0*pi))*costh*CMPLX(cosphi,sinphi)
320      dphi=-SQRT(3.d0/(8.d0*pi))*sinth*CMPLX(-sinphi,cosphi)
321    end if
322 
323  case (2)
324    if (ABS(im)==0) then
325      dth= -SQRT(5.d0/(16.d0*pi))*6.d0*costh*sinth
326      dphi= czero
327    else if (ABS(im)==1) then
328      dth=  -SQRT(15.d0/(8.d0*pi))*(costh**2-sinth**2)*CMPLX(cosphi,sinphi)
329      dphi= -SQRT(15.d0/(8.d0*pi))*costh*sinth*(0.d0,1.d0)*CMPLX(cosphi,sinphi)
330    else if (abs(im)==2) then
331      dth  = SQRT(15.d0/(32.d0*pi))*2.d0*costh*sinth*CMPLX(costwophi,sintwophi)
332      dphi = SQRT(15.d0/(32.d0*pi))*sinth**2*(0.d0,2.d0)*CMPLX(costwophi,sintwophi)
333    end if
334 
335  case (3)
336    if (ABS(im)==0) then
337      dth = SQRT(7.d0/(16*pi))*(-15.d0*costh**2*sinth + 3.d0**sinth)
338      dphi= czero
339    else if (ABS(im)==1) then
340      c = SQRT(21.d0/(64.d0*pi))
341      dth= -c*      (15.d0*costh**3-11.d0*costh)*            CMPLX(cosphi,sinphi)
342      dphi=-c*sinth*( 5.d0*costh**2-1          )*(0.d0,1.d0)*CMPLX(cosphi,sinphi)
343    else if (ABS(im)==2) then
344      c = SQRT(105.d0/(32.d0*pi))
345      dth =c*(2.d0*sinth*costh**2-sinth**3)   *CMPLX(costwophi,sintwophi)
346      dphi=c*(2.d0*sinth**2*costh)*(0.d0,1.d0)*CMPLX(costwophi,sintwophi)
347    else if (abs(im)==3) then
348      dth =-SQRT(35.d0/(64.d0*pi))*3.d0*sinth**2*costh*CMPLX(costhreephi,sinthreephi)
349      dphi=-SQRT(35.d0/(64.d0*pi))*sinth**3*(0.d0,3.d0)*CMPLX(costhreephi,sinthreephi)
350    end if
351 
352  case default
353    write(msg,'(2(a,i0))')' The maximum allowed value for l is,',LMAX,' however, l=',il
354    LIBPAW_ERROR(msg)
355  end select
356 !
357 !=== Treat the case im < 0 ===
358  if (im<0) then
359    ctmp = (-one)**(im)*CONJG(dth)
360    dth = ctmp
361    ctmp= (-one)**(im)*CONJG(dphi)
362    dphi= ctmp
363  end if
364 
365 end subroutine ylmcd

m_paw_sphharm/ys [ Functions ]

[ Top ] [ m_paw_sphharm ] [ Functions ]

NAME

 ys

FUNCTION

  Computes the matrix element <Y_(l2,m2)|S_(l1,m1)>

INPUTS

  integer :: l2,m2,l1,m1

OUTPUT

  complex(dpc) :: ys_val

NOTES

 Ylm is the standard complex-valued spherical harmonic, Slm is the real spherical harmonic
 used througout abinit. 

SOURCE

737 subroutine ys(l2,m2,l1,m1,ys_val)
738 
739 !Arguments ---------------------------------------------
740 !scalars
741  integer,intent(in) :: l1,l2,m1,m2
742  complex(dpc),intent(out) :: ys_val
743 
744 !Local variables ---------------------------------------
745  !scalars
746  integer :: mp1
747 
748 ! *********************************************************************
749 
750 ! See Blanco et al., J. Mol Struct. 419, 19-27 (1997) Eq. 19
751 ! <Y_l2,m2|S_l1,m1> is given by C^l_{m1,m2} where
752 ! l1 == l2 and |m1| == |m2|, 0 otherwise
753 
754  ys_val = czero
755 
756  if ( l2 /= l1 ) return
757  if ( abs(m2) /= abs(m1) ) return
758 
759  mp1=(-1)**abs(m1)
760  
761  if(m1.EQ.0) then
762    ys_val=cone
763  else if((m1.GT.0).AND.(m2.GT.0)) then
764    ys_val=mp1*sqrthalf
765  else if((m1.GT.0).AND.(m2.LT.0)) then
766    ys_val=sqrthalf
767  else if((m1.LT.0).AND.(m2.GT.0)) then
768    ys_val=-j_dpc*mp1*sqrthalf
769  else if((m1.LT.0).AND.(m2.LT.0)) then
770    ys_val=j_dpc*sqrthalf
771  else
772    ys_val=czero
773  end if
774 
775 end subroutine ys

m_pawsphharm/realgaunt [ Functions ]

[ Top ] [ Functions ]

NAME

 realgaunt

FUNCTION

 This routine compute "real Gaunt coefficients", i.e. gaunt
 coefficients according to "real spherical harmonics"
   RealGaunt(ilm,ilm_i,ilm_j) = Int[ S_lm Slm_i Slm_j]

INPUTS

  l_max= max. value of ang. momentum l+1;  Gaunt coeffs up to
          [(2*l_max-1,m),(l_max,m),(l_max,m)] are computed

OUTPUT

  gntselect((2*l_max-1)**2,l_max**2*(l_max**2+1)/2)=
          selection rules for Gaunt coefficients
          if Gaunt coeff. is zero, gntselect=0
          if Gaunt coeff. is non-zero, gntselect is the index of
                           the coeff. in realgnt(:) array
  ngnt= number of non-zero Gaunt coefficients
  realgnt((2*l_max-1)**2*l_max**4)= non-zero real Gaunt coefficients

 NOTE
 Second index of gntselect is in "upper triangle" format.
 Its formula is klm_ij = ilm_i*(ilm_i-1)/2 + ilm_j,
   corresponding to the two index pairs: (ilm_i,ilm_j) and (ilm_j,ilmj)

SOURCE

2772 subroutine realgaunt(l_max,ngnt,gntselect,realgnt)
2773 
2774 !Arguments ---------------------------------------------
2775 !scalars
2776  integer,intent(in) :: l_max
2777  integer,intent(out) :: ngnt
2778 !arrays
2779  integer,intent(out) :: gntselect(:,:)
2780  real(dp),intent(out) :: realgnt(:)
2781 
2782 !Local variables ------------------------------
2783 !scalars
2784  integer :: ilm1,ilm2,ilmp1,k0lm1,klm1,l1,l2,ll,lp1,m1,m2,mm,mm1,mm2,mm3,mp1
2785  real(dp) :: c11,c12,c21,c22,c31,c32,fact,realgnt_tmp
2786  character(len=500) :: msg
2787 !arrays
2788  integer,allocatable :: ssgn(:)
2789  type(coeff3_type), allocatable :: coeff(:)
2790 
2791 !************************************************************************
2792 
2793  if ( size(gntselect)<(2*l_max-1)**2*(l_max**2*(l_max**2+1))/2 .or. &
2794 &     size(realgnt)  <(2*l_max-1)**2*(l_max**2*(l_max**2+1))/2 ) then
2795    msg='Too small sizes for gntselect/realgnt!'
2796    LIBPAW_BUG(msg)
2797  end if
2798 
2799 !Initialize output arrays with zeros.
2800  gntselect = 0; realgnt = zero
2801 
2802 !Compute matrix cc where Sl=cc*Yl (Sl=real sph. harm.)
2803 !------------------------------------------------
2804  LIBPAW_DATATYPE_ALLOCATE(coeff,(4*l_max-3))
2805  do ll=1,4*l_max-3
2806    LIBPAW_ALLOCATE(coeff(ll)%value,(2,2*ll-1,2*ll-1))
2807    coeff(ll)%value(:,:,:)=zero
2808    coeff(ll)%value(1,ll,ll)=one
2809    do mm=1,ll-1
2810      coeff(ll)%value(1,ll+mm,ll+mm)= (-1._dp)**mm/sqrt(2._dp)
2811      coeff(ll)%value(1,ll-mm,ll+mm)= ( 1._dp)    /sqrt(2._dp)
2812      coeff(ll)%value(2,ll+mm,ll-mm)=-(-1._dp)**mm/sqrt(2._dp)
2813      coeff(ll)%value(2,ll-mm,ll-mm)= ( 1._dp)    /sqrt(2._dp)
2814    end do
2815  end do
2816 
2817  LIBPAW_ALLOCATE(ssgn,(l_max**2))
2818  ssgn(:)=1
2819  if (l_max>0) then
2820    do l1=1,l_max-1
2821      ilm1=1+l1**2+l1
2822      do m1=-l1,-1
2823        ssgn(ilm1+m1)=-1
2824      end do
2825    end do
2826  end if
2827 
2828  ngnt=0
2829 
2830 !Loop on (lp1,mp1)
2831 !------------------------------------------------
2832  do lp1=0,l_max-1
2833    do mp1=-lp1,lp1
2834      ilmp1=1+lp1**2+lp1+mp1
2835      k0lm1=ilmp1*(ilmp1-1)/2
2836 
2837 !    Loop on (l1,m1)<=(lp1,mp1)
2838 !    ------------------------------------------------
2839      do l1=0,l_max-1
2840        do m1=-l1,l1
2841          ilm1=1+l1**2+l1+m1
2842 
2843          if (ilm1<=ilmp1) then
2844 
2845            klm1=k0lm1+ilm1
2846            gntselect(:,klm1)=0
2847 
2848 !          Loop on (l2,m2)
2849 !          ------------------------------------------------
2850            do l2=abs(l1-lp1),l1+lp1,2
2851              do m2=-l2,l2
2852                ilm2=1+l2**2+l2+m2
2853 
2854 !              Real Gaunt coeffs selection rules
2855 !              ------------------------------------------------
2856                if ((l2<=l1+lp1).and.&
2857 &               (((m1== mp1).and.((m2==0).or.(m2==2*abs(mp1)))).or.&
2858 &               ((m1==-mp1).and.(m2==-abs(m1)-abs(mp1))).or.&
2859 &               ((abs(m1)/=(abs(mp1)).and.&
2860 &               ((m2==ssgn(ilm1)*ssgn(ilmp1)*   (abs(m1)+abs(mp1))).or.&
2861 &               (m2==ssgn(ilm1)*ssgn(ilmp1)*abs(abs(m1)-abs(mp1)))&
2862                ))))) then
2863 
2864 !                Compute selected real Gaunt coefficient
2865 !                ------------------------------------------------
2866                  realgnt_tmp=zero
2867                  do mm1=-l1,l1
2868                    c11=coeff(l1+1)%value(1,l1+mm1+1,l1+m1+1)
2869                    c12=coeff(l1+1)%value(2,l1+mm1+1,l1+m1+1)
2870                    do mm2= -lp1,lp1
2871                      c21=coeff(lp1+1)%value(1,lp1+mm2+1,lp1+mp1+1)
2872                      c22=coeff(lp1+1)%value(2,lp1+mm2+1,lp1+mp1+1)
2873                      do mm3= -l2,l2
2874                        c31=coeff(l2+1)%value(1,l2+mm3+1,l2+m2+1)
2875                        c32=coeff(l2+1)%value(2,l2+mm3+1,l2+m2+1)
2876                        fact=c11*c21*c31  -  c12*c22*c31&
2877 &                       -c11*c22*c32  -  c12*c21*c32
2878                        if((abs(fact)>=tol12).and.(mm3==-mm2-mm1)) &
2879 &                       realgnt_tmp=realgnt_tmp+fact*(-1)**mm2 &
2880 &                       *gaunt(l2,mm3,l1,mm1,lp1,-mm2)
2881                      end do
2882                    end do
2883                  end do
2884 
2885 !                Count additional non-zero real Gaunt coeffs
2886 !                ------------------------------------------------
2887                  if (abs(realgnt_tmp)>=tol12) then
2888                    ngnt=ngnt+1
2889                    gntselect(ilm2,klm1)=ngnt
2890                    realgnt(ngnt)=realgnt_tmp/sqrt(four_pi)
2891                  end if
2892 
2893 !                End loops
2894 !                ------------------------------------------------
2895                end if
2896              end do
2897            end do
2898          end if
2899        end do
2900      end do
2901    end do
2902  end do
2903 
2904 !Deallocate memory
2905 !------------------------------------------------
2906  do ll=1,4*l_max-3
2907    LIBPAW_DEALLOCATE(coeff(ll)%value)
2908  end do
2909  LIBPAW_DATATYPE_DEALLOCATE(coeff)
2910  LIBPAW_DEALLOCATE(ssgn)
2911 
2912 end subroutine realgaunt