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/dble_factorial
- m_paw_sphharm/initylmr
- m_paw_sphharm/lxyz.F90
- m_paw_sphharm/mat_mlms2jmj
- m_paw_sphharm/mat_slm2ylm
- m_paw_sphharm/mkeuler
- 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/setnabla_ylm
- m_paw_sphharm/setsym_ylm
- m_paw_sphharm/slxyzs
- m_paw_sphharm/ylm_cmplx
- m_paw_sphharm/ylmc
- m_paw_sphharm/ylmcd
- m_paw_sphharm/ys
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-2018 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
22 #include "libpaw.h" 23 24 MODULE m_paw_sphharm 25 26 USE_DEFS 27 USE_MSG_HANDLING 28 USE_MEMORY_PROFILING 29 30 implicit none 31 32 private 33 34 !public procedures. 35 public :: ylmc ! Complex Spherical harmonics for l<=3. 36 public :: ylmcd ! First derivative of complex Ylm wrt theta and phi up to l<=3 37 public :: ylm_cmplx ! All (complex) spherical harmonics for lx<=4 38 public :: initylmr ! Real Spherical Harmonics on a set of vectors 39 public :: ys ! Matrix element <Yl'm'|Slm> 40 public :: lxyz ! Matrix element <Yl'm'|L_idir|Ylm> 41 public :: slxyzs ! Matrix element <Sl'm'|L_idir|Slm> 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_d2theta ! d2(Plm (cos(theta)))/d(theta)2 (P_lm= ass. legendre polynomial) 45 public :: plm_dphi ! m*P_lm(x)/sqrt((1-x^2) (P_lm= ass. legendre polynomial) 46 public :: plm_dtheta ! -(1-x^2)^1/2*d/dx{P_lm(x)} (P_lm= ass. legendre polynomial) 47 public :: pl_deriv ! d2(Pl (x)))/d(x)2 where P_l is a legendre polynomial 48 public :: mkeuler ! For a given symmetry operation, determines the corresponding Euler angles 49 public :: dble_factorial ! Compute factorial of an integer; returns a double precision real 50 public :: dbeta ! Calculate the rotation matrix d^l_{m{\prim}m}(beta) 51 public :: phim ! Computes Phi_m[theta]=Sqrt[2] cos[m theta], if m>0 52 ! Sqrt[2] sin[Abs(m) theta], if m<0 53 ! 1 , if m=0 54 public :: mat_mlms2jmj ! Change a matrix from the Ylm basis to the J,M_J basis 55 public :: mat_slm2ylm ! Change a matrix from the Slm to the Ylm basis or from Ylm to Slm 56 public :: create_slm2ylm ! For a given angular momentum lcor, compute slm2ylm 57 public :: create_mlms2jmj ! For a given angular momentum lcor, give the rotation matrix msml2jmj 58 public :: setsym_ylm ! Compute rotation matrices expressed in the basis of real spherical harmonics 59 public :: setnabla_ylm ! Compute rotation matrices expressed in the basis of real spherical harmonics
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
PARENTS
CHILDREN
SOURCE
1116 function ass_leg_pol(l,m,xarg) 1117 1118 1119 !This section has been created automatically by the script Abilint (TD). 1120 !Do not modify the following lines by hand. 1121 #undef ABI_FUNC 1122 #define ABI_FUNC 'ass_leg_pol' 1123 !End of the abilint section 1124 1125 implicit none 1126 1127 !Arguments ------------------------------------ 1128 !scalars 1129 integer, intent(in) :: l,m 1130 real(dp), intent(in) :: xarg 1131 real(dp) :: ass_leg_pol 1132 1133 !Local variables------------------------------- 1134 !scalars 1135 integer :: i,ll 1136 real(dp) :: pll,polmm,tmp1,sqrx,x 1137 character(len=100) :: msg 1138 1139 ! ************************************************************************* 1140 1141 x=xarg 1142 if (m.lt.0.or.m.gt.l.or.abs(x).gt.1.d0) then 1143 if (m.lt.0.or.m.gt.l.or.abs(x).gt.1.d0+1.d-10) then 1144 msg='Bad choice of l, m or x !' 1145 MSG_BUG(msg) 1146 endif 1147 x=1.d0 1148 endif 1149 1150 polmm=1.d0 1151 if (m>0) then 1152 sqrx=sqrt(abs((1.d0-x)*(1.d0+x))) 1153 do i=1,m 1154 polmm=polmm*(1.0d0-2.0d0*i)*sqrx 1155 enddo 1156 endif 1157 1158 if (l==m) then 1159 ass_leg_pol=polmm 1160 else 1161 tmp1=x*(2.0d0*m+1.0d0)*polmm 1162 if (l==(m+1)) then 1163 ass_leg_pol=tmp1 1164 else 1165 do ll=m+2,l 1166 pll=(x*(2.0d0*ll-1.0d0)*tmp1-(ll+m-1.0d0)*polmm)/dble(ll-m) 1167 polmm=tmp1 1168 tmp1=pll 1169 enddo 1170 ass_leg_pol=pll 1171 endif 1172 endif 1173 1174 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
PARENTS
m_paw_sphharm
CHILDREN
SOURCE
2401 subroutine create_mlms2jmj(lcor,mlmstwojmj) 2402 2403 2404 !This section has been created automatically by the script Abilint (TD). 2405 !Do not modify the following lines by hand. 2406 #undef ABI_FUNC 2407 #define ABI_FUNC 'create_mlms2jmj' 2408 !End of the abilint section 2409 2410 implicit none 2411 2412 !Arguments --------------------------------------------- 2413 !scalars 2414 integer,intent(in) :: lcor 2415 !arrays 2416 complex(dpc),intent(out) :: mlmstwojmj(2*(2*lcor+1),2*(2*lcor+1)) 2417 2418 !Local variables --------------------------------------- 2419 !scalars 2420 integer :: jc1,jj,jm,ll,ml1,ms1 2421 real(dp) :: invsqrt2lp1,xj,xmj 2422 character(len=500) :: msg 2423 !arrays 2424 integer, allocatable :: ind_msml(:,:) 2425 complex(dpc),allocatable :: mat_mlms2(:,:) 2426 2427 !********************************************************************* 2428 2429 !--------------- Built indices + allocations 2430 ll=lcor 2431 mlmstwojmj=czero 2432 LIBPAW_ALLOCATE(ind_msml,(2,-ll:ll)) 2433 LIBPAW_ALLOCATE(mat_mlms2,(2*(2*lcor+1),2*(2*lcor+1))) 2434 mlmstwojmj=czero 2435 jc1=0 2436 do ms1=1,2 2437 do ml1=-ll,ll 2438 jc1=jc1+1 2439 ind_msml(ms1,ml1)=jc1 2440 end do 2441 end do 2442 2443 !--------------- built mlmstwojmj 2444 !do jj=ll,ll+1 ! the physical value of j are ll-0.5,ll+0.5 2445 !xj(jj)=jj-0.5 2446 if(ll==0)then 2447 msg=' ll should not be equal to zero !' 2448 MSG_BUG(msg) 2449 end if 2450 jc1=0 2451 invsqrt2lp1=one/sqrt(float(2*lcor+1)) 2452 do jj=ll,ll+1 2453 xj=float(jj)-half 2454 do jm=-jj,jj-1 2455 xmj=float(jm)+half 2456 jc1=jc1+1 2457 if(nint(xj+0.5)==ll+1) then 2458 if(nint(xmj+0.5)==ll+1) then 2459 mlmstwojmj(ind_msml(2,ll),jc1)=1.0 ! J=L+0.5 and m_J=L+0.5 2460 else if(nint(xmj-0.5)==-ll-1) then 2461 mlmstwojmj(ind_msml(1,-ll),jc1)=1.0 ! J=L+0.5 and m_J=-L-0.5 2462 else 2463 mlmstwojmj(ind_msml(2,nint(xmj-0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5)) 2464 mlmstwojmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)-xmj+0.5)) 2465 end if 2466 end if 2467 if(nint(xj+0.5)==ll) then 2468 mlmstwojmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5)) 2469 mlmstwojmj(ind_msml(2,nint(xmj-0.5)),jc1)=-invsqrt2lp1*(sqrt(float(ll)-xmj+0.5)) 2470 end if 2471 end do 2472 end do 2473 2474 LIBPAW_DEALLOCATE(ind_msml) 2475 LIBPAW_DEALLOCATE(mat_mlms2) 2476 2477 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
PARENTS
m_paw_sphharm
CHILDREN
SOURCE
2332 subroutine create_slm2ylm(lcor,slmtwoylm) 2333 2334 2335 !This section has been created automatically by the script Abilint (TD). 2336 !Do not modify the following lines by hand. 2337 #undef ABI_FUNC 2338 #define ABI_FUNC 'create_slm2ylm' 2339 !End of the abilint section 2340 2341 implicit none 2342 2343 !Arguments --------------------------------------------- 2344 !scalars 2345 integer,intent(in) :: lcor 2346 !arrays 2347 complex(dpc),intent(out) :: slmtwoylm(2*lcor+1,2*lcor+1) 2348 2349 !Local variables --------------------------------------- 2350 !scalars 2351 integer :: jm,ll,mm,im 2352 real(dp),parameter :: invsqrt2=one/sqrt2 2353 real(dp) :: onem 2354 !arrays 2355 2356 ! ********************************************************************* 2357 2358 ll=lcor 2359 slmtwoylm=czero 2360 do im=1,2*ll+1 2361 mm=im-ll-1;jm=-mm+ll+1 2362 onem=dble((-1)**mm) 2363 if (mm> 0) then 2364 slmtwoylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp) 2365 slmtwoylm(jm,im)= cmplx(invsqrt2, zero,kind=dp) 2366 end if 2367 if (mm==0) then 2368 slmtwoylm(im,im)=cone 2369 end if 2370 if (mm< 0) then 2371 slmtwoylm(im,im)= cmplx(zero, invsqrt2,kind=dp) 2372 slmtwoylm(jm,im)=-cmplx(zero,onem*invsqrt2,kind=dp) 2373 end if 2374 end do 2375 2376 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) 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
PARENTS
m_paw_sphharm
CHILDREN
SOURCE
1729 function dbeta(cosbeta,ll,mp,mm) 1730 1731 1732 !This section has been created automatically by the script Abilint (TD). 1733 !Do not modify the following lines by hand. 1734 #undef ABI_FUNC 1735 #define ABI_FUNC 'dbeta' 1736 !End of the abilint section 1737 1738 implicit none 1739 1740 !Arguments --------------------------------------------- 1741 !scalars 1742 integer,intent(in) :: ll,mm,mp 1743 real(dp) :: dbeta 1744 real(dp),intent(in) :: cosbeta 1745 1746 !Local variables ------------------------------ 1747 !scalars 1748 integer,parameter :: mxterms=200 1749 integer :: ii,ina,inb,inc,ml,ms 1750 real(dp) :: arg,cosbetab2,pref,sinbetab2,sum,tt 1751 1752 !************************************************************************ 1753 dbeta=zero 1754 1755 !Special cases 1756 if (abs(cosbeta-1._dp).lt.tol10) then 1757 if (mp.eq.mm) dbeta=1 1758 else if (abs(cosbeta+1._dp).lt.tol10) then 1759 if (mp.eq.-mm) dbeta=(-1)**(ll+mm) 1760 else 1761 1762 ! General case 1763 cosbetab2=sqrt((1+cosbeta)*0.5_dp) 1764 sinbetab2=sqrt((1-cosbeta)*0.5_dp) 1765 ml=max(mp,mm) 1766 ms=min(mp,mm) 1767 if (ml.ne.mp) sinbetab2=-sinbetab2 1768 tt=-(sinbetab2/cosbetab2)**2 1769 pref=sqrt((dble_factorial(ll-ms)*dble_factorial(ll+ml))& 1770 & /(dble_factorial(ll+ms)*dble_factorial(ll-ml)))& 1771 & /dble_factorial(ml-ms)*(cosbetab2**(2*ll+ms-ml))& 1772 & *((-sinbetab2)**(ml-ms)) 1773 sum=1._dp 1774 arg=1._dp 1775 ina=ml-ll 1776 inb=-ms-ll 1777 inc=ml-ms+1 1778 do ii=1,mxterms 1779 if (ina.eq.0.or.inb.eq.0) exit 1780 arg=(arg*ina*inb*tt)/(ii*inc) 1781 sum=sum+arg 1782 ina=ina+1 1783 inb=inb+1 1784 inc=inc+1 1785 end do 1786 dbeta=pref*sum 1787 end if 1788 1789 end function dbeta
m_paw_sphharm/dble_factorial [ Functions ]
[ Top ] [ m_paw_sphharm ] [ Functions ]
NAME
dble_factorial
FUNCTION
PRIVATE function Calculates N! as a double precision real.
INPUTS
nn=input integer
OUTPUT
factorial= N! (double precision)
PARENTS
CHILDREN
SOURCE
1666 elemental function dble_factorial(nn) 1667 1668 1669 !This section has been created automatically by the script Abilint (TD). 1670 !Do not modify the following lines by hand. 1671 #undef ABI_FUNC 1672 #define ABI_FUNC 'dble_factorial' 1673 !End of the abilint section 1674 1675 implicit none 1676 1677 !Arguments --------------------------------------------- 1678 !scalars 1679 integer,intent(in) :: nn 1680 real(dp) :: dble_factorial 1681 1682 !Local variables --------------------------------------- 1683 !scalars 1684 integer :: ii 1685 1686 ! ********************************************************************* 1687 1688 dble_factorial=one 1689 do ii=2,nn 1690 dble_factorial=dble_factorial*ii 1691 end do 1692 1693 end function dble_factorial
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}
PARENTS
debug_tools,denfgr,m_paw_finegrid,m_paw_pwaves_lmn,m_pawang mlwfovlp_ylmfar,posdoppler,pspnl_operat_rec,qijb_kk,smatrix_pawinit
CHILDREN
wrtout
SOURCE
573 subroutine initylmr(mpsang,normchoice,npts,nrm,option,rr,ylmr,ylmr_gr) 574 575 576 !This section has been created automatically by the script Abilint (TD). 577 !Do not modify the following lines by hand. 578 #undef ABI_FUNC 579 #define ABI_FUNC 'initylmr' 580 !End of the abilint section 581 582 implicit none 583 584 !Arguments ------------------------------------ 585 !scalars 586 integer,intent(in) :: mpsang,normchoice,npts,option 587 !arrays 588 real(dp),intent(in) :: nrm(npts),rr(3,npts) 589 real(dp),intent(out) :: ylmr(mpsang*mpsang,npts) 590 real(dp),optional,intent(out) :: ylmr_gr(3*(option/2)+6*(option/3),mpsang*mpsang,npts) 591 592 !Local variables ------------------------------ 593 !scalars 594 integer :: dimgr,ilang,inpt,l0,ll,mm 595 real(dp) :: cphi,ctheta,fact,onem,rnorm,sphi,stheta,work1,work2,ylmcst,ylmcst2 596 logical :: compute_ylm,compute_ylm2gr,compute_ylmgr 597 !arrays 598 real(dp) :: dphi(3),dtheta(3),iphase(mpsang-1),rphase(mpsang-1) 599 real(dp),allocatable :: blm(:,:) 600 601 !************************************************************************ 602 603 !What has to be computed ? 604 compute_ylm = (option==1.or.option==2.or.option==3) 605 compute_ylmgr =(( option==2.or.option==3).and.present(ylmr_gr)) 606 compute_ylm2gr=(( option==3).and.present(ylmr_gr)) 607 dimgr=3*(option/2)+6*(option/3) 608 609 !Initialisation of spherical harmonics 610 if (compute_ylm ) ylmr (: ,1:npts)=zero 611 if (compute_ylmgr) ylmr_gr(:,:,1:npts)=zero 612 613 !Special case for l=0 614 if (compute_ylm ) ylmr(1,1:npts)=1._dp/sqrt(four_pi) 615 if (compute_ylmgr) ylmr_gr(1:dimgr,1,1:npts)=zero 616 if (mpsang>1) then 617 618 ! Loop over all rr 619 do inpt=1,npts 620 621 ! Load module of rr 622 rnorm=one 623 if (normchoice==1) rnorm=nrm(inpt) 624 625 ! Continue only for r<>0 626 627 if (rnorm>tol10) then 628 629 ! Determine theta and phi 630 cphi=one 631 sphi=zero 632 ctheta=rr(3,inpt)/rnorm 633 ! MM030519 : abs is needed to prevent very small negative arg 634 stheta=sqrt(abs((one-ctheta)*(one+ctheta))) 635 if (stheta>tol10) then 636 cphi=rr(1,inpt)/(rnorm*stheta) 637 sphi=rr(2,inpt)/(rnorm*stheta) 638 end if 639 do mm=1,mpsang-1 640 rphase(mm)=dreal(dcmplx(cphi,sphi)**mm) 641 iphase(mm)=aimag(dcmplx(cphi,sphi)**mm) 642 end do 643 644 ! Determine gradients of theta and phi 645 if (compute_ylmgr) then 646 dtheta(1)=ctheta*cphi 647 dtheta(2)=ctheta*sphi 648 dtheta(3)=-stheta 649 dphi(1)=-sphi 650 dphi(2)=cphi 651 dphi(3)=zero 652 end if 653 654 ! COMPUTE Ylm(R) 655 if (compute_ylm) then 656 ! Loop over angular momentum l 657 do ilang=2,mpsang 658 ll=ilang-1 659 l0=ll**2+ll+1 660 fact=1._dp/real(ll*(ll+1),dp) 661 ylmcst=sqrt(real(2*ll+1,dp)/four_pi) 662 ! Special case m=0 663 ylmr(l0,inpt)=ylmcst*ass_leg_pol(ll,0,ctheta) 664 ! Compute for m>0 665 onem=one 666 do mm=1,ll 667 onem=-onem 668 work1=ylmcst*sqrt(fact)*onem*ass_leg_pol(ll,mm,ctheta)*sqrt(2._dp) 669 ylmr(l0+mm,inpt)=work1*rphase(mm) 670 ylmr(l0-mm,inpt)=work1*iphase(mm) 671 if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp) 672 end do ! End loop over m 673 end do ! End loop over l 674 end if 675 676 ! COMPUTE dYlm/dRi 677 if (compute_ylmgr) then 678 ! Loop over angular momentum l 679 do ilang=2,mpsang 680 ll=ilang-1 681 l0=ll**2+ll+1 682 fact=1._dp/real(ll*(ll+1),dp) 683 ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/rnorm 684 ! Special case m=0 685 work1=ylmcst*plm_dtheta(ll,0,ctheta) 686 ylmr_gr(1:3,l0,inpt)=work1*dtheta(1:3) 687 ! Compute for m>0 688 onem=one 689 do mm=1,ll 690 onem=-onem 691 work1=ylmcst*sqrt(fact)*onem*plm_dtheta(ll,mm,ctheta)*sqrt(2._dp) 692 work2=ylmcst*sqrt(fact)*onem*plm_dphi (ll,mm,ctheta)*sqrt(2._dp) 693 ylmr_gr(1:3,l0+mm,inpt)=rphase(mm)*work1*dtheta(1:3)-iphase(mm)*work2*dphi(1:3) 694 ylmr_gr(1:3,l0-mm,inpt)=iphase(mm)*work1*dtheta(1:3)+rphase(mm)*work2*dphi(1:3) 695 if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp) 696 end do ! End loop over m 697 end do ! End loop over l 698 end if 699 700 ! COMPUTE d2Ylm/dRidRj 701 if (compute_ylm2gr) then 702 LIBPAW_ALLOCATE(blm,(5,mpsang*mpsang)) 703 call plm_coeff(blm,mpsang,ctheta) 704 705 ! Loop over angular momentum l 706 do ilang=2,mpsang 707 ll=ilang-1 708 l0=ll**2+ll+1 709 fact=1._dp/real(ll*(ll+1),dp) 710 ylmcst=sqrt(real(2*ll+1,dp)/four_pi)/(rnorm**2) 711 ! Special case m=0 712 ylmr_gr(4,l0,inpt)=ylmcst*(-blm(3,l0)*sphi*sphi+blm(4,l0)*cphi*cphi) 713 ylmr_gr(5,l0,inpt)=ylmcst*(-blm(3,l0)*cphi*cphi+blm(4,l0)*sphi*sphi) 714 ylmr_gr(6,l0,inpt)=ylmcst*blm(1,l0) 715 ylmr_gr(7,l0,inpt)=ylmcst*blm(2,l0)*sphi 716 ylmr_gr(8,l0,inpt)=ylmcst*blm(2,l0)*cphi 717 ylmr_gr(9,l0,inpt)=ylmcst*(blm(3,l0)+blm(4,l0))*sphi*cphi 718 ! Compute for m>0 719 onem=one 720 do mm=1,ll 721 onem=-onem;ylmcst2=ylmcst*sqrt(fact)*sqrt(two) 722 ylmr_gr(4,l0+mm,inpt)=ylmcst2*((-blm(3,l0+mm)*sphi*sphi+blm(4,l0+mm)*cphi*cphi)*rphase(mm)-& 723 & blm(5,l0+mm)*2.d0*cphi*sphi*mm*iphase(mm)) 724 ylmr_gr(4,l0-mm,inpt)=ylmcst2*((-blm(3,l0+mm)*sphi*sphi+blm(4,l0+mm)*cphi*cphi)*iphase(mm)+& 725 & blm(5,l0+mm)*2.d0*cphi*sphi*mm*rphase(mm)) 726 ylmr_gr(5,l0+mm,inpt)=ylmcst2*((-blm(3,l0+mm)*cphi*cphi+blm(4,l0+mm)*sphi*sphi)*rphase(mm)+& 727 & blm(5,l0+mm)*2.d0*cphi*sphi*mm*iphase(mm)) 728 ylmr_gr(5,l0-mm,inpt)=ylmcst2*((-blm(3,l0+mm)*cphi*cphi+blm(4,l0+mm)*sphi*sphi)*iphase(mm)-& 729 & blm(5,l0+mm)*2.d0*cphi*sphi*mm*rphase(mm)) 730 ylmr_gr(6,l0+mm,inpt)=ylmcst2*blm(1,l0+mm)*rphase(mm) 731 ylmr_gr(6,l0-mm,inpt)=ylmcst2*blm(1,l0+mm)*iphase(mm) 732 ylmr_gr(7,l0+mm,inpt)=ylmcst2*(blm(2,l0+mm)*sphi*rphase(mm)+& 733 & mm*iphase(mm)*cphi*onem*plm_dtheta(ll,mm,ctheta)) 734 ylmr_gr(7,l0-mm,inpt)=ylmcst2*(blm(2,l0+mm)*sphi*iphase(mm)-& 735 & mm*rphase(mm)*cphi*onem*plm_dtheta(ll,mm,ctheta)) 736 ylmr_gr(8,l0+mm,inpt)=ylmcst2*(blm(2,l0+mm)*cphi*rphase(mm)-& 737 & mm*iphase(mm)*sphi*onem*plm_dtheta(ll,mm,ctheta)) 738 ylmr_gr(8,l0-mm,inpt)=ylmcst2*(blm(2,l0+mm)*cphi*iphase(mm)+& 739 & mm*rphase(mm)*sphi*onem*plm_dtheta(ll,mm,ctheta)) 740 ylmr_gr(9,l0+mm,inpt)=ylmcst2*((blm(3,l0+mm)+blm(4,l0+mm))*sphi*cphi*rphase(mm)-& 741 & blm(5,l0+mm)*(sphi*sphi-cphi*cphi)*mm*iphase(mm)) 742 ylmr_gr(9,l0-mm,inpt)=ylmcst2*((blm(3,l0+mm)+blm(4,l0+mm))*sphi*cphi*iphase(mm)+& 743 & blm(5,l0+mm)*(sphi*sphi-cphi*cphi)*mm*rphase(mm)) 744 if (mm/=ll) fact=fact/real((ll+mm+1)*(ll-mm),dp) 745 end do ! End loop over m 746 end do ! End loop over l 747 LIBPAW_DEALLOCATE(blm) 748 end if 749 750 ! End condition r<>0 751 end if 752 753 ! End loop over rr 754 end do 755 756 ! End condition l<>0 757 end if 758 759 end subroutine initylmr
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
PARENTS
m_paw_sphharm
CHILDREN
wrtout
SOURCE
865 subroutine lxyz(lp,mp,idir,ll,mm,lidir) 866 867 868 !This section has been created automatically by the script Abilint (TD). 869 !Do not modify the following lines by hand. 870 #undef ABI_FUNC 871 #define ABI_FUNC 'lxyz' 872 !End of the abilint section 873 874 implicit none 875 876 !Arguments --------------------------------------------- 877 !scalars 878 integer,intent(in) :: idir,ll,lp,mm,mp 879 complex(dpc),intent(out) :: lidir 880 881 !Local variables --------------------------------------- 882 !scalars 883 complex(dpc) :: jmme, jpme 884 885 ! ********************************************************************* 886 887 jpme=czero; jmme=czero 888 if (lp==ll) then 889 if (mp==mm+1) jpme=cone*sqrt((ll-mm)*(ll+mm+one)) 890 if (mp==mm-1) jmme=cone*sqrt((ll-mm+one)*(ll+mm)) 891 end if 892 893 lidir = czero 894 if (lp == ll) then 895 select case (idir) 896 case (1) ! Lx 897 lidir = cone*half*(jpme+jmme) 898 case (2) ! Ly 899 lidir = -(zero,one)*half*(jpme-jmme) 900 case (3) ! Lz 901 if (mp == mm) lidir = mm*cone 902 end select 903 end if 904 905 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
PARENTS
m_pawang,pawprt,setnoccmmp
CHILDREN
wrtout
SOURCE
1892 subroutine mat_mlms2jmj(lcor,mat_mlms,mat_jmj,ndij,option,optspin,prtvol,unitfi,wrt_mode) 1893 1894 1895 !This section has been created automatically by the script Abilint (TD). 1896 !Do not modify the following lines by hand. 1897 #undef ABI_FUNC 1898 #define ABI_FUNC 'mat_mlms2jmj' 1899 !End of the abilint section 1900 1901 implicit none 1902 1903 !Arguments --------------------------------------------- 1904 !scalars 1905 integer,intent(in) :: ndij,lcor,option,optspin,prtvol,unitfi 1906 character(len=4),intent(in) :: wrt_mode 1907 !arrays 1908 complex(dpc),intent(inout) :: mat_mlms(2*lcor+1,2*lcor+1,ndij) 1909 complex(dpc),intent(inout) :: mat_jmj(2*(2*lcor+1),2*(2*lcor+1)) 1910 1911 !Local variables --------------------------------------- 1912 !scalars 1913 integer :: ii,im,im1,im2,ispden,jc1,jc2,jj,jm,ll,ml1,ml2,ms1,ms2 1914 real(dp),parameter :: invsqrt2=one/sqrt2 1915 real(dp) :: invsqrt2lp1,xj,xmj 1916 complex(dpc) :: mat_tmp,tmp2 1917 character(len=9),parameter :: dspinold(6)=(/"up ","down ","up-up ","down-down","up-dn ","dn-up "/) 1918 character(len=9),parameter :: dspin(6)=(/"dn ","up ","dn-dn ","up-up ","dn-up ","up-dn "/) 1919 character(len=500) :: msg 1920 !arrays 1921 integer, allocatable :: ind_msml(:,:) 1922 complex(dpc),allocatable :: mat_mlms2(:,:),mlms2jmj(:,:) 1923 1924 !********************************************************************* 1925 1926 if(ndij/=4) then 1927 msg=" ndij/=4 !" 1928 MSG_BUG(msg) 1929 end if 1930 if (option/=1.and.option/=2) then 1931 msg=' option=/1 and =/2 !' 1932 MSG_BUG(msg) 1933 end if 1934 if (optspin/=1.and.optspin/=2) then 1935 msg=' optspin=/1 and =/2 !' 1936 MSG_BUG(msg) 1937 end if 1938 1939 if (unitfi/=-1) then 1940 if(option==1) then 1941 write(msg,'(3a)') ch10,& 1942 & "matrix in |l,s,m_l,m_s> basis is changed into |l,s,j,m_j> basis" 1943 call wrtout(unitfi,msg,wrt_mode) 1944 else if(option==2) then 1945 write(msg,'(3a)') ch10,& 1946 & "matrix in |l,s,j,m_j> basis is changed into |l,s,m_l,m_s> basis" 1947 call wrtout(unitfi,msg,wrt_mode) 1948 end if 1949 end if 1950 1951 if(option==1) then 1952 if(optspin==2) then 1953 if(abs(prtvol)>2.and.unitfi/=-1)& 1954 & write(msg,'(3a)') ch10,"assume spin dn is the first in the array" 1955 else if (optspin==1) then 1956 if(abs(prtvol)>2.and.unitfi/=-1)& 1957 & write(msg,'(3a)') ch10,"change array in order that spin dn is the first in the array" 1958 do ii=1,2*lcor+1 1959 do jj=1,2*lcor+1 1960 mat_tmp=mat_mlms(ii,jj,2) 1961 mat_mlms(ii,jj,2)=mat_mlms(ii,jj,1) 1962 mat_mlms(ii,jj,1)=mat_tmp 1963 mat_tmp=mat_mlms(ii,jj,4) 1964 mat_mlms(ii,jj,4)=mat_mlms(ii,jj,3) 1965 mat_mlms(ii,jj,3)=mat_tmp 1966 end do 1967 end do 1968 ! mat_tmp(:,:,1)=mat_mlms(:,:,2);mat_tmp(:,:,2)=mat_mlms(:,:,1) 1969 ! mat_tmp(:,:,3)=mat_mlms(:,:,4);mat_tmp(:,:,4)=mat_mlms(:,:,3) 1970 ! mat_mlms(:,:,:)=mat_tmp(:,:,:) 1971 end if 1972 if(abs(prtvol)>2.and.unitfi/=-1) then 1973 call wrtout(unitfi,msg,wrt_mode) 1974 end if 1975 end if 1976 1977 if(option==1.and.abs(prtvol)>2.and.unitfi/=-1) then 1978 do ispden=1,ndij 1979 write(msg,'(3a)') ch10,& 1980 & "Input matrix in the Ylm basis for component ",trim(dspin(ispden+2*(ndij/4))) 1981 call wrtout(unitfi,msg,wrt_mode) 1982 do im1=1,lcor*2+1 1983 write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))')& 1984 & (mat_mlms(im1,im2,ispden),im2=1,lcor*2+1) 1985 call wrtout(unitfi,msg,wrt_mode) 1986 end do 1987 end do 1988 end if ! option==1 1989 1990 !--------------- Built indices + allocations 1991 ll=lcor 1992 LIBPAW_ALLOCATE(mlms2jmj,(2*(2*ll+1),2*(2*ll+1))) 1993 mlms2jmj=czero 1994 LIBPAW_BOUND2_ALLOCATE(ind_msml,BOUNDS(1,2),BOUNDS(-ll,ll)) 1995 LIBPAW_ALLOCATE(mat_mlms2,(2*(2*lcor+1),2*(2*lcor+1))) 1996 mlms2jmj=czero 1997 jc1=0 1998 do ms1=1,2 1999 do ml1=-ll,ll 2000 jc1=jc1+1 2001 ind_msml(ms1,ml1)=jc1 2002 end do 2003 end do 2004 !--------------- Change representation of input matrix for ndij==4 2005 if(option==1) then 2006 jc1=0 2007 do ms1=1,2 2008 do ml1=1,2*ll+1 2009 jc1=jc1+1 2010 jc2=0 2011 do ms2=1,2 2012 do ml2=1,2*ll+1 2013 jc2=jc2+1 2014 if(ms1==ms2) mat_mlms2(jc1,jc2)=mat_mlms(ml1,ml2,ms1) 2015 if(ms1<ms2) mat_mlms2(jc1,jc2)=mat_mlms(ml1,ml2,3) 2016 if(ms1>ms2) mat_mlms2(jc1,jc2)=mat_mlms(ml1,ml2,4) 2017 end do 2018 end do 2019 end do 2020 end do 2021 if(abs(prtvol)>1.and.unitfi/=-1) then 2022 write(msg,'(3a)') ch10,"Input matrix in the lms basis for all component" 2023 call wrtout(unitfi,msg,wrt_mode) 2024 do im1=1,2*(lcor*2+1) 2025 write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))')& 2026 & (mat_mlms2(im1,im2),im2=1,2*(lcor*2+1)) 2027 call wrtout(unitfi,msg,wrt_mode) 2028 end do 2029 end if 2030 end if ! option==1 2031 2032 !--------------- built mlms2jmj 2033 !do jj=ll,ll+1 ! the physical value of j are ll-0.5,ll+0.5 2034 !xj(jj)=jj-0.5 2035 if(ll==0)then 2036 msg=' ll should not be equal to zero !' 2037 MSG_BUG(msg) 2038 end if 2039 jc1=0 2040 invsqrt2lp1=one/sqrt(float(2*lcor+1)) 2041 do jj=ll,ll+1 2042 xj=float(jj)-half 2043 do jm=-jj,jj-1 2044 xmj=float(jm)+half 2045 jc1=jc1+1 2046 if(nint(xj+0.5)==ll+1) then 2047 if(nint(xmj+0.5)==ll+1) then 2048 mlms2jmj(ind_msml(2,ll),jc1)=1.0 ! J=L+0.5 and m_J=L+0.5 2049 else if(nint(xmj-0.5)==-ll-1) then 2050 mlms2jmj(ind_msml(1,-ll),jc1)=1.0 ! J=L+0.5 and m_J=-L-0.5 2051 else 2052 mlms2jmj(ind_msml(2,nint(xmj-0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5)) 2053 mlms2jmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)-xmj+0.5)) 2054 end if 2055 end if 2056 if(nint(xj+0.5)==ll) then 2057 mlms2jmj(ind_msml(1,nint(xmj+0.5)),jc1)=invsqrt2lp1*(sqrt(float(ll)+xmj+0.5)) 2058 mlms2jmj(ind_msml(2,nint(xmj-0.5)),jc1)=-invsqrt2lp1*(sqrt(float(ll)-xmj+0.5)) 2059 end if 2060 end do 2061 end do 2062 if(abs(prtvol)>2.and.unitfi/=-1) then 2063 write(msg,'(3a)') ch10,"Matrix to go from |M_L,M_S> to |J,M_J>" 2064 call wrtout(unitfi,msg,wrt_mode) 2065 do im1=1,2*(lcor*2+1) 2066 write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') (mlms2jmj(im1,im2),im2=1,2*(lcor*2+1)) 2067 call wrtout(unitfi,msg,wrt_mode) 2068 end do 2069 end if 2070 2071 do jm=1,2*(2*ll+1) 2072 do im=1,2*(2*ll+1) 2073 tmp2=czero 2074 do ii=1,2*(2*ll+1) 2075 do jj=1,2*(2*ll+1) 2076 if(option==1) then 2077 tmp2=tmp2+mat_mlms2(ii,jj)*CONJG(mlms2jmj(ii,im))*(mlms2jmj(jj,jm)) 2078 else if(option==2) then 2079 tmp2=tmp2+mat_jmj(ii,jj)*(mlms2jmj(im,ii))*CONJG(mlms2jmj(jm,jj)) ! inv=t* 2080 end if 2081 end do 2082 end do 2083 if(option==1) then 2084 mat_jmj(im,jm)=tmp2 2085 else if(option==2) then 2086 mat_mlms2(im,jm)=tmp2 2087 end if 2088 end do 2089 end do 2090 if(option==1) then 2091 if (abs(prtvol)>=1.and.unitfi/=-1) then 2092 write(msg,'(3a)') ch10," Matrix in the J,M_J basis" 2093 call wrtout(unitfi,msg,wrt_mode) 2094 do im1=1,2*(lcor*2+1) 2095 write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') (mat_jmj(im1,im2),im2=1,2*(lcor*2+1)) 2096 call wrtout(unitfi,msg,wrt_mode) 2097 end do 2098 end if 2099 else if(option==2) then 2100 if (abs(prtvol)>=1.and.unitfi/=-1) then 2101 write(msg,'(3a)') ch10," Matrix in the m_s m_l basis" 2102 call wrtout(unitfi,msg,wrt_mode) 2103 do im1=1,2*(lcor*2+1) 2104 write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') (mat_mlms2(im1,im2),im2=1,2*(lcor*2+1)) 2105 call wrtout(unitfi,msg,wrt_mode) 2106 end do 2107 end if 2108 jc1=0 2109 do ms1=1,2 2110 do ml1=1,2*ll+1 2111 jc1=jc1+1 2112 jc2=0 2113 do ms2=1,2 2114 do ml2=1,2*ll+1 2115 jc2=jc2+1 2116 if(ms1==ms2) mat_mlms(ml1,ml2,ms1)=mat_mlms2(jc1,jc2) 2117 if(ms1<ms2) mat_mlms(ml1,ml2,3)=mat_mlms2(jc1,jc2) 2118 if(ms1>ms2) mat_mlms(ml1,ml2,4)=mat_mlms2(jc1,jc2) 2119 end do 2120 end do 2121 end do 2122 end do 2123 end if 2124 LIBPAW_DEALLOCATE(mlms2jmj) 2125 LIBPAW_DEALLOCATE(mat_mlms2) 2126 LIBPAW_DEALLOCATE(ind_msml) 2127 2128 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
PARENTS
m_pawang,pawprt,setnoccmmp
CHILDREN
wrtout
SOURCE
2167 subroutine mat_slm2ylm(lcor,mat_inp_c,mat_out_c,ndij,option,optspin,prtvol,unitfi,wrt_mode) 2168 2169 2170 !This section has been created automatically by the script Abilint (TD). 2171 !Do not modify the following lines by hand. 2172 #undef ABI_FUNC 2173 #define ABI_FUNC 'mat_slm2ylm' 2174 !End of the abilint section 2175 2176 implicit none 2177 2178 !Arguments --------------------------------------------- 2179 !scalars 2180 integer,intent(in) :: ndij,lcor,option,optspin,prtvol,unitfi 2181 character(len=4),intent(in) :: wrt_mode 2182 !arrays 2183 complex(dpc) :: mat_inp_c(2*lcor+1,2*lcor+1,ndij),mat_out(2*lcor+1,2*lcor+1,ndij) 2184 complex(dpc) :: mat_out_c(2*lcor+1,2*lcor+1,ndij) 2185 2186 !Local variables --------------------------------------- 2187 !scalars 2188 integer :: jm,ii,jj,ll,mm,ispden,im,im1,im2 2189 real(dp),parameter :: invsqrt2=one/sqrt2 2190 real(dp) :: onem 2191 complex(dpc) :: tmp2 2192 character(len=9),parameter :: dspinc(6)=(/"up ","down ","up-up ","down-down","up-dn ","dn-up "/)! optspin 1 2193 character(len=9),parameter :: dspinc2(6)=(/"up ","down ","dn-dn ","up-up ","dn-up ","up-dn "/)! optspin 2 2194 character(len=500) :: msg 2195 !arrays 2196 complex(dpc),allocatable :: slm2ylm(:,:) 2197 2198 ! ********************************************************************* 2199 2200 if(ndij/=4) then 2201 msg=' ndij:=4 !' 2202 MSG_BUG(msg) 2203 end if 2204 if (option/=1.and.option/=2.and.option/=3.and.option/=4) then 2205 msg=' option=/1 or 2 or 3 or 4 !' 2206 MSG_BUG(msg) 2207 end if 2208 2209 if(abs(prtvol)>2.and.unitfi/=-1) then 2210 write(msg,'(3a)') ch10, " mat_slm2ylm" 2211 call wrtout(unitfi,msg,wrt_mode) 2212 end if 2213 2214 if(abs(prtvol)>2.and.unitfi/=-1) then 2215 if(option==1.or.option==3) then 2216 write(msg,'(3a)') ch10,"matrix in Slm basis is changed into Ylm basis" 2217 call wrtout(unitfi,msg,wrt_mode) 2218 else if(option==2.or.option==4) then 2219 write(msg,'(3a)') ch10,"matrix in Ylm basis is changed into Slm basis" 2220 call wrtout(unitfi,msg,wrt_mode) 2221 end if 2222 end if 2223 2224 ll=lcor 2225 LIBPAW_ALLOCATE(slm2ylm,(2*ll+1,2*ll+1)) 2226 slm2ylm=czero 2227 mat_out=zero 2228 mat_out_c=czero 2229 do im=1,2*ll+1 2230 mm=im-ll-1;jm=-mm+ll+1 2231 onem=dble((-1)**mm) 2232 if (mm> 0) then 2233 slm2ylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp) 2234 slm2ylm(jm,im)= cmplx(invsqrt2, zero,kind=dp) 2235 end if 2236 if (mm==0) then 2237 slm2ylm(im,im)=cone 2238 end if 2239 if (mm< 0) then 2240 slm2ylm(im,im)= cmplx(zero, invsqrt2,kind=dp) 2241 slm2ylm(jm,im)=-cmplx(zero,onem*invsqrt2,kind=dp) 2242 end if 2243 end do 2244 if(abs(prtvol)>2.and.unitfi/=-1) then 2245 do ispden=1,ndij 2246 if(optspin==1) then 2247 if(option==1.or.option==3)& 2248 & write(msg,'(3a)') ch10,& 2249 & "Input matrix in the Slm basis for component ",trim(dspinc(ispden+2*(ndij/4))) 2250 if(option==2.or.option==3)& 2251 & write(msg,'(3a)') ch10,& 2252 & "Input matrix in the Ylm basis for component ",trim(dspinc(ispden+2*(ndij/4))) 2253 else 2254 if(option==1.or.option==3)& 2255 & write(msg,'(3a)') ch10,& 2256 & "Input matrix in the Slm basis for component ",trim(dspinc2(ispden+2*(ndij/4))) 2257 if(option==2.or.option==3)& 2258 & write(msg,'(3a)') ch10,& 2259 & "Input matrix in the Ylm basis for component ",trim(dspinc2(ispden+2*(ndij/4))) 2260 end if 2261 call wrtout(unitfi,msg,wrt_mode) 2262 do im1=1,lcor*2+1 2263 write(msg,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')& 2264 & (mat_inp_c(im1,im2,ispden),im2=1,lcor*2+1) 2265 call wrtout(unitfi,msg,wrt_mode) 2266 end do 2267 end do 2268 end if 2269 do ispden=1,ndij 2270 do jm=1,2*ll+1 2271 do im=1,2*ll+1 2272 tmp2=czero 2273 do ii=1,2*ll+1 2274 do jj=1,2*ll+1 2275 if(option==1) then 2276 tmp2=tmp2+mat_inp_c(ii,jj,ispden)*(slm2ylm(im,ii))*CONJG(slm2ylm(jm,jj)) 2277 else if(option==2) then 2278 tmp2=tmp2+mat_inp_c(ii,jj,ispden)*CONJG(slm2ylm(ii,im))*(slm2ylm(jj,jm)) 2279 end if 2280 end do 2281 end do 2282 mat_out_c(im,jm,ispden)=tmp2 2283 end do 2284 end do 2285 end do ! ispden 2286 do ii=1,2*ll+1 2287 do jj=1,2*ll+1 2288 mat_out(ii,jj,1)=real(mat_out_c(ii,jj,1)) 2289 mat_out(ii,jj,2)=real(mat_out_c(ii,jj,2)) 2290 mat_out(ii,jj,3)=real(mat_out_c(ii,jj,3)) 2291 mat_out(ii,jj,4)=aimag(mat_out_c(ii,jj,3)) 2292 ! check that n_{m,m'}^{alpha,beta}=conjg(n_{m',m"}^{beta,alpha}). 2293 if((abs(aimag(mat_out_c(ii,jj,3))+aimag(mat_out_c(jj,ii,4))).ge.0.0001).or. & 2294 & (abs(real(mat_out_c(ii,jj,3))-real(mat_out_c(jj,ii,4))).ge.0.0001)) then 2295 write(msg,'(a,4f10.4)') & 2296 & ' prb with mat_out_c ',mat_out_c(ii,jj,3),mat_out_c(ii,jj,4) 2297 MSG_BUG(msg) 2298 end if 2299 end do 2300 end do 2301 2302 LIBPAW_DEALLOCATE(slm2ylm) 2303 2304 end subroutine mat_slm2ylm
m_paw_sphharm/mkeuler [ Functions ]
[ Top ] [ m_paw_sphharm ] [ Functions ]
NAME
mkeuler
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 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
PARENTS
m_paw_sphharm
CHILDREN
SOURCE
1567 subroutine mkeuler(rot,cosbeta,cosalp,sinalp,cosgam,singam,isn) 1568 1569 1570 !This section has been created automatically by the script Abilint (TD). 1571 !Do not modify the following lines by hand. 1572 #undef ABI_FUNC 1573 #define ABI_FUNC 'mkeuler' 1574 !End of the abilint section 1575 1576 implicit none 1577 1578 !Arguments --------------------------------------------- 1579 !scalars 1580 integer,intent(out) :: isn 1581 real(dp),intent(out) :: cosalp,cosbeta,cosgam,sinalp,singam 1582 !arrays 1583 real(dp),intent(in) :: rot(3,3) 1584 1585 !Local variables --------------------------------------- 1586 !scalars 1587 integer :: ier 1588 real(dp) :: check,sinbeta 1589 character(len=500) :: msg 1590 1591 ! ********************************************************************* 1592 1593 do isn= -1,1,2 1594 cosbeta=real(isn)*rot(3,3) 1595 if(abs(1._dp-cosbeta*cosbeta)<tol10) then 1596 sinbeta=zero 1597 else 1598 sinbeta=sqrt(1._dp-cosbeta*cosbeta) 1599 end if 1600 if (abs(sinbeta).gt.tol10) then 1601 cosalp=isn*rot(3,1)/sinbeta 1602 sinalp=isn*rot(3,2)/sinbeta 1603 cosgam=-isn*rot(1,3)/sinbeta 1604 singam=isn*rot(2,3)/sinbeta 1605 else 1606 cosalp=isn*rot(1,1)/cosbeta 1607 sinalp=isn*rot(1,2)/cosbeta 1608 cosgam=one 1609 singam=zero 1610 end if 1611 1612 ! Check matrix: 1613 ier=0 1614 check=cosalp*cosbeta*cosgam-sinalp*singam 1615 if (abs(check-isn*rot(1,1))>tol8) ier=ier+1 1616 check=sinalp*cosbeta*cosgam+cosalp*singam 1617 if (abs(check-isn*rot(1,2))>tol8) ier=ier+1 1618 check=-sinbeta*cosgam 1619 if (abs(check-isn*rot(1,3))>tol8) ier=ier+1 1620 check=-cosalp*cosbeta*singam-sinalp*cosgam 1621 if (abs(check-isn*rot(2,1))>tol8) ier=ier+1 1622 check=-sinalp*cosbeta*singam+cosalp*cosgam 1623 if (abs(check-isn*rot(2,2))>tol8) ier=ier+1 1624 check=sinbeta*singam 1625 if (abs(check-isn*rot(2,3))>tol8) ier=ier+1 1626 check=cosalp*sinbeta 1627 if (abs(check-isn*rot(3,1))>tol8) ier=ier+1 1628 check=sinalp*sinbeta 1629 if (abs(check-isn*rot(3,2))>tol8) ier=ier+1 1630 if (ier.eq.0) return 1631 end do 1632 1633 isn=0 1634 write(msg, '(7a)' )& 1635 & 'Error during determination of symetries!',ch10,& 1636 & 'Action: check your input file:',ch10,& 1637 & 'unit cell vectors and/or atoms positions',ch10,& 1638 & 'have to be given with a better precision.' 1639 MSG_ERROR(msg) 1640 1641 end subroutine mkeuler
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
PARENTS
m_paw_sphharm
CHILDREN
SOURCE
1822 pure function phim(costheta,sintheta,mm) 1823 1824 1825 !This section has been created automatically by the script Abilint (TD). 1826 !Do not modify the following lines by hand. 1827 #undef ABI_FUNC 1828 #define ABI_FUNC 'phim' 1829 !End of the abilint section 1830 1831 implicit none 1832 1833 !Arguments --------------------------------------------- 1834 !scalars 1835 integer,intent(in) :: mm 1836 real(dp) :: phim 1837 real(dp),intent(in) :: costheta,sintheta 1838 1839 ! ********************************************************************* 1840 1841 if (mm==0) phim=one 1842 if (mm==1) phim=sqrt2*costheta 1843 if (mm==-1) phim=sqrt2*sintheta 1844 if (mm==2) phim=sqrt2*(costheta*costheta-sintheta*sintheta) 1845 if (mm==-2) phim=sqrt2*two*sintheta*costheta 1846 if (mm==3) phim=sqrt2*& 1847 & (costheta*(costheta*costheta-sintheta*sintheta)& 1848 & -sintheta*two*sintheta*costheta) 1849 if (mm==-3) phim=sqrt2*& 1850 & (sintheta*(costheta*costheta-sintheta*sintheta)& 1851 & +costheta*two*sintheta*costheta) 1852 1853 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)
PARENTS
m_paw_sphharm
CHILDREN
wrtout
SOURCE
1486 subroutine pl_deriv(mpsang,pl_d2,xx) 1487 1488 1489 !This section has been created automatically by the script Abilint (TD). 1490 !Do not modify the following lines by hand. 1491 #undef ABI_FUNC 1492 #define ABI_FUNC 'pl_deriv' 1493 !End of the abilint section 1494 1495 implicit none 1496 1497 !Arguments --------------------------------------------- 1498 !scalars 1499 integer,intent(in) :: mpsang 1500 real(dp),intent(in) :: xx 1501 !arrays 1502 real(dp),intent(out) :: pl_d2(mpsang) 1503 1504 !Local variables --------------------------------------- 1505 !scalars 1506 integer :: il,ilm 1507 real(dp) :: il_,il_m1,il_2m1 1508 character(len=500) :: msg 1509 !arrays 1510 real(dp) :: pl(mpsang),pl_d1(mpsang) 1511 1512 ! ********************************************************************* 1513 1514 if (abs(xx).gt.1.d0) then 1515 msg = 'pl_deriv : xx > 1 !' 1516 MSG_ERROR(msg) 1517 end if 1518 1519 pl_d2=zero; pl_d1=zero; pl=zero 1520 pl(1)=one; pl(2)=xx 1521 pl_d1(1)=zero; pl_d1(2)=one 1522 pl_d2(1)=zero; pl_d2(2)=zero 1523 if (mpsang>2) then 1524 do il=2,mpsang-1 1525 il_=dble(il);il_m1=dble(il-1);il_2m1=dble(2*il-1) 1526 ilm=il+1 1527 pl(ilm)=(il_2m1*xx*pl(ilm-1)-il_m1*pl(ilm-2))/il_ 1528 pl_d1(ilm)=(il_2m1*(xx*pl_d1(ilm-1)+pl(ilm-1))-il_m1*pl_d1(ilm-2))/il_ 1529 pl_d2(ilm)=(il_2m1*(xx*pl_d2(ilm-1)+two*pl_d1(ilm-1))-il_m1*pl_d2(ilm-2))/il_ 1530 end do 1531 end if 1532 1533 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
PARENTS
initylmg,m_paw_sphharm
CHILDREN
wrtout
SOURCE
1002 subroutine plm_coeff(blm,mpsang,xx) 1003 1004 1005 !This section has been created automatically by the script Abilint (TD). 1006 !Do not modify the following lines by hand. 1007 #undef ABI_FUNC 1008 #define ABI_FUNC 'plm_coeff' 1009 !End of the abilint section 1010 1011 implicit none 1012 1013 !Arguments --------------------------------------------- 1014 !scalars 1015 integer,intent(in) :: mpsang 1016 real(dp),intent(in) :: xx 1017 !arrays 1018 real(dp),intent(out) :: blm(5,mpsang*mpsang) 1019 1020 !Local variables --------------------------------------- 1021 !scalars 1022 integer :: il,ilm,ilm0,ilm1,im 1023 real(dp) :: dplm_dt,d2plm_dt2,llp1,onemx2,plm,sqrx,xsqrx,xx2,yy 1024 logical :: is_one 1025 character(len=500) :: msg 1026 !arrays 1027 real(dp) :: pl_d2(mpsang),plm_d2t(mpsang*mpsang) 1028 1029 !************************************************************************ 1030 1031 if (abs(xx).gt.1.d0) then 1032 msg = ' plm_coeff : xx > 1 !' 1033 MSG_ERROR(msg) 1034 end if 1035 1036 blm=zero 1037 is_one=(abs(abs(xx)-one)<=tol12) 1038 xx2=xx**2 1039 onemx2=abs(one-xx2) 1040 sqrx=sqrt(onemx2) 1041 xsqrx=xx*sqrt(onemx2) 1042 1043 call plm_d2theta(mpsang,plm_d2t,xx) 1044 if (is_one) then 1045 yy=sign(one,xx) 1046 call pl_deriv(mpsang,pl_d2,yy) 1047 end if 1048 1049 do il=0,mpsang-1 1050 llp1=dble(il*(il+1)) 1051 ilm0=il*il+il+1 1052 do im=0,il 1053 ilm=ilm0+im;ilm1=ilm0-im 1054 1055 plm =(-1)**im*ass_leg_pol(il,im,xx) 1056 dplm_dt =(-1)**im*plm_dtheta(il,im,xx) 1057 d2plm_dt2= plm_d2t(ilm) 1058 1059 blm(1,ilm)= two*xsqrx *dplm_dt+onemx2*d2plm_dt2 1060 blm(2,ilm)= (one-two*xx2)*dplm_dt-xsqrx *d2plm_dt2 1061 blm(3,ilm)=llp1*plm+ d2plm_dt2 1062 blm(4,ilm)= -two*xsqrx *dplm_dt+xx2 *d2plm_dt2 1063 1064 1065 if (is_one) then 1066 if (im==1) then 1067 blm(5,ilm)=llp1*plm+d2plm_dt2 1068 end if 1069 if (im==2) then 1070 blm(5,ilm)=d2plm_dt2-three*pl_d2(il+1) 1071 end if 1072 else 1073 if(im>0) then 1074 blm(5,ilm)=plm/onemx2-dplm_dt*xx/sqrx 1075 end if 1076 end if 1077 1078 if (im>0) then 1079 blm(1,ilm1)=blm(1,ilm) 1080 blm(2,ilm1)=blm(2,ilm) 1081 blm(3,ilm1)=blm(3,ilm) 1082 blm(4,ilm1)=blm(4,ilm) 1083 blm(5,ilm1)=blm(5,ilm) 1084 end if 1085 1086 end do 1087 end do 1088 1089 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)
PARENTS
m_paw_sphharm
CHILDREN
wrtout
SOURCE
1201 subroutine plm_d2theta(mpsang,plm_d2t,xx) 1202 1203 1204 !This section has been created automatically by the script Abilint (TD). 1205 !Do not modify the following lines by hand. 1206 #undef ABI_FUNC 1207 #define ABI_FUNC 'plm_d2theta' 1208 !End of the abilint section 1209 1210 implicit none 1211 1212 !Arguments --------------------------------------------- 1213 !scalars 1214 integer,intent(in) :: mpsang 1215 real(dp),intent(in) :: xx 1216 !arrays 1217 real(dp),intent(out) :: plm_d2t(mpsang*mpsang) 1218 1219 !Local variables --------------------------------------- 1220 !scalars 1221 integer :: il,ilm,ilmm1,ilmm2,im 1222 real(dp) :: sqrx 1223 character(len=500) :: msg 1224 1225 !************************************************************************ 1226 if (abs(xx).gt.1.d0) then 1227 msg = 'plm_d2theta : xx > 1 !' 1228 MSG_ERROR(msg) 1229 end if 1230 1231 plm_d2t=zero 1232 if (mpsang>1) then 1233 sqrx=sqrt(abs((1.d0-xx)*(1.d0+xx))) 1234 1235 do il=1,mpsang-1 1236 ilm=il*il+2*il+1 1237 ilmm1=(il-1)*(il-1)+2*(il-1)+1 1238 ! terme d2(Pll)/dtet2 1239 plm_d2t(ilm)=(2*il-1)*(sqrx*(plm_d2t(ilmm1)-(-1)**(il-1)*ass_leg_pol(il-1,il-1,xx))+& 1240 & 2.d0*xx*(-1)**(il-1)*plm_dtheta(il-1,il-1,xx)) 1241 plm_d2t(ilm-2*il)=plm_d2t(ilm) 1242 ! terme d2(Pl(l-1))/dtet2 1243 plm_d2t(ilm-1)=(2*il-1)*(xx*(plm_d2t(ilmm1)-(-1)**(il-1)*ass_leg_pol(il-1,il-1,xx))-& 1244 & 2.d0*sqrx*(-1)**(il-1)*plm_dtheta(il-1,il-1,xx)) 1245 if(il>1) plm_d2t(il*il+2)=plm_d2t(ilm-1) 1246 end do 1247 ! terme d2(Plm)/dtet2 1248 if(mpsang>2) then 1249 do il=2,mpsang-1 1250 do im=0,il-2 1251 ilm=il*il+il+1+im 1252 ilmm1=(il-1)*(il-1)+il+im 1253 ilmm2=(il-2)*(il-2)+il-1+im 1254 plm_d2t(ilm)=dble(2*il-1)/dble(il-im)*(xx*(plm_d2t(ilmm1)-(-1)**im*ass_leg_pol(il-1,im,xx))-& 1255 & 2.d0*sqrx*(-1)**im*plm_dtheta(il-1,im,xx))-& 1256 & dble(il+im-1)/dble(il-im)*plm_d2t(ilmm2) 1257 plm_d2t(ilm-2*im)=plm_d2t(ilm) 1258 end do 1259 end do 1260 end if 1261 end if 1262 1263 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))
PARENTS
CHILDREN
SOURCE
1293 function plm_dphi(ll,mm,xx) 1294 1295 1296 !This section has been created automatically by the script Abilint (TD). 1297 !Do not modify the following lines by hand. 1298 #undef ABI_FUNC 1299 #define ABI_FUNC 'plm_dphi' 1300 !End of the abilint section 1301 1302 implicit none 1303 1304 !Arguments --------------------------------------------- 1305 !scalars 1306 integer,intent(in) :: ll,mm 1307 real(dp) :: plm_dphi 1308 real(dp),intent(in) :: xx 1309 1310 !Local variables --------------------------------------- 1311 !scalars 1312 integer :: il,im 1313 real(dp) :: dosomx2,fact,pll,pmm,pmmp1,somx2 1314 character(len=500) :: msg 1315 1316 ! ********************************************************************* 1317 1318 if (mm.lt.0.or.mm.gt.ll.or.abs(xx).gt.1.d0) then 1319 msg = 'plm_dphi : mm < 0 or mm > ll or xx > 1 !' 1320 MSG_ERROR(msg) 1321 end if 1322 1323 plm_dphi=zero 1324 if (mm==0) return 1325 1326 pmm=one 1327 dosomx2=one 1328 if (mm > 0) then 1329 somx2=sqrt((1-xx)*(1+xx)) 1330 fact=one 1331 do im=1,mm 1332 pmm=-pmm*fact 1333 fact=fact+2 1334 end do 1335 if (mm > 1) then 1336 do im=2,mm 1337 dosomx2=somx2*dosomx2 1338 end do 1339 end if 1340 pmm=pmm*dosomx2 !due to one more term (-1^M) 1341 end if 1342 if(ll==mm) then 1343 plm_dphi=pmm*mm 1344 else 1345 pmmp1=xx*(2*mm+1)*pmm 1346 if(ll==mm+1) then 1347 plm_dphi=pmmp1*mm 1348 else if(ll>=mm+2) then 1349 do il=mm+2,ll 1350 pll=(xx*(2*il-1)*pmmp1-(il+mm-1)*pmm)/(il-mm) 1351 pmm=pmmp1 1352 pmmp1=pll 1353 end do 1354 plm_dphi=pll*mm 1355 end if 1356 end if 1357 1358 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))
PARENTS
CHILDREN
SOURCE
1388 function plm_dtheta(ll,mm,xx) 1389 1390 1391 !This section has been created automatically by the script Abilint (TD). 1392 !Do not modify the following lines by hand. 1393 #undef ABI_FUNC 1394 #define ABI_FUNC 'plm_dtheta' 1395 !End of the abilint section 1396 1397 implicit none 1398 1399 !Arguments --------------------------------------------- 1400 !scalars 1401 integer,intent(in) :: ll,mm 1402 real(dp) :: plm_dtheta 1403 real(dp),intent(in) :: xx 1404 1405 !Local variables --------------------------------------- 1406 !scalars 1407 integer :: il,im 1408 real(dp) :: dosomx2,dpll,dpmm,dpmmp1,fact,pll,pmm,pmmp1,somx2 1409 character(len=500) :: msg 1410 1411 ! ********************************************************************* 1412 1413 if (mm.lt.0.or.mm.gt.ll.or.abs(xx).gt.1.d0) then 1414 msg = 'plm_dtheta : mm < 0 or mm > ll or xx > 1 !' 1415 MSG_ERROR(msg) 1416 end if 1417 1418 plm_dtheta=zero 1419 pmm=one 1420 dpmm=one 1421 dosomx2=one 1422 somx2=sqrt((1-xx)*(1+xx)) 1423 if(mm==0)then 1424 dpmm=zero 1425 elseif (mm > 0) then 1426 fact=one 1427 do im=1,mm 1428 pmm=-pmm*fact*somx2 1429 dpmm=-dpmm*fact 1430 fact=fact+2 1431 end do 1432 if(mm>1)then 1433 do im=2,mm 1434 dosomx2=dosomx2*somx2 1435 end do 1436 end if 1437 dpmm= dpmm*mm*xx*dosomx2 1438 end if 1439 if(ll==mm)then 1440 plm_dtheta=dpmm 1441 else 1442 pmmp1=xx*(2*mm+1)*pmm 1443 dpmmp1=-(2*mm+1)*somx2*pmm+xx*(2*mm+1)*dpmm 1444 if(ll==mm+1) then 1445 plm_dtheta=dpmmp1 1446 else if(ll>=mm+2)then 1447 do il=mm+2,ll 1448 pll=(xx*(2*il-1)*pmmp1-(il+mm-1)*pmm)/(il-mm) 1449 dpll=(-somx2*(2*il-1)*pmmp1+(xx*(2*il-1)*dpmmp1-(il+mm-1)*dpmm))/(il-mm) 1450 pmm=pmmp1 1451 pmmp1=pll 1452 dpmm=dpmmp1 1453 dpmmp1=dpll 1454 end do 1455 plm_dtheta=dpll 1456 end if 1457 end if 1458 1459 end function plm_dtheta
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
PARENTS
CHILDREN
SOURCE
2687 subroutine setnabla_ylm(ang_phipphj,mpsang) 2688 2689 2690 !This section has been created automatically by the script Abilint (TD). 2691 !Do not modify the following lines by hand. 2692 #undef ABI_FUNC 2693 #define ABI_FUNC 'setnabla_ylm' 2694 !End of the abilint section 2695 2696 implicit none 2697 2698 !Arguments ------------------------------------ 2699 !scalars 2700 integer,intent(in) :: mpsang 2701 !arrays 2702 real(dp),intent(out) :: ang_phipphj(mpsang**2,mpsang**2,8) 2703 2704 !Local variables------------------------------- 2705 character(len=500) :: msg 2706 real(dp) :: ang_phipphj_tmp(16,16,8) 2707 2708 ! ************************************************************************ 2709 2710 if (mpsang>4) then 2711 msg=' Not designed for angular momentum greater than 3 !' 2712 MSG_ERROR(msg) 2713 end if 2714 2715 !8 angular integrals for l=0..3, m=-l..+l 2716 !ang_phipphj(1,4,1)=\frac{1}{\sqrt{3}} 2717 !ang_phipphj(2,5,1)=\frac{1}{\sqrt{5}} 2718 !ang_phipphj(3,8,1)=\frac{1}{\sqrt{5}} 2719 !ang_phipphj(4,1,1)=\frac{1}{\sqrt{3}} 2720 !ang_phipphj(4,7,1)=-\frac{1}{\sqrt{15}} 2721 !ang_phipphj(4,9,1)=\frac{1}{\sqrt{5}} 2722 !ang_phipphj(5,2,1)=\frac{1}{\sqrt{5}} 2723 !ang_phipphj(5,10,1)=\sqrt{\frac{3}{14}} 2724 !ang_phipphj(5,12,1)=-\frac{1}{\sqrt{70}} 2725 !ang_phipphj(6,11,1)=\frac{1}{\sqrt{7}} 2726 !ang_phipphj(7,4,1)=-\frac{1}{\sqrt{15}} 2727 !ang_phipphj(7,14,1)=\sqrt{\frac{6}{35}} 2728 !ang_phipphj(8,3,1)=\frac{1}{\sqrt{5}} 2729 !ang_phipphj(8,13,1)=-\sqrt{\frac{3}{35}} 2730 !ang_phipphj(8,15,1)=\frac{1}{\sqrt{7}} 2731 !ang_phipphj(9,4,1)=\frac{1}{\sqrt{5}} 2732 !ang_phipphj(9,14,1)=-\frac{1}{\sqrt{70}} 2733 !ang_phipphj(9,16,1)=\sqrt{\frac{3}{14}} 2734 !ang_phipphj(10,5,1)=\sqrt{\frac{3}{14}} 2735 !ang_phipphj(11,6,1)=\frac{1}{\sqrt{7}} 2736 !ang_phipphj(12,5,1)=-\frac{1}{\sqrt{70}} 2737 !ang_phipphj(13,8,1)=-\sqrt{\frac{3}{35}} 2738 !ang_phipphj(14,7,1)=\sqrt{\frac{6}{35}} 2739 !ang_phipphj(14,9,1)=-\frac{1}{\sqrt{70}} 2740 !ang_phipphj(15,8,1)=\frac{1}{\sqrt{7}} 2741 !ang_phipphj(16,9,1)=\sqrt{\frac{3}{14}} 2742 !ang_phipphj(1,4,2)=\frac{1}{2 \sqrt{3}} 2743 !ang_phipphj(1,14,2)=-\frac{\sqrt{\frac{7}{6}}}{2} 2744 !ang_phipphj(2,5,2)=\frac{1}{2 \sqrt{5}} 2745 !ang_phipphj(3,8,2)=\frac{1}{2 \sqrt{5}} 2746 !ang_phipphj(4,7,2)=-\sqrt{\frac{3}{5}} 2747 !ang_phipphj(4,9,2)=\frac{1}{2 \sqrt{5}} 2748 !ang_phipphj(5,2,2)=\frac{1}{4 \sqrt{5}} 2749 !ang_phipphj(5,10,2)=\frac{\sqrt{\frac{3}{14}}}{2} 2750 !ang_phipphj(5,12,2)=-2 \sqrt{\frac{2}{35}} 2751 !ang_phipphj(6,11,2)=\frac{1}{2 \sqrt{7}} 2752 !ang_phipphj(7,4,2)=\frac{1}{\sqrt{15}} 2753 !ang_phipphj(7,14,2)=\frac{13}{2 \sqrt{210}} 2754 !ang_phipphj(8,3,2)=-\frac{1}{\sqrt{5}} 2755 !ang_phipphj(8,13,2)=-4 \sqrt{\frac{3}{35}} 2756 !ang_phipphj(8,15,2)=\frac{1}{2 \sqrt{7}} 2757 !ang_phipphj(9,4,2)=\frac{1}{4 \sqrt{5}} 2758 !ang_phipphj(9,14,2)=-2 \sqrt{\frac{2}{35}} 2759 !ang_phipphj(9,16,2)=\frac{\sqrt{\frac{3}{14}}}{2} 2760 !ang_phipphj(10,5,2)=\frac{1}{\sqrt{42}} 2761 !ang_phipphj(11,6,2)=-\frac{1}{4 \sqrt{7}} 2762 !ang_phipphj(12,5,2)=\sqrt{\frac{2}{35}} 2763 !ang_phipphj(13,8,2)=2 \sqrt{\frac{3}{35}} 2764 !ang_phipphj(14,7,2)=-2 \sqrt{\frac{6}{35}} 2765 !ang_phipphj(14,9,2)=\sqrt{\frac{2}{35}} 2766 !ang_phipphj(15,8,2)=-\frac{1}{4 \sqrt{7}} 2767 !ang_phipphj(16,9,2)=\frac{1}{\sqrt{42}} 2768 !ang_phipphj(1,4,3)=\frac{\sqrt{3}}{2} 2769 !ang_phipphj(1,14,3)=\frac{\sqrt{\frac{7}{6}}}{2} 2770 !ang_phipphj(2,5,3)=\frac{\sqrt{5}}{2} 2771 !ang_phipphj(3,8,3)=\frac{\sqrt{5}}{2} 2772 !ang_phipphj(4,9,3)=\frac{\sqrt{5}}{2} 2773 !ang_phipphj(5,2,3)=-\frac{\sqrt{5}}{4} 2774 !ang_phipphj(5,10,3)=\frac{\sqrt{\frac{21}{2}}}{2} 2775 !ang_phipphj(6,11,3)=\frac{\sqrt{7}}{2} 2776 !ang_phipphj(7,14,3)=\frac{\sqrt{\frac{35}{6}}}{2} 2777 !ang_phipphj(8,15,3)=\frac{\sqrt{7}}{2} 2778 !ang_phipphj(9,4,3)=-\frac{\sqrt{5}}{4} 2779 !ang_phipphj(9,16,3)=\frac{\sqrt{\frac{21}{2}}}{2} 2780 !ang_phipphj(10,5,3)=-\sqrt{\frac{7}{6}} 2781 !ang_phipphj(11,6,3)=-\frac{\sqrt{7}}{4} 2782 !ang_phipphj(15,8,3)=-\frac{\sqrt{7}}{4} 2783 !ang_phipphj(16,9,3)=-\sqrt{\frac{7}{6}} 2784 !ang_phipphj(1,2,4)=\frac{1}{\sqrt{3}} 2785 !ang_phipphj(2,1,4)=\frac{1}{\sqrt{3}} 2786 !ang_phipphj(2,7,4)=-\frac{1}{\sqrt{15}} 2787 !ang_phipphj(2,9,4)=-\frac{1}{\sqrt{5}} 2788 !ang_phipphj(3,6,4)=\frac{1}{\sqrt{5}} 2789 !ang_phipphj(4,5,4)=\frac{1}{\sqrt{5}} 2790 !ang_phipphj(5,4,4)=\frac{1}{\sqrt{5}} 2791 !ang_phipphj(5,14,4)=-\frac{1}{\sqrt{70}} 2792 !ang_phipphj(5,16,4)=-\sqrt{\frac{3}{14}} 2793 !ang_phipphj(6,3,4)=\frac{1}{\sqrt{5}} 2794 !ang_phipphj(6,13,4)=-\sqrt{\frac{3}{35}} 2795 !ang_phipphj(6,15,4)=-\frac{1}{\sqrt{7}} 2796 !ang_phipphj(7,2,4)=-\frac{1}{\sqrt{15}} 2797 !ang_phipphj(7,12,4)=\sqrt{\frac{6}{35}} 2798 !ang_phipphj(8,11,4)=\frac{1}{\sqrt{7}} 2799 !ang_phipphj(9,2,4)=-\frac{1}{\sqrt{5}} 2800 !ang_phipphj(9,10,4)=\sqrt{\frac{3}{14}} 2801 !ang_phipphj(9,12,4)=\frac{1}{\sqrt{70}} 2802 !ang_phipphj(10,9,4)=\sqrt{\frac{3}{14}} 2803 !ang_phipphj(11,8,4)=\frac{1}{\sqrt{7}} 2804 !ang_phipphj(12,7,4)=\sqrt{\frac{6}{35}} 2805 !ang_phipphj(12,9,4)=\frac{1}{\sqrt{70}} 2806 !ang_phipphj(13,6,4)=-\sqrt{\frac{3}{35}} 2807 !ang_phipphj(14,5,4)=-\frac{1}{\sqrt{70}} 2808 !ang_phipphj(15,6,4)=-\frac{1}{\sqrt{7}} 2809 !ang_phipphj(16,5,4)=-\sqrt{\frac{3}{14}} 2810 !ang_phipphj(1,2,5)=\frac{1}{2 \sqrt{3}} 2811 !ang_phipphj(1,12,5)=-\frac{\sqrt{\frac{7}{6}}}{2} 2812 !ang_phipphj(2,7,5)=-\sqrt{\frac{3}{5}} 2813 !ang_phipphj(2,9,5)=-\frac{1}{2 \sqrt{5}} 2814 !ang_phipphj(3,6,5)=\frac{1}{2 \sqrt{5}} 2815 !ang_phipphj(4,5,5)=\frac{1}{2 \sqrt{5}} 2816 !ang_phipphj(5,4,5)=\frac{1}{4 \sqrt{5}} 2817 !ang_phipphj(5,14,5)=-2 \sqrt{\frac{2}{35}} 2818 !ang_phipphj(5,16,5)=-\frac{\sqrt{\frac{3}{14}}}{2} 2819 !ang_phipphj(6,3,5)=-\frac{1}{\sqrt{5}} 2820 !ang_phipphj(6,13,5)=-4 \sqrt{\frac{3}{35}} 2821 !ang_phipphj(6,15,5)=-\frac{1}{2 \sqrt{7}} 2822 !ang_phipphj(7,2,5)=\frac{1}{\sqrt{15}} 2823 !ang_phipphj(7,12,5)=\frac{13}{2 \sqrt{210}} 2824 !ang_phipphj(8,11,5)=\frac{1}{2 \sqrt{7}} 2825 !ang_phipphj(9,2,5)=-\frac{1}{4 \sqrt{5}} 2826 !ang_phipphj(9,10,5)=\frac{\sqrt{\frac{3}{14}}}{2} 2827 !ang_phipphj(9,12,5)=2 \sqrt{\frac{2}{35}} 2828 !ang_phipphj(10,9,5)=\frac{1}{\sqrt{42}} 2829 !ang_phipphj(11,8,5)=-\frac{1}{4 \sqrt{7}} 2830 !ang_phipphj(12,7,5)=-2 \sqrt{\frac{6}{35}} 2831 !ang_phipphj(12,9,5)=-\sqrt{\frac{2}{35}} 2832 !ang_phipphj(13,6,5)=2 \sqrt{\frac{3}{35}} 2833 !ang_phipphj(14,5,5)=\sqrt{\frac{2}{35}} 2834 !ang_phipphj(15,6,5)=\frac{1}{4 \sqrt{7}} 2835 !ang_phipphj(16,5,5)=-\frac{1}{\sqrt{42}} 2836 !ang_phipphj(1,2,6)=\frac{\sqrt{3}}{2} 2837 !ang_phipphj(1,12,6)=\frac{\sqrt{\frac{7}{6}}}{2} 2838 !ang_phipphj(2,9,6)=-\frac{\sqrt{5}}{2} 2839 !ang_phipphj(3,6,6)=\frac{\sqrt{5}}{2} 2840 !ang_phipphj(4,5,6)=\frac{\sqrt{5}}{2} 2841 !ang_phipphj(5,4,6)=-\frac{\sqrt{5}}{4} 2842 !ang_phipphj(5,16,6)=-\frac{\sqrt{\frac{21}{2}}}{2} 2843 !ang_phipphj(6,15,6)=-\frac{\sqrt{7}}{2} 2844 !ang_phipphj(7,12,6)=\frac{\sqrt{\frac{35}{6}}}{2} 2845 !ang_phipphj(8,11,6)=\frac{\sqrt{7}}{2} 2846 !ang_phipphj(9,2,6)=\frac{\sqrt{5}}{4} 2847 !ang_phipphj(9,10,6)=\frac{\sqrt{\frac{21}{2}}}{2} 2848 !ang_phipphj(10,9,6)=-\sqrt{\frac{7}{6}} 2849 !ang_phipphj(11,8,6)=-\frac{\sqrt{7}}{4} 2850 !ang_phipphj(15,6,6)=\frac{\sqrt{7}}{4} 2851 !ang_phipphj(16,5,6)=\sqrt{\frac{7}{6}} 2852 !ang_phipphj(1,3,7)=\frac{1}{\sqrt{3}} 2853 !ang_phipphj(2,6,7)=\frac{1}{\sqrt{5}} 2854 !ang_phipphj(3,1,7)=\frac{1}{\sqrt{3}} 2855 !ang_phipphj(3,7,7)=\frac{2}{\sqrt{15}} 2856 !ang_phipphj(4,8,7)=\frac{1}{\sqrt{5}} 2857 !ang_phipphj(5,11,7)=\frac{1}{\sqrt{7}} 2858 !ang_phipphj(6,2,7)=\frac{1}{\sqrt{5}} 2859 !ang_phipphj(6,12,7)=2 \sqrt{\frac{2}{35}} 2860 !ang_phipphj(7,3,7)=\frac{2}{\sqrt{15}} 2861 !ang_phipphj(7,13,7)=\frac{3}{\sqrt{35}} 2862 !ang_phipphj(8,4,7)=\frac{1}{\sqrt{5}} 2863 !ang_phipphj(8,14,7)=2 \sqrt{\frac{2}{35}} 2864 !ang_phipphj(9,15,7)=\frac{1}{\sqrt{7}} 2865 !ang_phipphj(11,5,7)=\frac{1}{\sqrt{7}} 2866 !ang_phipphj(12,6,7)=2 \sqrt{\frac{2}{35}} 2867 !ang_phipphj(13,7,7)=\frac{3}{\sqrt{35}} 2868 !ang_phipphj(14,8,7)=2 \sqrt{\frac{2}{35}} 2869 !ang_phipphj(15,9,7)=\frac{1}{\sqrt{7}} 2870 !ang_phipphj(1,3,8)=\frac{2}{\sqrt{3}} 2871 !ang_phipphj(2,6,8)=\frac{3}{\sqrt{5}} 2872 !ang_phipphj(3,7,8)=2 \sqrt{\frac{3}{5}} 2873 !ang_phipphj(4,8,8)=\frac{3}{\sqrt{5}} 2874 !ang_phipphj(5,11,8)=\frac{4}{\sqrt{7}} 2875 !ang_phipphj(6,2,8)=-\frac{1}{\sqrt{5}} 2876 !ang_phipphj(6,12,8)=8 \sqrt{\frac{2}{35}} 2877 !ang_phipphj(7,3,8)=-\frac{2}{\sqrt{15}} 2878 !ang_phipphj(7,13,8)=\frac{12}{\sqrt{35}} 2879 !ang_phipphj(8,4,8)=-\frac{1}{\sqrt{5}} 2880 !ang_phipphj(8,14,8)=8 \sqrt{\frac{2}{35}} 2881 !ang_phipphj(9,15,8)=\frac{4}{\sqrt{7}} 2882 !ang_phipphj(11,5,8)=-\frac{2}{\sqrt{7}} 2883 !ang_phipphj(12,6,8)=-4 \sqrt{\frac{2}{35}} 2884 !ang_phipphj(13,7,8)=-\frac{6}{\sqrt{35}} 2885 !ang_phipphj(14,8,8)=-4 \sqrt{\frac{2}{35}} 2886 !ang_phipphj(15,9,8)=-\frac{2}{\sqrt{7}} 2887 2888 2889 ang_phipphj_tmp=zero 2890 ! 2891 ang_phipphj_tmp(1,4,1)=0.57735026918962576451_dp 2892 ang_phipphj_tmp(2,5,1)=0.44721359549995793928_dp 2893 ang_phipphj_tmp(3,8,1)=0.44721359549995793928_dp 2894 ang_phipphj_tmp(4,1,1)=0.57735026918962576451_dp 2895 ang_phipphj_tmp(4,7,1)=-0.25819888974716112568_dp 2896 ang_phipphj_tmp(4,9,1)=0.44721359549995793928_dp 2897 ang_phipphj_tmp(5,2,1)=0.44721359549995793928_dp 2898 ang_phipphj_tmp(5,10,1)=0.46291004988627573078_dp 2899 ang_phipphj_tmp(5,12,1)=-0.11952286093343936400_dp 2900 ang_phipphj_tmp(6,11,1)=0.37796447300922722721_dp 2901 ang_phipphj_tmp(7,4,1)=-0.25819888974716112568_dp 2902 ang_phipphj_tmp(7,14,1)=0.41403933560541253068_dp 2903 ang_phipphj_tmp(8,3,1)=0.44721359549995793928_dp 2904 ang_phipphj_tmp(8,13,1)=-0.29277002188455995381_dp 2905 ang_phipphj_tmp(8,15,1)=0.37796447300922722721_dp 2906 ang_phipphj_tmp(9,4,1)=0.44721359549995793928_dp 2907 ang_phipphj_tmp(9,14,1)=-0.11952286093343936400_dp 2908 ang_phipphj_tmp(9,16,1)=0.46291004988627573078_dp 2909 ang_phipphj_tmp(10,5,1)=0.46291004988627573078_dp 2910 ang_phipphj_tmp(11,6,1)=0.37796447300922722721_dp 2911 ang_phipphj_tmp(12,5,1)=-0.11952286093343936400_dp 2912 ang_phipphj_tmp(13,8,1)=-0.29277002188455995381_dp 2913 ang_phipphj_tmp(14,7,1)=0.41403933560541253068_dp 2914 ang_phipphj_tmp(14,9,1)=-0.11952286093343936400_dp 2915 ang_phipphj_tmp(15,8,1)=0.37796447300922722721_dp 2916 ang_phipphj_tmp(16,9,1)=0.46291004988627573078_dp 2917 ! 2918 ang_phipphj_tmp(1,4,2)=0.28867513459481288225_dp 2919 ang_phipphj_tmp(1,14,2)=-0.54006172486732168591_dp 2920 ang_phipphj_tmp(2,5,2)=0.22360679774997896964_dp 2921 ang_phipphj_tmp(3,8,2)=0.22360679774997896964_dp 2922 ang_phipphj_tmp(4,7,2)=-0.77459666924148337704_dp 2923 ang_phipphj_tmp(4,9,2)=0.22360679774997896964_dp 2924 ang_phipphj_tmp(5,2,2)=0.11180339887498948482_dp 2925 ang_phipphj_tmp(5,10,2)=0.23145502494313786539_dp 2926 ang_phipphj_tmp(5,12,2)=-0.47809144373375745599_dp 2927 ang_phipphj_tmp(6,11,2)=0.18898223650461361361_dp 2928 ang_phipphj_tmp(7,4,2)=0.25819888974716112568_dp 2929 ang_phipphj_tmp(7,14,2)=0.44854261357253024157_dp 2930 ang_phipphj_tmp(8,3,2)=-0.44721359549995793928_dp 2931 ang_phipphj_tmp(8,13,2)=-1.1710800875382398152_dp 2932 ang_phipphj_tmp(8,15,2)=0.18898223650461361361_dp 2933 ang_phipphj_tmp(9,4,2)=0.11180339887498948482_dp 2934 ang_phipphj_tmp(9,14,2)=-0.47809144373375745599_dp 2935 ang_phipphj_tmp(9,16,2)=0.23145502494313786539_dp 2936 ang_phipphj_tmp(10,5,2)=0.15430334996209191026_dp 2937 ang_phipphj_tmp(11,6,2)=-0.094491118252306806804_dp 2938 ang_phipphj_tmp(12,5,2)=0.23904572186687872799_dp 2939 ang_phipphj_tmp(13,8,2)=0.58554004376911990761_dp 2940 ang_phipphj_tmp(14,7,2)=-0.82807867121082506136_dp 2941 ang_phipphj_tmp(14,9,2)=0.23904572186687872799_dp 2942 ang_phipphj_tmp(15,8,2)=-0.094491118252306806804_dp 2943 ang_phipphj_tmp(16,9,2)=0.15430334996209191026_dp 2944 ! 2945 ang_phipphj_tmp(1,4,3)=0.86602540378443864676_dp 2946 ang_phipphj_tmp(1,14,3)=0.54006172486732168591_dp 2947 ang_phipphj_tmp(2,5,3)=1.1180339887498948482_dp 2948 ang_phipphj_tmp(3,8,3)=1.1180339887498948482_dp 2949 ang_phipphj_tmp(4,9,3)=1.1180339887498948482_dp 2950 ang_phipphj_tmp(5,2,3)=-0.55901699437494742410_dp 2951 ang_phipphj_tmp(5,10,3)=1.6201851746019650577_dp 2952 ang_phipphj_tmp(6,11,3)=1.3228756555322952953_dp 2953 ang_phipphj_tmp(7,14,3)=1.2076147288491198811_dp 2954 ang_phipphj_tmp(8,15,3)=1.3228756555322952953_dp 2955 ang_phipphj_tmp(9,4,3)=-0.55901699437494742410_dp 2956 ang_phipphj_tmp(9,16,3)=1.6201851746019650577_dp 2957 ang_phipphj_tmp(10,5,3)=-1.0801234497346433718_dp 2958 ang_phipphj_tmp(11,6,3)=-0.66143782776614764763_dp 2959 ang_phipphj_tmp(15,8,3)=-0.66143782776614764763_dp 2960 ang_phipphj_tmp(16,9,3)=-1.0801234497346433718_dp 2961 ! 2962 ang_phipphj_tmp(1,2,4)=0.57735026918962576451_dp 2963 ang_phipphj_tmp(2,1,4)=0.57735026918962576451_dp 2964 ang_phipphj_tmp(2,7,4)=-0.25819888974716112568_dp 2965 ang_phipphj_tmp(2,9,4)=-0.44721359549995793928_dp 2966 ang_phipphj_tmp(3,6,4)=0.44721359549995793928_dp 2967 ang_phipphj_tmp(4,5,4)=0.44721359549995793928_dp 2968 ang_phipphj_tmp(5,4,4)=0.44721359549995793928_dp 2969 ang_phipphj_tmp(5,14,4)=-0.11952286093343936400_dp 2970 ang_phipphj_tmp(5,16,4)=-0.46291004988627573078_dp 2971 ang_phipphj_tmp(6,3,4)=0.44721359549995793928_dp 2972 ang_phipphj_tmp(6,13,4)=-0.29277002188455995381_dp 2973 ang_phipphj_tmp(6,15,4)=-0.37796447300922722721_dp 2974 ang_phipphj_tmp(7,2,4)=-0.25819888974716112568_dp 2975 ang_phipphj_tmp(7,12,4)=0.41403933560541253068_dp 2976 ang_phipphj_tmp(8,11,4)=0.37796447300922722721_dp 2977 ang_phipphj_tmp(9,2,4)=-0.44721359549995793928_dp 2978 ang_phipphj_tmp(9,10,4)=0.46291004988627573078_dp 2979 ang_phipphj_tmp(9,12,4)=0.11952286093343936400_dp 2980 ang_phipphj_tmp(10,9,4)=0.46291004988627573078_dp 2981 ang_phipphj_tmp(11,8,4)=0.37796447300922722721_dp 2982 ang_phipphj_tmp(12,7,4)=0.41403933560541253068_dp 2983 ang_phipphj_tmp(12,9,4)=0.11952286093343936400_dp 2984 ang_phipphj_tmp(13,6,4)=-0.29277002188455995381_dp 2985 ang_phipphj_tmp(14,5,4)=-0.11952286093343936400_dp 2986 ang_phipphj_tmp(15,6,4)=-0.37796447300922722721_dp 2987 ang_phipphj_tmp(16,5,4)=-0.46291004988627573078_dp 2988 ! 2989 ang_phipphj_tmp(1,2,5)=0.28867513459481288225_dp 2990 ang_phipphj_tmp(1,12,5)=-0.54006172486732168591_dp 2991 ang_phipphj_tmp(2,7,5)=-0.77459666924148337704_dp 2992 ang_phipphj_tmp(2,9,5)=-0.22360679774997896964_dp 2993 ang_phipphj_tmp(3,6,5)=0.22360679774997896964_dp 2994 ang_phipphj_tmp(4,5,5)=0.22360679774997896964_dp 2995 ang_phipphj_tmp(5,4,5)=0.11180339887498948482_dp 2996 ang_phipphj_tmp(5,14,5)=-0.47809144373375745599_dp 2997 ang_phipphj_tmp(5,16,5)=-0.23145502494313786539_dp 2998 ang_phipphj_tmp(6,3,5)=-0.44721359549995793928_dp 2999 ang_phipphj_tmp(6,13,5)=-1.1710800875382398152_dp 3000 ang_phipphj_tmp(6,15,5)=-0.18898223650461361361_dp 3001 ang_phipphj_tmp(7,2,5)=0.25819888974716112568_dp 3002 ang_phipphj_tmp(7,12,5)=0.44854261357253024157_dp 3003 ang_phipphj_tmp(8,11,5)=0.18898223650461361361_dp 3004 ang_phipphj_tmp(9,2,5)=-0.11180339887498948482_dp 3005 ang_phipphj_tmp(9,10,5)=0.23145502494313786539_dp 3006 ang_phipphj_tmp(9,12,5)=0.47809144373375745599_dp 3007 ang_phipphj_tmp(10,9,5)=0.15430334996209191026_dp 3008 ang_phipphj_tmp(11,8,5)=-0.094491118252306806804_dp 3009 ang_phipphj_tmp(12,7,5)=-0.82807867121082506136_dp 3010 ang_phipphj_tmp(12,9,5)=-0.23904572186687872799_dp 3011 ang_phipphj_tmp(13,6,5)=0.58554004376911990761_dp 3012 ang_phipphj_tmp(14,5,5)=0.23904572186687872799_dp 3013 ang_phipphj_tmp(15,6,5)=0.094491118252306806804_dp 3014 ang_phipphj_tmp(16,5,5)=-0.15430334996209191026_dp 3015 ! 3016 ang_phipphj_tmp(1,2,6)=0.86602540378443864676_dp 3017 ang_phipphj_tmp(1,12,6)=0.54006172486732168591_dp 3018 ang_phipphj_tmp(2,9,6)=-1.1180339887498948482_dp 3019 ang_phipphj_tmp(3,6,6)=1.1180339887498948482_dp 3020 ang_phipphj_tmp(4,5,6)=1.1180339887498948482_dp 3021 ang_phipphj_tmp(5,4,6)=-0.55901699437494742410_dp 3022 ang_phipphj_tmp(5,16,6)=-1.6201851746019650577_dp 3023 ang_phipphj_tmp(6,15,6)=-1.3228756555322952953_dp 3024 ang_phipphj_tmp(7,12,6)=1.2076147288491198811_dp 3025 ang_phipphj_tmp(8,11,6)=1.3228756555322952953_dp 3026 ang_phipphj_tmp(9,2,6)=0.55901699437494742410_dp 3027 ang_phipphj_tmp(9,10,6)=1.6201851746019650577_dp 3028 ang_phipphj_tmp(10,9,6)=-1.0801234497346433718_dp 3029 ang_phipphj_tmp(11,8,6)=-0.66143782776614764763_dp 3030 ang_phipphj_tmp(15,6,6)=0.66143782776614764763_dp 3031 ang_phipphj_tmp(16,5,6)=1.0801234497346433718_dp 3032 ! 3033 ang_phipphj_tmp(1,3,7)=0.57735026918962576451_dp 3034 ang_phipphj_tmp(2,6,7)=0.44721359549995793928_dp 3035 ang_phipphj_tmp(3,1,7)=0.57735026918962576451_dp 3036 ang_phipphj_tmp(3,7,7)=0.51639777949432225136_dp 3037 ang_phipphj_tmp(4,8,7)=0.44721359549995793928_dp 3038 ang_phipphj_tmp(5,11,7)=0.37796447300922722721_dp 3039 ang_phipphj_tmp(6,2,7)=0.44721359549995793928_dp 3040 ang_phipphj_tmp(6,12,7)=0.47809144373375745599_dp 3041 ang_phipphj_tmp(7,3,7)=0.51639777949432225136_dp 3042 ang_phipphj_tmp(7,13,7)=0.50709255283710994651_dp 3043 ang_phipphj_tmp(8,4,7)=0.44721359549995793928_dp 3044 ang_phipphj_tmp(8,14,7)=0.47809144373375745599_dp 3045 ang_phipphj_tmp(9,15,7)=0.37796447300922722721_dp 3046 ang_phipphj_tmp(11,5,7)=0.37796447300922722721_dp 3047 ang_phipphj_tmp(12,6,7)=0.47809144373375745599_dp 3048 ang_phipphj_tmp(13,7,7)=0.50709255283710994651_dp 3049 ang_phipphj_tmp(14,8,7)=0.47809144373375745599_dp 3050 ang_phipphj_tmp(15,9,7)=0.37796447300922722721_dp 3051 ! 3052 ang_phipphj_tmp(1,3,8)=1.1547005383792515290_dp 3053 ang_phipphj_tmp(2,6,8)=1.3416407864998738178_dp 3054 ang_phipphj_tmp(3,7,8)=1.5491933384829667541_dp 3055 ang_phipphj_tmp(4,8,8)=1.3416407864998738178_dp 3056 ang_phipphj_tmp(5,11,8)=1.5118578920369089089_dp 3057 ang_phipphj_tmp(6,2,8)=-0.44721359549995793928_dp 3058 ang_phipphj_tmp(6,12,8)=1.9123657749350298240_dp 3059 ang_phipphj_tmp(7,3,8)=-0.51639777949432225136_dp 3060 ang_phipphj_tmp(7,13,8)=2.0283702113484397860_dp 3061 ang_phipphj_tmp(8,4,8)=-0.44721359549995793928_dp 3062 ang_phipphj_tmp(8,14,8)=1.9123657749350298240_dp 3063 ang_phipphj_tmp(9,15,8)=1.5118578920369089089_dp 3064 ang_phipphj_tmp(11,5,8)=-0.75592894601845445443_dp 3065 ang_phipphj_tmp(12,6,8)=-0.95618288746751491198_dp 3066 ang_phipphj_tmp(13,7,8)=-1.0141851056742198930_dp 3067 ang_phipphj_tmp(14,8,8)=-0.95618288746751491198_dp 3068 ang_phipphj_tmp(15,9,8)=-0.75592894601845445443_dp 3069 3070 ang_phipphj(:,:,:)=ang_phipphj_tmp(1:mpsang**2,1:mpsang**2,:) 3071 3072 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
PARENTS
CHILDREN
SOURCE
2527 subroutine setsym_ylm(gprimd,lmax,nsym,pawprtvol,rprimd,sym,zarot) 2528 2529 2530 !This section has been created automatically by the script Abilint (TD). 2531 !Do not modify the following lines by hand. 2532 #undef ABI_FUNC 2533 #define ABI_FUNC 'setsym_ylm' 2534 !End of the abilint section 2535 2536 implicit none 2537 2538 !Arguments --------------------------------------------- 2539 !scalars 2540 integer,intent(in) :: lmax,nsym,pawprtvol 2541 !arrays 2542 integer,intent(in) :: sym(3,3,nsym) 2543 real(dp),intent(in) :: gprimd(3,3),rprimd(3,3) 2544 real(dp),intent(out) :: zarot(2*lmax+1,2*lmax+1,lmax+1,nsym) 2545 2546 !Local variables ------------------------------ 2547 !scalars 2548 integer :: i1,ii,il,irot,isn,j1,jj,k1,ll,mm,mp 2549 real(dp) :: cosalp,cosbeta,cosgam,sinalp,singam 2550 character(len=1000) :: msg 2551 !arrays 2552 real(dp) :: prod(3,3),rot(3,3) 2553 !************************************************************************ 2554 2555 if (abs(pawprtvol)>=3) then 2556 write(msg,'(8a,i4)') ch10,& 2557 & ' PAW TEST:',ch10,& 2558 & ' ==== setsym_ylm: rotation matrices in the basis ============',ch10,& 2559 & ' ==== of real spherical harmonics ============',ch10,& 2560 & ' > Number of symmetries (nsym)=',nsym 2561 call wrtout(std_out,msg,'COLL') 2562 end if 2563 2564 zarot=zero 2565 2566 do irot=1,nsym 2567 2568 if (abs(pawprtvol)>=3) then 2569 write(msg,'(a,i2,a,9i2,a)') ' >For symmetry ',irot,' (',sym(:,:,irot),')' 2570 call wrtout(std_out,msg,'COLL') 2571 end if 2572 2573 ! === l=0 case === 2574 zarot(1,1,1,irot)=one 2575 2576 ! === l>0 case === 2577 if (lmax>0) then 2578 ! Calculate the rotations in the cartesian basis 2579 rot=zero;prod=zero 2580 do k1=1,3 2581 do j1=1,3 2582 do i1=1,3 2583 prod(i1,j1)=prod(i1,j1)+sym(i1,k1,irot)*rprimd(j1,k1) 2584 end do 2585 end do 2586 end do 2587 do j1=1,3 2588 do i1=1,3 2589 do k1=1,3 2590 rot(i1,j1)=rot(i1,j1)+gprimd(i1,k1)*prod(k1,j1) 2591 end do 2592 if(abs(rot(i1,j1))<tol10) rot(i1,j1)=zero 2593 end do 2594 end do 2595 call mkeuler(rot,cosbeta,cosalp,sinalp,cosgam,singam,isn) 2596 do ll=1,lmax 2597 il=(isn)**ll 2598 do mp=-ll,ll 2599 jj=mp+ll+1 2600 do mm=-ll,ll 2601 ii=mm+ll+1 2602 2603 ! Formula (47) from the paper of Blanco et al 2604 zarot(ii,jj,ll+1,irot)=il& 2605 & *(phim(cosalp,sinalp,mm)*phim(cosgam,singam,mp)*sign(1,mp)& 2606 *(dbeta(cosbeta,ll,abs(mp),abs(mm))& 2607 & +(-1._dp)**mm*dbeta(cosbeta,ll,abs(mm),-abs(mp)))*half& 2608 & -phim(cosalp,sinalp,-mm)*phim(cosgam,singam,-mp)*sign(1,mm)& 2609 *(dbeta(cosbeta,ll,abs(mp),abs(mm))& 2610 & -(-1._dp)**mm*dbeta(cosbeta,ll,abs(mm),-abs(mp)))*half) 2611 end do 2612 end do 2613 end do 2614 end if ! lmax case 2615 2616 if (abs(pawprtvol)>=3) then 2617 if(lmax>0) then 2618 write(msg,'(2a,3(3(2x,f7.3),a))') & 2619 & ' Rotation matrice for l=1:',ch10,& 2620 & (zarot(1,jj,2,irot),jj=1,3),ch10,& 2621 & (zarot(2,jj,2,irot),jj=1,3),ch10,& 2622 & (zarot(3,jj,2,irot),jj=1,3) 2623 call wrtout(std_out,msg,'COLL') 2624 end if 2625 if(lmax>1) then 2626 write(msg,'(2a,5(5(2x,f7.3),a))') & 2627 & ' Rotation matrice for l=2:',ch10,& 2628 & (zarot(1,jj,3,irot),jj=1,5),ch10,& 2629 & (zarot(2,jj,3,irot),jj=1,5),ch10,& 2630 & (zarot(3,jj,3,irot),jj=1,5),ch10,& 2631 & (zarot(4,jj,3,irot),jj=1,5),ch10,& 2632 & (zarot(5,jj,3,irot),jj=1,5) 2633 call wrtout(std_out,msg,'COLL') 2634 end if 2635 if(lmax>2) then 2636 write(msg,'(2a,7(7(2x,f7.3),a))') & 2637 & ' Rotation matrice for l=3:',ch10,& 2638 & (zarot(1,jj,4,irot),jj=1,7),ch10,& 2639 & (zarot(2,jj,4,irot),jj=1,7),ch10,& 2640 & (zarot(3,jj,4,irot),jj=1,7),ch10,& 2641 & (zarot(4,jj,4,irot),jj=1,7),ch10,& 2642 & (zarot(5,jj,4,irot),jj=1,7),ch10,& 2643 & (zarot(6,jj,4,irot),jj=1,7),ch10,& 2644 & (zarot(7,jj,4,irot),jj=1,7) 2645 call wrtout(std_out,msg,'COLL') 2646 end if 2647 end if 2648 2649 end do ! isym loop 2650 2651 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>
PARENTS
m_pawdij
CHILDREN
wrtout
SOURCE
936 subroutine slxyzs(lp,mp,idir,ll,mm,sls_val) 937 938 939 !This section has been created automatically by the script Abilint (TD). 940 !Do not modify the following lines by hand. 941 #undef ABI_FUNC 942 #define ABI_FUNC 'slxyzs' 943 !End of the abilint section 944 945 implicit none 946 947 !Arguments --------------------------------------------- 948 !scalars 949 integer,intent(in) :: idir,ll,lp,mm,mp 950 complex(dpc),intent(out) :: sls_val 951 952 !Local variables --------------------------------------- 953 !scalars 954 integer :: lpp,lppp,mpp,mppp 955 complex(dpc) :: lidir,sy_val,ys_val 956 957 ! ********************************************************************* 958 959 sls_val = czero 960 961 if (lp == ll) then 962 lpp = ll 963 lppp = ll 964 do mpp = -lpp, lpp 965 call ys(lpp,mpp,lp,mp,sy_val) 966 do mppp = -lppp, lppp 967 call lxyz(lpp,mpp,idir,lppp,mppp,lidir) 968 call ys(lppp,mppp,ll,mm,ys_val) 969 sls_val = sls_val + conjg(sy_val)*lidir*ys_val 970 end do 971 end do 972 end if 973 974 end subroutine slxyzs
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
PARENTS
mlwfovlp_proj,mlwfovlp_ylmfac
CHILDREN
wrtout
SOURCE
414 subroutine ylm_cmplx(lx,ylm,xx,yy,zz) 415 416 417 !This section has been created automatically by the script Abilint (TD). 418 !Do not modify the following lines by hand. 419 #undef ABI_FUNC 420 #define ABI_FUNC 'ylm_cmplx' 421 !End of the abilint section 422 423 implicit none 424 425 !Arguments ------------------------------------ 426 !scalars 427 integer,intent(in) :: lx 428 real(dp),intent(in) :: xx,yy,zz 429 !arrays 430 complex(dpc),intent(out) :: ylm((lx+1)*(lx+1)) 431 432 !Local variables------------------------------- 433 !scalars 434 integer :: ii,l1,m1,nc,nn 435 real(dp) :: dc,dl,dm,ds,rr,rrs,rs,sq2,w,x,xs,ya,yi,yr 436 !arrays 437 real(dp) :: cosa(lx+1),fact(2*(lx+1)),plm(lx+2,lx+2),qlm(lx+2,lx+2),sgn(lx+1) 438 real(dp) :: sina(lx+1) 439 440 ! ************************************************************************* 441 442 !normalization coefficients 443 sq2=sqrt(2.0d0) 444 fact(1)=1.0d0 445 do ii=2,2*lx+1 446 fact(ii)=(ii-1)*fact(ii-1) 447 end do 448 do l1=1,lx+1 449 sgn(l1)=(-1.d0)**(l1-1) 450 do m1=1,l1 451 qlm(l1,m1)=sqrt((2*l1-1)*fact(l1-m1+1)/& 452 & (four_pi*fact(l1+m1-1))) 453 end do 454 end do 455 456 !legendre polynomials 457 rs=xx**2 + yy**2 + zz**2 458 if(rs > tol8) then 459 xs=zz**2/rs 460 x=zz/sqrt(rs) 461 w=sqrt(abs(1.0d0 - xs)) 462 else 463 x=0.0d0 464 465 w=1.0d0 466 end if 467 plm(1,1)=1.0d0 468 plm(2,1)=x 469 plm(2,2)=w 470 plm(3,2)=3.0d0*x*w 471 do m1=1,lx 472 dm=m1-1 473 if(m1 > 1) then 474 plm(m1+1,m1)=x*plm(m1,m1) + 2*dm*w*plm(m1,m1-1) 475 end if 476 if(m1 < lx) then 477 do l1=m1+2,lx+1 478 dl=l1-1 479 plm(l1,m1)=((2*dl-1)*x*plm(l1-1,m1)& 480 & - (dl+dm-1)*plm(l1-2,m1))/(dl-dm) 481 end do 482 end if 483 plm(m1+1,m1+1)=(2*dm+1)*w*plm(m1,m1) 484 end do 485 486 !azimuthal angle phase factors 487 rrs=xx**2 + yy**2 488 if(rrs > tol8) then 489 rr=sqrt(rrs) 490 dc=xx/rr 491 ds=yy/rr 492 else 493 dc=1.0d0 494 ds=0.0d0 495 end if 496 cosa(1)=1.0d0 497 sina(1)=0.0d0 498 do m1=2,lx+1 499 cosa(m1)=dc*cosa(m1-1) - ds*sina(m1-1) 500 sina(m1)=ds*cosa(m1-1) + dc*sina(m1-1) 501 end do 502 503 !combine factors 504 do l1=1,lx+1 505 do m1=2,l1 506 nn=(l1-1)**2 + (l1-1) + (m1-1) + 1 507 nc=(l1-1)**2 + (l1-1) - (m1-1) + 1 508 ! note that we are supressing the so-called Condon-Shortley phase 509 ! ya=sgn(m1)*qlm(l1,m1)*plm(l1,m1) 510 ya=qlm(l1,m1)*plm(l1,m1) 511 yr=ya*cosa(m1) 512 yi=ya*sina(m1) 513 ylm(nc)=sgn(m1)*cmplx(yr,-yi) 514 ylm(nn)=cmplx(yr,yi) 515 end do 516 end do 517 do l1=1,lx+1 518 nn=(l1-1)**2 + (l1-1) + 1 519 ya=qlm(l1,1)*plm(l1,1) 520 ylm(nn)=cmplx(ya,0.d0) 521 end do 522 523 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.
PARENTS
CHILDREN
SOURCE
93 function ylmc(il,im,kcart) 94 95 96 !This section has been created automatically by the script Abilint (TD). 97 !Do not modify the following lines by hand. 98 #undef ABI_FUNC 99 #define ABI_FUNC 'ylmc' 100 !End of the abilint section 101 102 implicit none 103 104 !Arguments ------------------------------------ 105 !scalars 106 integer,intent(in) :: il,im 107 complex(dpc) :: ylmc 108 !arrays 109 real(dp),intent(in) :: kcart(3) 110 111 !Local variables------------------------------- 112 !scalars 113 integer,parameter :: LMAX=3 114 real(dp),parameter :: PPAD=tol8 115 real(dp) :: cosphi,costh,costhreephi,costwophi,r,rxy,sinphi,sinth,sinthreephi,sintwophi 116 !complex(dpc) :: new_ylmc 117 character(len=500) :: msg 118 119 ! ************************************************************************* 120 121 if (ABS(im)>ABS(il)) then 122 write(msg,'(3(a,i0))') 'm is,',im,' however it should be between ',-il,' and ',il 123 MSG_ERROR(msg) 124 end if 125 126 ylmc = czero 127 128 r=SQRT(kcart(1)**2+kcart(2)**2+kcart(3)**2) 129 if (r<PPAD) r=r+PPAD 130 !$if (r<tol10) RETURN 131 132 rxy=SQRT(kcart(1)**2+kcart(2)**2) 133 if (rxy<PPAD)rxy=r+PPAD 134 ! 135 ! Determine theta and phi 136 costh= kcart(3)/r 137 138 #if 1 139 ! old buggy coding 140 sinth= rxy/r 141 cosphi= kcart(1)/rxy 142 sinphi= kcart(2)/rxy 143 #else 144 sinth=sqrt(abs((one-costh)*(one+costh))) ! abs is needed to prevent very small negative arg 145 cosphi=one 146 sinphi=zero 147 if (sinth>tol10) then 148 cosphi=kcart(1)/(r*sinth) 149 sinphi=kcart(2)/(r*sinth) 150 end if 151 #endif 152 153 costwophi= two*cosphi**2 - one 154 sintwophi= two*sinphi*cosphi 155 costhreephi=cosphi*costwophi-sinphi*sintwophi 156 sinthreephi=cosphi*sintwophi+sinphi*costwophi 157 158 select case (il) 159 160 case (0) 161 ylmc= one/SQRT(four_pi) 162 163 case (1) 164 if (ABS(im)==0) then 165 ylmc = SQRT(three/(four_pi))*costh 166 else if (ABS(im)==1) then 167 ylmc = -SQRT(three/(8._dp*pi))*sinth*CMPLX(cosphi,sinphi) 168 else 169 msg='wrong im' 170 MSG_ERROR(msg) 171 end if 172 173 case (2) 174 if (ABS(im)==0) then 175 ylmc = SQRT(5.d0/(16.d0*pi))*(three*costh**2-one) 176 else if (ABS(im)==1) then 177 ylmc = -SQRT(15.d0/(8.d0*pi))*sinth*costh*cmplx(cosphi,sinphi) 178 else if (ABS(im)==2) then 179 ylmc = SQRT(15.d0/(32.d0*pi))*(sinth)**2*CMPLX(costwophi,sintwophi) 180 else 181 msg='wrong im' 182 MSG_ERROR(msg) 183 end if 184 185 case (3) 186 if (ABS(im)==0) then 187 ylmc= SQRT(7.d0/(16.d0*pi))*(5.d0*costh**3 -3.d0*costh) 188 else if (ABS(im)==1) then 189 ylmc= -SQRT(21.d0/(64.d0*pi))*sinth*(5.d0*costh**2-one)*CMPLX(cosphi,sinphi) 190 else if (ABS(im)==2) then 191 ylmc= SQRT(105.d0/(32.d0*pi))*sinth**2*costh*CMPLX(costwophi,sintwophi) 192 else if (ABS(im)==3) then 193 ylmc=-SQRT(35.d0/(64.d0*pi))*sinth**3*CMPLX(costhreephi,sinthreephi) 194 else 195 msg='wrong im' 196 MSG_ERROR(msg) 197 end if 198 199 case default 200 !write(msg,'(a,i6,a,i6)')' The maximum allowed value for l is,',LMAX,' however l=',il 201 !MSG_ERROR(msg) 202 end select 203 ! 204 !=== Treat the case im < 0 === 205 if (im < 0) ylmc=(-one)**(im)*CONJG(ylmc) 206 207 ! FIXME: Use the piece of code below as it works for arbitrary (l,m) 208 ! the implementation above is buggy when the vector is along z! 209 ! 210 #if 0 211 ! Remember the expression of complex spherical harmonics: 212 ! $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}$ 213 new_ylmc = SQRT((2*il+1)*dble_factorial(il-ABS(im))/(dble_factorial(il+ABS(im))*four_pi)) * & 214 & ass_leg_pol(il,ABS(im),costh) * CMPLX(cosphi,sinphi)**ABS(im) 215 if (im<0) new_ylmc=(-one)**(im)*CONJG(new_ylmc) 216 217 if (ABS(new_ylmc-ylmc)>tol6) then 218 !MSG_WARNING("Check new_ylmc") 219 !write(std_out,*)"il,im,new_ylmc, ylmc",il,im,new_ylmc,ylmc 220 !write(std_out,*)"fact",SQRT((2*il+1)*dble_factorial(il-ABS(im))/(dble_factorial(il+ABS(im))*four_pi)) 221 !write(std_out,*)"costh,sinth,ass_leg_pol",costh,sinth,ass_leg_pol(il,ABS(im),costh) 222 !write(std_out,*)"cosphi,sinphi,e^{imphi}",cosphi,sinphi,CMPLX(cosphi,sinphi)**ABS(im) 223 end if 224 ylmc = new_ylmc 225 #endif 226 227 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.
PARENTS
m_vkbr
CHILDREN
wrtout
SOURCE
262 subroutine ylmcd(il,im,kcart,dth,dphi) 263 264 265 !This section has been created automatically by the script Abilint (TD). 266 !Do not modify the following lines by hand. 267 #undef ABI_FUNC 268 #define ABI_FUNC 'ylmcd' 269 !End of the abilint section 270 271 implicit none 272 273 !Arguments ------------------------------------ 274 !scalars 275 integer,intent(in) :: il,im 276 complex(dpc),intent(out) :: dphi,dth 277 !arrays 278 real(dp),intent(in) :: kcart(3) 279 280 !Local variables------------------------------- 281 !scalars 282 integer,parameter :: LMAX=3 283 real(dp),parameter :: PPAD=tol8 284 real(dp) :: cosphi,costh,costhreephi,costwophi,r,rxy,sinphi,sinth,sinthreephi,sintwophi,c 285 character(len=500) :: msg 286 287 ! ************************************************************************* 288 289 if (ABS(im)>ABS(il))then 290 write(msg,'(3(a,i0))')' m is,',im,' however it should be between ',-il,' and ',il 291 MSG_ERROR(msg) 292 end if 293 294 dphi=czero; dth=czero 295 296 r=SQRT(kcart(1)**2+kcart(2)**2+kcart(3)**2) 297 if (r<PPAD) r=r+PPAD 298 !$if (r<tol10) RETURN 299 300 rxy=SQRT(kcart(1)**2+kcart(2)**2) 301 if (rxy<PPAD) rxy=r+PPAD 302 303 ! Determine theta and phi 304 costh= kcart(3)/r 305 #if 1 306 ! old buggy coding 307 sinth= rxy/r 308 cosphi= kcart(1)/rxy 309 sinphi= kcart(2)/rxy 310 #else 311 sinth=sqrt(abs((one-costh)*(one+costh))) ! abs is needed to prevent very small negative arg 312 cosphi=one 313 sinphi=zero 314 if (sinth>tol10) then 315 cosphi=kcart(1)/(r*sinth) 316 sinphi=kcart(2)/(r*sinth) 317 end if 318 #endif 319 320 costwophi= two*cosphi**2 - one 321 sintwophi= two*sinphi*cosphi 322 costhreephi=cosphi*costwophi-sinphi*sintwophi 323 sinthreephi=cosphi*sintwophi+sinphi*costwophi 324 325 select case (il) 326 327 case (0) 328 dth = czero 329 dphi = czero 330 331 case (1) 332 if (ABS(im)==0) then 333 dth= -SQRT(three/(four_pi))*sinth 334 dphi= czero 335 else if (abs(im)==1) then 336 dth= -SQRT(3.d0/(8.d0*pi))*costh*CMPLX(cosphi,sinphi) 337 dphi=-SQRT(3.d0/(8.d0*pi))*sinth*CMPLX(-sinphi,cosphi) 338 end if 339 340 case (2) 341 if (ABS(im)==0) then 342 dth= -SQRT(5.d0/(16.d0*pi))*6.d0*costh*sinth 343 dphi= czero 344 else if (ABS(im)==1) then 345 dth= -SQRT(15.d0/(8.d0*pi))*(costh**2-sinth**2)*CMPLX(cosphi,sinphi) 346 dphi= -SQRT(15.d0/(8.d0*pi))*costh*sinth*(0.d0,1.d0)*CMPLX(cosphi,sinphi) 347 else if (abs(im)==2) then 348 dth = SQRT(15.d0/(32.d0*pi))*2.d0*costh*sinth*CMPLX(costwophi,sintwophi) 349 dphi = SQRT(15.d0/(32.d0*pi))*sinth**2*(0.d0,2.d0)*CMPLX(costwophi,sintwophi) 350 end if 351 352 case (3) 353 if (ABS(im)==0) then 354 dth = SQRT(7.d0/(16*pi))*(-15.d0*costh**2*sinth + 3.d0**sinth) 355 dphi= czero 356 else if (ABS(im)==1) then 357 c = SQRT(21.d0/(64.d0*pi)) 358 dth= -c* (15.d0*costh**3-11.d0*costh)* CMPLX(cosphi,sinphi) 359 dphi=-c*sinth*( 5.d0*costh**2-1 )*(0.d0,1.d0)*CMPLX(cosphi,sinphi) 360 else if (ABS(im)==2) then 361 c = SQRT(105.d0/(32.d0*pi)) 362 dth =c*(2.d0*sinth*costh**2-sinth**3) *CMPLX(costwophi,sintwophi) 363 dphi=c*(2.d0*sinth**2*costh)*(0.d0,1.d0)*CMPLX(costwophi,sintwophi) 364 else if (abs(im)==3) then 365 dth =-SQRT(35.d0/(64.d0*pi))*3.d0*sinth**2*costh*CMPLX(costhreephi,sinthreephi) 366 dphi=-SQRT(35.d0/(64.d0*pi))*sinth**3*(0.d0,3.d0)*CMPLX(costhreephi,sinthreephi) 367 end if 368 369 case default 370 write(msg,'(2(a,i0))')' The maximum allowed value for l is,',LMAX,' however, l=',il 371 MSG_ERROR(msg) 372 end select 373 ! 374 !=== Treat the case im < 0 === 375 if (im<0) then 376 dth = (-one)**(im)*CONJG(dth) 377 dphi= (-one)**(im)*CONJG(dphi) 378 end if 379 380 end subroutine ylmcd
m_paw_sphharm/ys [ Functions ]
[ Top ] [ m_paw_sphharm ] [ Functions ]
NAME
ys
FUNCTION
Computes the matrix element <Yl'm'|Slm>
INPUTS
integer :: l',m',l,m
OUTPUT
complex(dpc) :: ys_val
NOTES
Ylm is the standard complex-valued spherical harmonic, Slm is the real spherical harmonic used througout abinit. <Yl'm'|Slm> is their overlap.
PARENTS
m_epjdos,m_paw_sphharm
CHILDREN
wrtout
SOURCE
789 subroutine ys(lp,mp,ll,mm,ys_val) 790 791 792 !This section has been created automatically by the script Abilint (TD). 793 !Do not modify the following lines by hand. 794 #undef ABI_FUNC 795 #define ABI_FUNC 'ys' 796 !End of the abilint section 797 798 implicit none 799 800 !Arguments --------------------------------------------- 801 !scalars 802 integer,intent(in) :: ll,lp,mm,mp 803 complex(dpc),intent(out) :: ys_val 804 805 !Local variables --------------------------------------- 806 !scalars 807 complex(dpc) :: dmpmm,dmpmmm,m1mm 808 809 ! ********************************************************************* 810 811 812 ys_val = czero 813 814 if(lp==ll .AND. (mp==mm .OR. mp==-mm) ) then 815 ! (-1)**mm 816 m1mm=cone; if(abs(mod(mm,2))==1) m1mm=-m1mm 817 818 ! delta(mp,mm) 819 dmpmm=czero; if(mp==mm) dmpmm=cone 820 821 ! delta(mp,-mm) 822 dmpmmm=czero; if(mp==-mm) dmpmmm=cone 823 824 select case (mm) 825 case (0) ! case for S_l0 826 ys_val = dmpmm 827 case (:-1) ! case for S_lm with m < 0 828 ys_val = -(zero,one)*m1mm*sqrthalf*(dmpmmm-m1mm*dmpmm) 829 case (1:) ! case for S_lm with m > 0 830 ys_val = m1mm*sqrthalf*(dmpmm+m1mm*dmpmmm) 831 end select 832 833 end if 834 835 end subroutine ys