TABLE OF CONTENTS
- ABINIT/m_pawang
- m_pawang/gauleg
- m_pawang/gaunt
- m_pawang/initang
- m_pawang/pawang_free
- m_pawang/pawang_init
- m_pawang/pawang_lsylm
- m_pawang/pawang_type
- m_pawang/perms
- m_pawang/realgaunt
- m_pawang/rfactorial
ABINIT/m_pawang [ Modules ]
NAME
m_pawang
FUNCTION
This module contains the definition of the pawang_type structured datatype, as well as related functions and methods. pawang_type variables define ANGular mesh discretization of PAW augmentation regions and related data.
COPYRIGHT
Copyright (C) 2013-2018 ABINIT group (MT, FJ, BA) 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
24 #include "libpaw.h" 25 26 MODULE m_pawang 27 28 USE_DEFS 29 USE_MSG_HANDLING 30 USE_MPI_WRAPPERS 31 USE_MEMORY_PROFILING 32 33 use m_paw_sphharm, only : initylmr, mat_mlms2jmj, mat_slm2ylm 34 35 implicit none 36 37 private 38 39 !public procedures. 40 public :: pawang_init ! Constructor 41 public :: pawang_free ! Free memory 42 public :: pawang_lsylm ! Compute the LS operator in the real spherical harmonics basis 43 public :: initang ! Initialize angular mesh for PAW calculations 44 45 ! MGPAW: Private? 46 public :: realgaunt ! compute "real Gaunt coefficients" with "real spherical harmonics" 47 public :: gaunt ! Gaunt coeffients for complex Yml 48 public :: gauleg 49 public :: mat_mlms2jmj 50 public :: mat_slm2ylm
m_pawang/gauleg [ Functions ]
[ Top ] [ m_pawang ] [ Functions ]
NAME
gauleg
FUNCTION
Compute the coefficients (supports and weights) for Gauss-Legendre integration
INPUTS
xmin=lower bound of integration xmax=upper bound of integration n=order of integration
OUTPUT
x(n)=array of support points weights(n)=array of integration weights
PARENTS
m_pawang
CHILDREN
SOURCE
983 subroutine gauleg(xmin,xmax,x,weights,n) 984 985 986 !This section has been created automatically by the script Abilint (TD). 987 !Do not modify the following lines by hand. 988 #undef ABI_FUNC 989 #define ABI_FUNC 'gauleg' 990 !End of the abilint section 991 992 implicit none 993 994 !Arguments --------------------------------------------- 995 !scalars 996 integer,intent(in) :: n 997 real(dp),intent(in) :: xmax,xmin 998 !arrays 999 real(dp),intent(out) :: weights(n),x(n) 1000 1001 !Local variables ------------------------------ 1002 !scalars 1003 integer :: ii,jj 1004 real(dp),parameter :: tol=1.d-13 1005 real(dp) :: p1,p2,p3,pi,xl,pp,xmean,z,z1 1006 !arrays 1007 1008 !************************************************************************ 1009 1010 pi=4._dp*atan(1._dp) 1011 xl=(xmax-xmin)*0.5_dp 1012 xmean=(xmax+xmin)*0.5_dp 1013 1014 do ii=1,(n+1)/2 1015 z=cos(pi*(ii-0.25_dp)/(n+0.5_dp)) 1016 do 1017 p1=1._dp 1018 p2=0._dp 1019 do jj=1,n 1020 p3=p2 1021 p2=p1 1022 p1=((2._dp*jj-1._dp)*z*p2-(jj-1._dp)*p3)/jj 1023 end do 1024 pp=n*(p2-z*p1)/(1._dp-z**2) 1025 z1=z 1026 z=z1-p1/pp 1027 if(abs(z-z1) < tol) exit 1028 end do 1029 x(ii)=xmean-xl*z 1030 x(n+1-ii)=xmean+xl*z 1031 weights(ii)=2._dp*xl/((1._dp-z**2)*pp**2) 1032 weights(n+1-ii)=weights(ii) 1033 end do 1034 1035 end subroutine gauleg
m_pawang/gaunt [ Functions ]
[ Top ] [ m_pawang ] [ 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
CHILDREN
SOURCE
856 function gaunt(ll,mm,l1,m1,l2,m2) 857 858 859 !This section has been created automatically by the script Abilint (TD). 860 !Do not modify the following lines by hand. 861 #undef ABI_FUNC 862 #define ABI_FUNC 'gaunt' 863 !End of the abilint section 864 865 implicit none 866 867 !Arguments --------------------------------------------- 868 !scalars 869 integer,intent(in) :: l1,l2,ll,m1,m2,mm 870 real(dp) :: gaunt 871 872 !Local variables ------------------------------ 873 !scalars 874 integer :: i1,i2,j1,j1half,j2,j2half,j3,j3half,j_half,jj,k1,k2,n1,n2 875 real(dp) :: argument,sign,sum,xx,yy 876 logical :: ok 877 878 !************************************************************************ 879 880 gaunt=zero;sum=zero;ok =.true. 881 882 if((-m1-mm+m2) /= 0) ok = .false. 883 if(abs(m1) > l1) ok = .false. 884 if(abs(mm) > ll) ok = .false. 885 if(abs(m2) > l2) ok = .false. 886 887 jj = l1 + ll + l2 888 if (mod(jj,2)/=0) ok = .false. 889 j1 = jj-2*l2 890 j2 = jj-2*ll 891 j3 = jj-2*l1 892 893 if (j1<0 .or. j2<0 .or. j3<0) ok = .false. 894 895 if (ok) then 896 897 xx = (2 * l1 + 1) * (2 * ll + 1) * (2 * l2 + 1) 898 899 j1half = j1/2 900 j2half = j2/2 901 j3half = j3/2 902 j_half = jj/2 903 904 gaunt = (-1)**j1half * sqrt(xx) 905 gaunt = gaunt * rfactorial(j2)*rfactorial(j3)/rfactorial(jj+1) 906 gaunt = gaunt * rfactorial(j_half)/(rfactorial(j1half)& 907 & * rfactorial(j2half)*rfactorial(j3half)) 908 909 yy = rfactorial(l2 + m2) * rfactorial(l2 - m2) 910 911 if (mm>=0) then 912 yy = yy * perms(ll+mm,2*mm) 913 else 914 yy = yy / perms(ll-mm,-2*mm) 915 end if 916 917 if (m1>=0) then 918 yy = yy / perms(l1+m1,2*m1) 919 else 920 yy = yy * perms(l1-m1,-2*m1) 921 end if 922 923 gaunt = gaunt * sqrt(yy) 924 925 i1 = l2 - ll - m1 926 i2 = l2 - l1 + mm 927 k1 = -min(0, i1, i2) 928 n1 = l1 + m1 929 n2 = ll - mm 930 k2 = min(j1, n1, n2) 931 932 sign = 1._dp 933 if(k1>0) sign = (-1._dp)**k1 934 935 argument = sign * perms(n1,k1)/rfactorial(k1) 936 argument = argument * perms(n2,k1)/rfactorial(i1 + k1) 937 argument = argument * perms(j1,k1)/rfactorial(i2 + k1) 938 sum = sum + argument 939 940 sign = -sign 941 k1 = k1 + 1 942 do while(k1 <= k2) 943 argument = sign * perms(n1, k1)/rfactorial(k1) 944 argument = argument * perms(n2, k1)/rfactorial(i1 + k1) 945 argument = argument * perms(j1, k1)/rfactorial(i2 + k1) 946 sum = sum + argument 947 sign = -sign 948 k1 = k1 + 1 949 end do 950 951 end if 952 953 gaunt = gaunt * sum 954 955 end function gaunt
m_pawang/initang [ Functions ]
[ Top ] [ m_pawang ] [ Functions ]
NAME
initang
FUNCTION
Initialize angular mesh for PAW calculations
INPUTS
pawang <type(pawang_type)>=paw angular mesh and related data pawang%angl_size - Total number of sample points in the angular mesh pawang%ntheta - Number of sample points in the theta dir pawang%nphi - Number of sample points in the phi dir
OUTPUT
pawang <type(pawang_type)>=paw angular mesh and related data pawang%anginit - (3 x angl_size) array, the ntheta*nphi dimensional arrays ax, ay, and az pawang%angwgth - (angl_size) array, the weight factor of the point (ax, ay, az)
PARENTS
m_pawang
CHILDREN
SOURCE
589 subroutine initang(pawang) 590 591 592 !This section has been created automatically by the script Abilint (TD). 593 !Do not modify the following lines by hand. 594 #undef ABI_FUNC 595 #define ABI_FUNC 'initang' 596 !End of the abilint section 597 598 implicit none 599 600 !Arguments ------------------------------------ 601 !scalars 602 type(pawang_type),intent(inout) :: pawang 603 604 !Local variables------------------------------- 605 !scalars 606 integer :: ip,it,npoints 607 real(dp) :: ang,con,cos_phi,cos_theta,sin_phi,sin_theta 608 character(len=500) :: msg 609 !arrays 610 real(dp) :: th(pawang%ntheta),wth(pawang%ntheta) 611 612 ! *********************************************************************** 613 614 if (pawang%angl_size==0) return 615 616 !Initializations 617 npoints=0 618 con=two_pi / pawang%nphi 619 call gauleg(-one,one,th,wth,pawang%ntheta) 620 621 !We now open two nested do-loops. The first loops through the number 622 !of theta angles, the second through the number of phi angles (?). 623 !The two together initialize anginit. 624 625 do it = 1, pawang%ntheta 626 627 cos_theta = th(it) 628 sin_theta = sqrt(one - cos_theta*cos_theta) 629 630 do ip = 1, pawang%nphi 631 632 ang = con * (ip-1) 633 cos_phi = cos(ang) 634 sin_phi = sin(ang) 635 636 npoints = npoints + 1 637 638 pawang%anginit(1, npoints) = sin_theta * cos_phi 639 pawang%anginit(2, npoints) = sin_theta * sin_phi 640 pawang%anginit(3, npoints) = cos_theta 641 642 ! Normalization required 643 pawang%angwgth(npoints) = wth(it) / (2 * pawang%nphi) 644 645 end do 646 end do 647 648 !The following is an error statement that will be generated 649 !if npoints exceeds nang... 650 if (npoints > pawang%angl_size) then 651 write(msg, '(a,i4,a,a,i4)' ) & 652 & ' anginit%npoints =',npoints,ch10,& 653 & ' angl_size =',pawang%angl_size 654 MSG_BUG(msg) 655 end if 656 657 end subroutine initang
m_pawang/pawang_free [ Functions ]
[ Top ] [ m_pawang ] [ Functions ]
NAME
pawang_free
FUNCTION
Free all dynamic memory and reset all flags stored in a pawang datastructure
SIDE EFFECTS
Pawang <type(pawang_type)>=ANGular mesh discretization and related data
PARENTS
dfpt_looppert,driver,pawinit
CHILDREN
SOURCE
303 subroutine pawang_free(Pawang) 304 305 306 !This section has been created automatically by the script Abilint (TD). 307 !Do not modify the following lines by hand. 308 #undef ABI_FUNC 309 #define ABI_FUNC 'pawang_free' 310 !End of the abilint section 311 312 implicit none 313 314 !Arguments ------------------------------------ 315 !scalars 316 type(Pawang_type),intent(inout) :: Pawang 317 318 ! ************************************************************************* 319 320 !@Pawang_type 321 if (allocated(pawang%angwgth)) then 322 LIBPAW_DEALLOCATE(pawang%angwgth) 323 end if 324 if (allocated(pawang%anginit)) then 325 LIBPAW_DEALLOCATE(pawang%anginit) 326 end if 327 if (allocated(pawang%zarot)) then 328 LIBPAW_DEALLOCATE(pawang%zarot) 329 end if 330 if (allocated(pawang%gntselect)) then 331 LIBPAW_DEALLOCATE(pawang%gntselect) 332 end if 333 if (allocated(pawang%realgnt)) then 334 LIBPAW_DEALLOCATE(pawang%realgnt) 335 end if 336 if (allocated(pawang%ylmr)) then 337 LIBPAW_DEALLOCATE(pawang%ylmr) 338 end if 339 if (allocated(pawang%ylmrgr)) then 340 LIBPAW_DEALLOCATE(pawang%ylmrgr) 341 end if 342 if (allocated(pawang%ls_ylm)) then 343 LIBPAW_DEALLOCATE(pawang%ls_ylm) 344 end if 345 346 pawang%angl_size =0 347 pawang%ylm_size =0 348 pawang%use_ls_ylm=0 349 pawang%l_max=-1 350 pawang%l_size_max=-1 351 pawang%gnt_option=-1 352 pawang%ngnt=0 353 354 end subroutine pawang_free
m_pawang/pawang_init [ Functions ]
[ Top ] [ m_pawang ] [ Functions ]
NAME
pawang_init
FUNCTION
Initialize a pawang datastructure
INPUTS
gnt_option=flag activated if pawang%gntselect and pawang%realgnt have to be allocated also determine the size of these pointers lmax=maximum value of angular momentum l nphi,ntheta=dimensions of paw angular mesh nsym=number of symetries pawxcdev=choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments) use_ls_ylm=flag activated if pawang%ls_ylm has to be allocated use_ylm=flag activated if pawang%ylmr and pawang%ylmrgr have to be allocated xclevel=XC functional level (1=LDA, 2=GGA)
OUTPUT
Pawang <type(pawang_type)>=ANGular mesh discretization and related data
PARENTS
dfpt_looppert,pawinit
CHILDREN
SOURCE
183 subroutine pawang_init(Pawang,gnt_option,lmax,nphi,nsym,ntheta,pawxcdev,use_ls_ylm,use_ylm,xclevel) 184 185 186 !This section has been created automatically by the script Abilint (TD). 187 !Do not modify the following lines by hand. 188 #undef ABI_FUNC 189 #define ABI_FUNC 'pawang_init' 190 !End of the abilint section 191 192 implicit none 193 194 !Arguments ------------------------------------ 195 !scalars 196 integer,intent(in) :: gnt_option,lmax,nphi,nsym,ntheta 197 integer,intent(in) :: pawxcdev,use_ls_ylm,use_ylm,xclevel 198 type(Pawang_type),intent(inout) :: Pawang 199 200 !Local variables------------------------------- 201 !scalars 202 integer :: ll,sz1,sz2,sz3 203 !arrays 204 real(dp),allocatable :: rgnt_tmp(:) 205 206 ! ************************************************************************* 207 208 !@Pawang_type 209 210 Pawang%l_max=lmax+1 211 Pawang%l_size_max=2*Pawang%l_max-1 212 Pawang%nsym=nsym 213 214 if (pawxcdev==0) then 215 Pawang%nphi=nphi 216 Pawang%ntheta=ntheta 217 Pawang%angl_size=ntheta*nphi 218 else 219 Pawang%nphi=0 220 Pawang%ntheta=0 221 Pawang%angl_size=0 222 end if 223 224 if (Pawang%angl_size>0) then 225 LIBPAW_ALLOCATE(Pawang%anginit,(3,Pawang%angl_size)) 226 LIBPAW_ALLOCATE(Pawang%angwgth,(Pawang%angl_size)) 227 call initang(pawang) 228 end if 229 230 if (use_ylm>0.and.Pawang%angl_size>0) then 231 if (xclevel>=0) ll=Pawang%l_size_max 232 if (xclevel>=2) ll=ll+1 233 Pawang%ylm_size=ll**2 234 LIBPAW_ALLOCATE(Pawang%ylmr,(Pawang%ylm_size,Pawang%angl_size)) 235 if (xclevel==2) then 236 LIBPAW_ALLOCATE(Pawang%ylmrgr,(3,Pawang%ylm_size,Pawang%angl_size)) 237 call initylmr(ll,0,pawang%angl_size,pawang%angwgth,2,pawang%anginit,pawang%ylmr,& 238 & ylmr_gr=pawang%ylmrgr) 239 else 240 call initylmr(ll,0,pawang%angl_size,pawang%angwgth,1,pawang%anginit,pawang%ylmr) 241 end if 242 else 243 Pawang%ylm_size=0 244 end if 245 246 Pawang%gnt_option=gnt_option 247 if (Pawang%gnt_option==1.or.Pawang%gnt_option==2) then 248 if (Pawang%gnt_option==1) then 249 sz1=(2*Pawang%l_max-1)**2*(Pawang%l_max)**4 250 sz2=(Pawang%l_size_max)**2 251 sz3=(Pawang%l_max**2)*(Pawang%l_max**2+1)/2 252 LIBPAW_ALLOCATE(rgnt_tmp,(sz1)) 253 LIBPAW_ALLOCATE(pawang%gntselect,(sz2,sz3)) 254 call realgaunt(Pawang%l_max,Pawang%ngnt,Pawang%gntselect,rgnt_tmp) 255 else if (Pawang%gnt_option==2) then 256 sz1=(4*Pawang%l_max-3)**2*(2*Pawang%l_max-1)**4 257 sz2=(2*Pawang%l_size_max-1)**2 258 sz3=((2*Pawang%l_max-1)**2)*((2*Pawang%l_max-1)**2+1)/2 259 LIBPAW_ALLOCATE(rgnt_tmp,(sz1)) 260 LIBPAW_ALLOCATE(pawang%gntselect,(sz2,sz3)) 261 call realgaunt(2*Pawang%l_max-1,Pawang%ngnt,Pawang%gntselect,rgnt_tmp) 262 end if 263 if (allocated(pawang%realgnt)) then 264 LIBPAW_DEALLOCATE(pawang%realgnt) 265 end if 266 LIBPAW_ALLOCATE(Pawang%realgnt,(Pawang%ngnt)) 267 Pawang%realgnt(1:Pawang%ngnt)=rgnt_tmp(1:Pawang%ngnt) 268 LIBPAW_DEALLOCATE(rgnt_tmp) 269 end if 270 271 Pawang%use_ls_ylm=use_ls_ylm 272 if (use_ls_ylm>0) then 273 LIBPAW_ALLOCATE(pawang%ls_ylm,(2,Pawang%l_max**2*(Pawang%l_max**2+1)/2,2)) 274 call pawang_lsylm(pawang) 275 end if 276 277 if (nsym>0) then 278 LIBPAW_ALLOCATE(Pawang%zarot,(Pawang%l_size_max,Pawang%l_size_max,Pawang%l_max,nsym)) 279 end if 280 281 end subroutine pawang_init
m_pawang/pawang_lsylm [ Functions ]
[ Top ] [ m_pawang ] [ Functions ]
NAME
pawang_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
pawang <type(pawang_type)>=paw angular mesh and related data
OUTPUT
pawang%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>
PARENTS
m_pawang
CHILDREN
SOURCE
389 subroutine pawang_lsylm(pawang) 390 391 392 !This section has been created automatically by the script Abilint (TD). 393 !Do not modify the following lines by hand. 394 #undef ABI_FUNC 395 #define ABI_FUNC 'pawang_lsylm' 396 !End of the abilint section 397 398 implicit none 399 400 !Arguments --------------------------------------------- 401 !scalars 402 type(pawang_type),intent(inout) :: pawang 403 404 !Local variables --------------------------------------- 405 !scalars 406 integer :: ii,ilm,im,j0lm,jj,jlm,jm,klm,l_max,ll,lm0,mm,ispden 407 real(dp),parameter :: invsqrt2=one/sqrt2 408 real(dp) :: onem 409 character(len=500) :: msg 410 logical,parameter :: tso=.false. ! use true to Test Spin Orbit and 411 ! write the matrix of L.S in different basis 412 !arrays 413 complex(dpc) :: tmp(2) 414 complex(dpc),allocatable :: ls_cplx(:,:,:),slm2ylm(:,:) 415 complex(dpc),allocatable :: mat_inp_c(:,:,:),mat_out_c(:,:,:) 416 complex(dpc),allocatable :: mat_ls_ylm(:,:,:),mat_jmj(:,:) 417 character(len=9),parameter :: dspin2(2)=(/"up-up ","up-dn "/) 418 character(len=9),parameter :: dspin6(6)=(/"dn ","up ","dn-dn ","up-up ","dn-up ","up-dn "/) 419 character(len=9),parameter :: dspinm(6)=(/"dn ","up ","n ","mx ","my ","mz "/) 420 421 ! ************************************************************************* 422 423 if (pawang%use_ls_ylm==0) then 424 msg=' ls_ylm pointer is not allocated !' 425 MSG_BUG(msg) 426 end if 427 428 !Initialization 429 pawang%ls_ylm=zero 430 l_max=pawang%l_max-1 431 432 !Nothing to do if lmax=0 433 if (l_max<=0) return 434 435 !Loop on l quantum number 436 do ll=1,l_max 437 438 ! Transformation matrixes: real->complex spherical harmonics 439 LIBPAW_ALLOCATE(slm2ylm,(2*ll+1,2*ll+1)) 440 slm2ylm=czero 441 do im=1,2*ll+1 442 mm=im-ll-1;jm=-mm+ll+1 443 onem=dble((-1)**mm) 444 if (mm> 0) then 445 slm2ylm(im,im)= cmplx(onem*invsqrt2,zero,kind=dp) 446 slm2ylm(jm,im)= cmplx(invsqrt2, zero,kind=dp) 447 end if 448 if (mm==0) then 449 slm2ylm(im,im)=cone 450 end if 451 if (mm< 0) then 452 slm2ylm(im,im)= cmplx(zero, invsqrt2,kind=dp) 453 slm2ylm(jm,im)=-cmplx(zero,onem*invsqrt2,kind=dp) 454 end if 455 end do 456 457 ! Compute <sigma, Y_lm1|L.S|Y_lm2, sigma_prime> (Y_lm=complex spherical harmonics) 458 ! 1= <up|L.S|up> ; 2= <up|L.S|dn> 459 LIBPAW_ALLOCATE(ls_cplx,(2*ll+1,2*ll+1,2)) 460 ls_cplx=czero 461 if(tso) then 462 LIBPAW_ALLOCATE(mat_ls_ylm,(2*ll+1,2*ll+1,4)) 463 if(tso) mat_ls_ylm=czero 464 end if 465 if(tso) then 466 LIBPAW_ALLOCATE(mat_jmj,(2*(2*ll+1),2*(2*ll+1))) 467 if(tso) mat_jmj=czero 468 end if 469 do im=1,2*ll+1 470 mm=im-ll-1 471 ls_cplx(im,im,1)=half*mm 472 if(tso) mat_ls_ylm(im,im,1)=-half*mm ! dn dn 473 if(tso) mat_ls_ylm(im,im,2)=half*mm ! up up 474 if ((mm+1)<= ll) then 475 ls_cplx(im,im+1,2)=half*sqrt(real((ll-mm)*(ll+mm+1),kind=dp)) 476 if(tso) mat_ls_ylm(im,im+1,4)=half*sqrt(real((ll-mm)*(ll+mm+1),kind=dp)) ! up dn 477 if(tso) mat_ls_ylm(im+1,im,3)=half*sqrt(real((ll-mm)*(ll+mm+1),kind=dp)) ! dn up 478 end if 479 if ((mm-1)>=-ll) then 480 ls_cplx(im-1,im,2)=half*sqrt(real((ll+mm)*(ll-mm+1),kind=dp)) 481 if(tso) mat_ls_ylm(im-1,im,4)=half*sqrt(real((ll+mm)*(ll-mm+1),kind=dp)) ! up dn 482 if(tso) mat_ls_ylm(im,im-1,3)=half*sqrt(real((ll+mm)*(ll-mm+1),kind=dp)) ! dn up 483 end if 484 end do 485 486 ! test : print LS in J,M_J basis 487 if(tso) then 488 do ispden=1,4 489 write(msg,'(3a)') ch10,"value of LS in the Ylm basis for " ,trim(dspin6(ispden+2*(4/4))) 490 call wrtout(std_out,msg,'COLL') 491 do im=1,ll*2+1 492 write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') (mat_ls_ylm(im,jm,ispden),jm=1,ll*2+1) 493 call wrtout(std_out,msg,'COLL') 494 end do 495 end do 496 call mat_mlms2jmj(ll,mat_ls_ylm,mat_jmj,4,1,2,3,std_out,'COLL') ! optspin=2 : dn spin are first 497 end if 498 499 ! Compute <sigma, S_lm1|L.S|S_lm2, sigma_prime> (S_lm=real spherical harmonics) 500 ! 1= <up|L.S|up> ; 2= <up|L.S|dn> 501 if(tso) then 502 LIBPAW_ALLOCATE(mat_inp_c,(2*ll+1,2*ll+1,4)) 503 LIBPAW_ALLOCATE(mat_out_c,(2*ll+1,2*ll+1,4)) 504 end if 505 lm0=ll**2 506 do jm=1,2*ll+1 507 jlm=lm0+jm;j0lm=jlm*(jlm-1)/2 508 do im=1,jm 509 ilm=lm0+im;klm=j0lm+ilm 510 tmp(:)=czero 511 do ii=1,2*ll+1 512 do jj=1,2*ll+1 513 tmp(:)=tmp(:)+ls_cplx(ii,jj,:)*CONJG(slm2ylm(ii,im))*slm2ylm(jj,jm) 514 end do 515 end do 516 pawang%ls_ylm(1,klm,:)=REAL(tmp(:),kind=dp) 517 pawang%ls_ylm(2,klm,:)=AIMAG(tmp(:)) 518 end do 519 end do 520 521 ! Test: print LS in Slm basis 522 if(tso) then 523 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 524 do ispden=1,4 525 write(msg,'(3a)') ch10,"value of LS in the Slm basis for " ,trim(dspin6(ispden+2*(4/4))) 526 call wrtout(std_out,msg,'COLL') 527 do im=1,ll*2+1 528 write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') (mat_inp_c(im,jm,ispden),jm=1,ll*2+1) 529 call wrtout(std_out,msg,'COLL') 530 end do 531 end do 532 ! change into n,m basis 533 mat_ls_ylm(:,:,1)=(mat_inp_c(:,:,1)+mat_inp_c(:,:,2)) 534 mat_ls_ylm(:,:,2)=(mat_inp_c(:,:,3)+mat_inp_c(:,:,4)) 535 mat_ls_ylm(:,:,3)=-cmplx(0.d0,1.d0)*(mat_inp_c(:,:,4)-mat_inp_c(:,:,3)) 536 mat_ls_ylm(:,:,4)=(mat_inp_c(:,:,1)-mat_inp_c(:,:,2)) 537 do ispden=1,4 538 write(msg,'(3a)') ch10,"value of LS in the Slm basis for " ,trim(dspinm(ispden+2*(4/4))) 539 call wrtout(std_out,msg,'COLL') 540 do im=1,ll*2+1 541 write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') (mat_ls_ylm(im,jm,ispden),jm=1,ll*2+1) 542 call wrtout(std_out,msg,'COLL') 543 end do 544 end do 545 LIBPAW_DEALLOCATE(mat_inp_c) 546 LIBPAW_DEALLOCATE(mat_ls_ylm) 547 LIBPAW_DEALLOCATE(mat_jmj) 548 LIBPAW_DEALLOCATE(mat_out_c) 549 end if ! tso 550 551 LIBPAW_DEALLOCATE(ls_cplx) 552 LIBPAW_DEALLOCATE(slm2ylm) 553 554 ! End loop on l 555 end do 556 557 end subroutine pawang_lsylm
m_pawang/pawang_type [ Types ]
[ Top ] [ m_pawang ] [ Types ]
NAME
pawang_type
FUNCTION
For PAW: ANGular mesh discretization of PAW augmentation regions and related data.
SOURCE
66 type,public :: pawang_type 67 68 !Integer scalars 69 70 integer :: angl_size=0 71 ! Dimension of paw angular mesh 72 ! angl_size=ntheta*nphi 73 74 integer :: l_max=-1 75 ! Maximum value of angular momentum l+1 76 77 integer :: l_size_max=-1 78 ! Maximum value of angular momentum +1 79 ! leading to non-zero Gaunt coefficients 80 ! l_size_max = 2*l_max-1 81 82 integer :: ngnt=0 83 ! Number of non zero Gaunt coefficients 84 85 integer :: ntheta, nphi 86 ! Dimensions of paw angular mesh 87 88 integer :: nsym 89 ! Number of symmetry elements in space group 90 91 integer :: gnt_option=-1 92 ! Option for Gaunt coefficients: 93 ! gnt_option==0, Gaunt coeffs are not computed (and not allocated) 94 ! gnt_option==1, Gaunt coeffs are computed up to 2*l_size_max-1 95 ! gnt_option==2, Gaunt coeffs are computed up to l_size_max 96 97 integer :: use_ls_ylm=0 98 ! Flag: use_ls_ylm=1 if pawang%ls_ylm is allocated 99 100 integer :: ylm_size=0 101 ! Size of ylmr/ylmrgr arrays 102 103 !Integer arrays 104 105 integer, allocatable :: gntselect(:,:) 106 ! gntselect(l_size_max**2,l_max**2*(l_max**2+1)/2) 107 ! Selection rules for Gaunt coefficients stored as (LM,ij) where ij is in packed form. 108 ! (if gntselect>0, Gaunt coeff. is non-zero) 109 110 !Real (real(dp)) arrays 111 112 real(dp), allocatable :: anginit(:,:) 113 ! anginit(3,angl_size) 114 ! For each point of the angular mesh, gives the coordinates 115 ! of the corresponding point on an unitary sphere 116 ! Not used in present version (5.3) 117 118 real(dp), allocatable :: angwgth(:) 119 ! angwgth(angl_size) 120 ! For each point of the angular mesh, gives the weight 121 ! of the corresponding point on an unitary sphere 122 123 real(dp), allocatable :: ls_ylm(:,:,:) 124 ! ls_ylm(2,l_max**2*(l_max**2+1)/2,2) 125 ! LS operator in the real spherical harmonics basis 126 ! ls_ylm(ilm1m2,ispin)= <sigma, y_lm1| LS |y_lm2, sigma_prime> 127 128 real(dp), allocatable :: realgnt(:) 129 ! realgnt(ngnt) 130 ! Non zero real Gaunt coefficients 131 132 real(dp), allocatable :: ylmr(:,:) 133 ! ylmr(ylm_size,angl_size) 134 ! Real Ylm calculated in real space 135 136 real(dp), allocatable :: ylmrgr(:,:,:) 137 ! ylmrgr(1:3,ylm_size,angl_size) 138 ! First gradients of real Ylm calculated in real space (cart. coordinates) 139 140 real(dp), allocatable :: zarot(:,:,:,:) 141 ! zarot(l_size_max,l_size_max,l_max,nsym) 142 ! Coeffs of the transformation of real spherical 143 ! harmonics under the symmetry operations symrec. 144 145 end type pawang_type
m_pawang/perms [ Functions ]
[ Top ] [ m_pawang ] [ 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)!
PARENTS
CHILDREN
SOURCE
1113 function perms(nn,kk) 1114 1115 1116 !This section has been created automatically by the script Abilint (TD). 1117 !Do not modify the following lines by hand. 1118 #undef ABI_FUNC 1119 #define ABI_FUNC 'perms' 1120 !End of the abilint section 1121 1122 implicit none 1123 1124 !Arguments --------------------------------------------- 1125 !scalars 1126 integer,intent(in) :: kk,nn 1127 real(dp) :: perms 1128 1129 !Local variables --------------------------------------- 1130 !scalars 1131 integer :: ii 1132 real(dp) :: pp 1133 1134 ! ********************************************************************* 1135 1136 if (nn>=0.and.nn>=kk) then 1137 pp=1._dp 1138 do ii=nn-kk+1,nn 1139 pp=pp*ii 1140 end do 1141 else 1142 pp=0._dp 1143 end if 1144 1145 perms=pp 1146 1147 end function perms
m_pawang/realgaunt [ Functions ]
[ Top ] [ m_pawang ] [ Functions ]
NAME
realgaunt
FUNCTION
This routine compute "real Gaunt coefficients", i.e. gaunt coefficients according to "real spherical harmonics"
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
PARENTS
m_paw_slater,m_pawang,m_pawpwij
CHILDREN
SOURCE
690 subroutine realgaunt(l_max,ngnt,gntselect,realgnt) 691 692 693 !This section has been created automatically by the script Abilint (TD). 694 !Do not modify the following lines by hand. 695 #undef ABI_FUNC 696 #define ABI_FUNC 'realgaunt' 697 !End of the abilint section 698 699 implicit none 700 701 !Arguments --------------------------------------------- 702 !scalars 703 integer,intent(in) :: l_max 704 integer,intent(out) :: ngnt 705 !arrays 706 integer,intent(out) :: gntselect((2*l_max-1)**2,l_max**2*(l_max**2+1)/2) 707 real(dp),intent(out) :: realgnt((2*l_max-1)**2*(l_max)**4) 708 709 !Local variables ------------------------------ 710 !scalars 711 integer :: ilm1,ilm2,ilmp1,k0lm1,klm1,l1,l2,ll,lp1,m1,m2,mm,mm1,mm2,mm3,mp1 712 real(dp) :: c11,c12,c21,c22,c31,c32,fact,realgnt_tmp 713 !arrays 714 integer,allocatable :: ssgn(:) 715 type(coeff3_type), allocatable :: coeff(:) 716 717 !************************************************************************ 718 719 ! Initialize output arrays with zeros. 720 gntselect = 0; realgnt = zero 721 722 !Compute matrix cc where Sl=cc*Yl (Sl=real sph. harm.) 723 !------------------------------------------------ 724 LIBPAW_DATATYPE_ALLOCATE(coeff,(4*l_max-3)) 725 do ll=1,4*l_max-3 726 LIBPAW_ALLOCATE(coeff(ll)%value,(2,2*ll-1,2*ll-1)) 727 coeff(ll)%value(:,:,:)=zero 728 coeff(ll)%value(1,ll,ll)=one 729 do mm=1,ll-1 730 coeff(ll)%value(1,ll+mm,ll+mm)= (-1._dp)**mm/sqrt(2._dp) 731 coeff(ll)%value(1,ll-mm,ll+mm)= ( 1._dp) /sqrt(2._dp) 732 coeff(ll)%value(2,ll+mm,ll-mm)=-(-1._dp)**mm/sqrt(2._dp) 733 coeff(ll)%value(2,ll-mm,ll-mm)= ( 1._dp) /sqrt(2._dp) 734 end do 735 end do 736 737 LIBPAW_ALLOCATE(ssgn,(l_max**2)) 738 ssgn(:)=1 739 if (l_max>0) then 740 do l1=1,l_max-1 741 ilm1=1+l1**2+l1 742 do m1=-l1,-1 743 ssgn(ilm1+m1)=-1 744 end do 745 end do 746 end if 747 748 ngnt=0 749 750 !Loop on (lp1,mp1) 751 !------------------------------------------------ 752 do lp1=0,l_max-1 753 do mp1=-lp1,lp1 754 ilmp1=1+lp1**2+lp1+mp1 755 k0lm1=ilmp1*(ilmp1-1)/2 756 757 ! Loop on (l1,m1)<=(lp1,mp1) 758 ! ------------------------------------------------ 759 do l1=0,l_max-1 760 do m1=-l1,l1 761 ilm1=1+l1**2+l1+m1 762 763 if (ilm1<=ilmp1) then 764 765 klm1=k0lm1+ilm1 766 gntselect(:,klm1)=0 767 768 ! Loop on (l2,m2) 769 ! ------------------------------------------------ 770 do l2=abs(l1-lp1),l1+lp1,2 771 do m2=-l2,l2 772 ilm2=1+l2**2+l2+m2 773 774 ! Real Gaunt coeffs selection rules 775 ! ------------------------------------------------ 776 if ((l2<=l1+lp1).and.& 777 & (((m1== mp1).and.((m2==0).or.(m2==2*abs(mp1)))).or.& 778 & ((m1==-mp1).and.(m2==-abs(m1)-abs(mp1))).or.& 779 & ((abs(m1)/=(abs(mp1)).and.& 780 & ((m2==ssgn(ilm1)*ssgn(ilmp1)* (abs(m1)+abs(mp1))).or.& 781 & (m2==ssgn(ilm1)*ssgn(ilmp1)*abs(abs(m1)-abs(mp1)))& 782 ))))) then 783 784 ! Compute selected real Gaunt coefficient 785 ! ------------------------------------------------ 786 realgnt_tmp=zero 787 do mm1=-l1,l1 788 c11=coeff(l1+1)%value(1,l1+mm1+1,l1+m1+1) 789 c12=coeff(l1+1)%value(2,l1+mm1+1,l1+m1+1) 790 do mm2= -lp1,lp1 791 c21=coeff(lp1+1)%value(1,lp1+mm2+1,lp1+mp1+1) 792 c22=coeff(lp1+1)%value(2,lp1+mm2+1,lp1+mp1+1) 793 do mm3= -l2,l2 794 c31=coeff(l2+1)%value(1,l2+mm3+1,l2+m2+1) 795 c32=coeff(l2+1)%value(2,l2+mm3+1,l2+m2+1) 796 fact=c11*c21*c31 - c12*c22*c31& 797 & -c11*c22*c32 - c12*c21*c32 798 if((abs(fact)>=tol12).and.(mm3==-mm2-mm1)) & 799 & realgnt_tmp=realgnt_tmp+fact*(-1)**mm2 & 800 & *gaunt(l2,mm3,l1,mm1,lp1,-mm2) 801 end do 802 end do 803 end do 804 805 ! Count additional non-zero real Gaunt coeffs 806 ! ------------------------------------------------ 807 if (abs(realgnt_tmp)>=tol12) then 808 ngnt=ngnt+1 809 gntselect(ilm2,klm1)=ngnt 810 realgnt(ngnt)=realgnt_tmp/sqrt(four_pi) 811 end if 812 813 ! End loops 814 ! ------------------------------------------------ 815 end if 816 end do 817 end do 818 end if 819 end do 820 end do 821 end do 822 end do 823 824 !Deallocate memory 825 !------------------------------------------------ 826 do ll=1,4*l_max-3 827 LIBPAW_DEALLOCATE(coeff(ll)%value) 828 end do 829 LIBPAW_DATATYPE_DEALLOCATE(coeff) 830 LIBPAW_DEALLOCATE(ssgn) 831 832 end subroutine realgaunt
m_pawang/rfactorial [ Functions ]
[ Top ] [ m_pawang ] [ Functions ]
NAME
rfactorial
FUNCTION
Private function Calculates N! as a double precision real.
INPUTS
nn=number to use
OUTPUT
factorial= n! (real)
PARENTS
CHILDREN
SOURCE
1060 elemental function rfactorial(nn) 1061 1062 1063 !This section has been created automatically by the script Abilint (TD). 1064 !Do not modify the following lines by hand. 1065 #undef ABI_FUNC 1066 #define ABI_FUNC 'rfactorial' 1067 !End of the abilint section 1068 1069 implicit none 1070 1071 !Arguments --------------------------------------------- 1072 !scalars 1073 integer,intent(in) :: nn 1074 real(dp) :: rfactorial 1075 1076 !Local variables --------------------------------------- 1077 !scalars 1078 integer :: ii 1079 1080 ! ********************************************************************* 1081 1082 rfactorial=one 1083 do ii=2,nn 1084 rfactorial=rfactorial*ii 1085 end do 1086 1087 end function rfactorial