TABLE OF CONTENTS
- ABINIT/m_paw_sphharm
- m_paw_sphharm/ass_leg_pol
- m_paw_sphharm/create_mlms2jmj
- m_paw_sphharm/create_slm2ylm
- m_paw_sphharm/dbeta
- m_paw_sphharm/gauleg
- m_paw_sphharm/gaunt
- m_paw_sphharm/initylmr
- m_paw_sphharm/lsylm
- m_paw_sphharm/lxyz.F90
- m_paw_sphharm/mat_mlms2jmj
- m_paw_sphharm/mat_slm2ylm
- m_paw_sphharm/mkeuler
- m_paw_sphharm/nablarealgaunt
- m_paw_sphharm/perms
- m_paw_sphharm/phim
- m_paw_sphharm/pl_deriv
- m_paw_sphharm/plm_coeff
- m_paw_sphharm/plm_d2theta
- m_paw_sphharm/plm_dphi
- m_paw_sphharm/plm_dtheta
- m_paw_sphharm/rfactorial
- m_paw_sphharm/setnabla_ylm
- m_paw_sphharm/setsym_ylm
- m_paw_sphharm/slxyzs
- m_paw_sphharm/ylm_angular_mesh
- m_paw_sphharm/ylm_cmplx
- m_paw_sphharm/ylmc
- m_paw_sphharm/ylmcd
- m_paw_sphharm/ys
- m_pawsphharm/realgaunt
ABINIT/m_paw_sphharm [ 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-2024 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 ]
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