TABLE OF CONTENTS


ABINIT/m_pawang [ Modules ]

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