TABLE OF CONTENTS


ABINIT/GAMMA_FUNCTION [ Functions ]

[ Top ] [ Functions ]

NAME

  GAMMA_FUNCTION

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      pspnl_hgh_rec,pspnl_operat_rec,vso_realspace_local

CHILDREN

      gsl_f90_sf_gamma

SOURCE

1031 subroutine GAMMA_FUNCTION(X,GA)
1032 
1033 
1034 !This section has been created automatically by the script Abilint (TD).
1035 !Do not modify the following lines by hand.
1036 #undef ABI_FUNC
1037 #define ABI_FUNC 'GAMMA_FUNCTION'
1038 !End of the abilint section
1039 
1040   implicit none
1041 
1042 #ifdef HAVE_GSL
1043 ! in case we have gsl, no need to use explicit function, just wrap the
1044 !  call to the GSL C function in 01_gsl_ext/
1045 
1046   ! arguments
1047 
1048   real(dp),intent(in) :: x
1049   real(dp),intent(out) :: ga
1050 
1051   call gsl_f90_sf_gamma(x,ga)
1052 
1053 #else
1054 !       ====================================================
1055 !       Purpose: This program computes the gamma function
1056 !                Gamma(x) using subroutine GAMMA
1057 !       Examples:
1058 !                   x          Gamma(x)
1059 !                ----------------------------
1060 !                  1/3       2.678938534708
1061 !                  0.5       1.772453850906
1062 !                 -0.5      -3.544907701811
1063 !                 -1.5       2.363271801207
1064 !                  5.0      24.000000000000
1065 !       ====================================================
1066 !
1067 !  This routine was downloaded from UIUC:
1068 !  http://jin.ece.uiuc.edu/routines/routines.html
1069 !
1070 !  The programs appear to accompany a book "Computation of Special
1071 !  Functions" (1996) John Wiley and Sons, but are distributed online
1072 !  by the authors. Exact copyright should be checked.
1073 !
1074 !  Authors / copyright:
1075 !     Shanjie Zhang and Jianming Jin
1076 !     Proposed contact is:  j-jin1@uiuc.edu
1077 !
1078 !  20 October 2008:
1079 !     Incorporated into ABINIT by M. Verstraete
1080 !
1081 !
1082 !
1083 !       ==================================================
1084 !       Purpose: Compute the gamma function Gamma(x)
1085 !       Input :  x  --- Argument of Gamma(x)
1086 !                       ( x is not equal to 0,-1,-2, etc )
1087 !       Output:  GA --- Gamma(x)
1088 !       ==================================================
1089 !
1090 
1091   ! arguments
1092 
1093   real(dp),intent(in) :: x
1094   real(dp),intent(out) :: ga
1095 
1096   ! local variables
1097   integer :: k,m
1098   real(dp) :: m1,z,r,gr
1099   real(dp) :: G(26)
1100 
1101   ! source code:
1102 
1103   ! initialization of reference data
1104   G=(/1.0D0,0.5772156649015329D0, &
1105      &  -0.6558780715202538D0, -0.420026350340952D-1, &
1106      &   0.1665386113822915D0,-.421977345555443D-1, &
1107      &  -.96219715278770D-2, .72189432466630D-2, &
1108      &  -.11651675918591D-2, -.2152416741149D-3, &
1109      &   .1280502823882D-3, -.201348547807D-4, &
1110      &  -.12504934821D-5, .11330272320D-5, &
1111      &  -.2056338417D-6, .61160950D-8, &
1112      &   .50020075D-8, -.11812746D-8, &
1113      &   .1043427D-9, .77823D-11, &
1114      &  -.36968D-11, .51D-12, &
1115      &  -.206D-13, -.54D-14, .14D-14, .1D-15/)
1116 
1117 
1118   ! for the integer case, do explicit factorial
1119   if (X==int(X)) then
1120     if (X > 0.0D0) then
1121       GA=1.0D0
1122       M1=X-1
1123       do K=2,int(M1)
1124         GA=GA*K
1125       end do
1126     else
1127       GA=1.0D+300
1128     end if
1129   ! for the integer case, do explicit factorial
1130   else
1131     if (abs(X) > 1.0D0) then
1132       Z=abs(X)
1133       M=int(Z)
1134       R=1.0D0
1135       do K=1,M
1136         R=R*(Z-K)
1137       end do
1138       Z=Z-M
1139     else
1140       Z=X
1141     end if
1142     GR=G(26)
1143     do K=25,1,-1
1144       GR=GR*Z+G(K)
1145     end do
1146     GA=1.0D0/(GR*Z)
1147     if (abs(X) > 1.0D0) then
1148       GA=GA*R
1149       if (X < 0.0D0) GA=-PI/(X*GA*SIN(PI*X))
1150     end if
1151   end if
1152   return
1153 
1154 #endif
1155 !  end preproc for presence of GSL
1156 
1157 end subroutine GAMMA_FUNCTION

ABINIT/m_special_funcs [ Modules ]

[ Top ] [ Modules ]

NAME

 m_special_funcs

FUNCTION

 This module contains routines and functions used to
 evaluate special functions frequently needed in Abinit.

COPYRIGHT

 Copyright (C) 2008-2018 ABINIT group (MG, MT, FB, XG, MVer, FJ, NH, GZ, DRH)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 module m_special_funcs
25 
26  use defs_basis
27  use m_abicore
28  use m_errors
29  use m_splines
30 
31  use m_fstrings,        only : sjoin, ftoa
32  use m_numeric_tools,   only : arth, simpson
33 
34  implicit none
35 
36  private
37 
38  public :: clp               ! x-1, if x>1/2, x+1, if x<-1/2
39  public :: factorial         ! Calculates N! returning a real.
40  public :: permutations      ! Returns N!/(N-k) if N>=0 and N-k>0 else 0.
41  public :: binomcoeff        ! Binominal coefficient n!/(n-k)!
42  public :: laguerre          ! Laguerre Polynomial(x,n,a).
43  public :: RadFnH            ! Atomic radial function(r,n,l,Z).
44  public :: iradfnh           ! Norm of atomic radial function(a,b,n,l,Z).
45  public :: dirac_delta       ! Approximate Dirac delta with normal distribution.
46  public :: gaussian          ! Normalized Gaussian distribution.
47  public :: abi_derf          ! Evaluates the error function in real(dp).
48  public :: abi_derfc         ! Evaluates the complementary error function in real(dp).
49  public :: gamma_function    ! Computes the gamma function
50  public :: besjm             ! Spherical bessel function of order nn. Handles nn=0,1,2,3,4, or 5 only.
51  public :: sbf8              ! Computes set of spherical bessel functions using accurate algorithm
52  public :: k_fermi           ! Fermi wave vector corresponding to the local value of the real space density rhor.
53  public :: k_thfermi         ! Thomas-Fermi wave vector corresponding to the local value of the real space density rhor

m_special_funcs/abi_derf [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 abi_derf

FUNCTION

 Evaluates the error function in real(dp).
 Same implementation as imsl.
 Simple mod of derfc.F90

INPUTS

 yy

OUTPUT

 derf_yy= error function of yy

PARENTS

      evdw_wannier,wvl_wfs_set

CHILDREN

SOURCE

721 elemental function abi_derf(yy) result(derf_yy)
722 
723 
724 !This section has been created automatically by the script Abilint (TD).
725 !Do not modify the following lines by hand.
726 #undef ABI_FUNC
727 #define ABI_FUNC 'abi_derf'
728 !End of the abilint section
729 
730  implicit none
731 
732 !Arguments ------------------------------------
733 !scalars
734  real(dp),intent(in) :: yy
735  real(dp) :: derf_yy
736 
737 !Local variables-------------------------------
738  integer          ::  done,ii,isw
739 ! coefficients for 0.0 <= yy < .477
740  real(dp), parameter :: &
741 &  pp(5)=(/ 113.8641541510502e0_dp, 377.4852376853020e0_dp,  &
742 &           3209.377589138469e0_dp, .1857777061846032e0_dp,  &
743 &           3.161123743870566e0_dp /)
744  real(dp), parameter :: &
745 &  qq(4)=(/ 244.0246379344442e0_dp, 1282.616526077372e0_dp,  &
746 &           2844.236833439171e0_dp, 23.60129095234412e0_dp/)
747 ! coefficients for .477 <= yy <= 4.0
748  real(dp), parameter :: &
749 &  p1(9)=(/ 8.883149794388376e0_dp, 66.11919063714163e0_dp,  &
750 &           298.6351381974001e0_dp, 881.9522212417691e0_dp,  &
751 &           1712.047612634071e0_dp, 2051.078377826071e0_dp,  &
752 &           1230.339354797997e0_dp, 2.153115354744038e-8_dp, &
753 &           .5641884969886701e0_dp /)
754  real(dp), parameter :: &
755 &  q1(8)=(/ 117.6939508913125e0_dp, 537.1811018620099e0_dp,  &
756 &           1621.389574566690e0_dp, 3290.799235733460e0_dp,  &
757 &           4362.619090143247e0_dp, 3439.367674143722e0_dp,  &
758 &           1230.339354803749e0_dp, 15.74492611070983e0_dp/)
759   ! coefficients for 4.0 < y,
760  real(dp), parameter :: &
761 &  p2(6)=(/ -3.603448999498044e-01_dp, -1.257817261112292e-01_dp,   &
762 &           -1.608378514874228e-02_dp, -6.587491615298378e-04_dp,   &
763 &           -1.631538713730210e-02_dp, -3.053266349612323e-01_dp/)
764  real(dp), parameter :: &
765 &  q2(5)=(/ 1.872952849923460e0_dp   , 5.279051029514284e-01_dp,    &
766 &           6.051834131244132e-02_dp , 2.335204976268692e-03_dp,    &
767 &           2.568520192289822e0_dp /)
768  real(dp), parameter :: &
769 &  sqrpi=.5641895835477563e0_dp, xbig=13.3e0_dp, xlarge=6.375e0_dp, xmin=1.0e-10_dp
770  real(dp) ::  res,xden,xi,xnum,xsq,xx
771 
772 ! ******************************************************************
773 
774  xx = yy
775  isw = 1
776 !Here change the sign of xx, and keep track of it thanks to isw
777  if (xx<0.0e0_dp) then
778    isw = -1
779    xx = -xx
780  end if
781 
782  done=0
783 
784 !Residual value, if yy < -6.375e0_dp
785  res=-1.0e0_dp
786 
787 !abs(yy) < .477, evaluate approximation for erfc
788  if (xx<0.477e0_dp) then
789 !  xmin is a very small number
790    if (xx<xmin) then
791      res = xx*pp(3)/qq(3)
792    else
793      xsq = xx*xx
794      xnum = pp(4)*xsq+pp(5)
795      xden = xsq+qq(4)
796      do ii = 1,3
797        xnum = xnum*xsq+pp(ii)
798        xden = xden*xsq+qq(ii)
799      end do
800      res = xx*xnum/xden
801    end if
802    if (isw==-1) res = -res
803    done=1
804  end if
805 
806 !.477 < abs(yy) < 4.0 , evaluate approximation for erfc
807  if (xx<=4.0e0_dp .and. done==0 ) then
808    xsq = xx*xx
809    xnum = p1(8)*xx+p1(9)
810    xden = xx+q1(8)
811    do ii=1,7
812      xnum = xnum*xx+p1(ii)
813      xden = xden*xx+q1(ii)
814    end do
815    res = xnum/xden
816    res = res* exp(-xsq)
817    if (isw.eq.-1) then
818      res = res-1.0e0_dp
819    else
820      res=1.0e0_dp-res
821    end if
822    done=1
823  end if
824 
825 !y > 13.3e0_dp
826  if (isw > 0 .and. xx > xbig .and. done==0 ) then
827    res = 1.0e0_dp
828    done=1
829  end if
830 
831 !4.0 < yy < 13.3e0_dp  .or. -6.375e0_dp < yy < -4.0
832 !evaluate minimax approximation for erfc
833  if ( ( isw > 0 .or. xx < xlarge ) .and. done==0 ) then
834    xsq = xx*xx
835    xi = 1.0e0_dp/xsq
836    xnum= p2(5)*xi+p2(6)
837    xden = xi+q2(5)
838    do ii = 1,4
839      xnum = xnum*xi+p2(ii)
840      xden = xden*xi+q2(ii)
841    end do
842    res = (sqrpi+xi*xnum/xden)/xx
843    res = res* exp(-xsq)
844    if (isw.eq.-1) then
845      res = res-1.0e0_dp
846    else
847      res=1.0e0_dp-res
848    end if
849  end if
850 
851 !All cases have been investigated
852  derf_yy = res
853 
854 end function abi_derf

m_special_funcs/abi_derfc [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 abi_derfc

FUNCTION

 Evaluates the complementary error function in real(dp).
 Same implementation as imsl.

INPUTS

 yy

OUTPUT

 derfc_yy=complementary error function of yy

PARENTS

      ewald,ewald2,dfpt_ewald,elt_ewald,ewald9,make_efg_ion,psp2lo,wvl_wfs_set

CHILDREN

SOURCE

 880 elemental function abi_derfc(yy) result(derfc_yy)
 881 
 882 
 883 !This section has been created automatically by the script Abilint (TD).
 884 !Do not modify the following lines by hand.
 885 #undef ABI_FUNC
 886 #define ABI_FUNC 'abi_derfc'
 887 !End of the abilint section
 888 
 889  implicit none
 890 
 891 !Arguments ------------------------------------
 892 !scalars
 893  real(dp),intent(in) :: yy
 894  real(dp) :: derfc_yy
 895 
 896 !Local variables-------------------------------
 897  integer          ::  done,ii,isw
 898 ! coefficients for 0.0 <= yy < .477
 899  real(dp), parameter :: &
 900 &  pp(5)=(/ 113.8641541510502e0_dp, 377.4852376853020e0_dp,  &
 901 &           3209.377589138469e0_dp, .1857777061846032e0_dp,  &
 902 &           3.161123743870566e0_dp /)
 903  real(dp), parameter :: &
 904 &  qq(4)=(/ 244.0246379344442e0_dp, 1282.616526077372e0_dp,  &
 905 &           2844.236833439171e0_dp, 23.60129095234412e0_dp/)
 906 ! coefficients for .477 <= yy <= 4.0
 907  real(dp), parameter :: &
 908 &  p1(9)=(/ 8.883149794388376e0_dp, 66.11919063714163e0_dp,  &
 909 &           298.6351381974001e0_dp, 881.9522212417691e0_dp,  &
 910 &           1712.047612634071e0_dp, 2051.078377826071e0_dp,  &
 911 &           1230.339354797997e0_dp, 2.153115354744038e-8_dp, &
 912 &           .5641884969886701e0_dp /)
 913  real(dp), parameter :: &
 914 &  q1(8)=(/ 117.6939508913125e0_dp, 537.1811018620099e0_dp,  &
 915 &           1621.389574566690e0_dp, 3290.799235733460e0_dp,  &
 916 &           4362.619090143247e0_dp, 3439.367674143722e0_dp,  &
 917 &           1230.339354803749e0_dp, 15.74492611070983e0_dp/)
 918  ! coefficients for 4.0 < y,
 919  real(dp), parameter :: &
 920 &  p2(6)=(/ -3.603448999498044e-01_dp, -1.257817261112292e-01_dp,   &
 921 &           -1.608378514874228e-02_dp, -6.587491615298378e-04_dp,   &
 922 &           -1.631538713730210e-02_dp, -3.053266349612323e-01_dp/)
 923  real(dp), parameter :: &
 924 &  q2(5)=(/ 1.872952849923460e0_dp   , 5.279051029514284e-01_dp,    &
 925 &           6.051834131244132e-02_dp , 2.335204976268692e-03_dp,    &
 926 &           2.568520192289822e0_dp /)
 927  real(dp), parameter :: &
 928 &  sqrpi=.5641895835477563e0_dp, xbig=13.3e0_dp, xlarge=6.375e0_dp, xmin=1.0e-10_dp
 929  real(dp) ::  res,xden,xi,xnum,xsq,xx
 930 
 931 !******************************************************************
 932 
 933  xx = yy
 934  isw = 1
 935 !Here change the sign of xx, and keep track of it thanks to isw
 936  if (xx<0.0e0_dp) then
 937    isw = -1
 938    xx = -xx
 939  end if
 940 
 941  done=0
 942 
 943 !Residual value, if yy < -6.375e0_dp
 944  res=2.0e0_dp
 945 
 946 !abs(yy) < .477, evaluate approximation for erfc
 947  if (xx<0.477e0_dp) then
 948 !  xmin is a very small number
 949    if (xx<xmin) then
 950      res = xx*pp(3)/qq(3)
 951    else
 952      xsq = xx*xx
 953      xnum = pp(4)*xsq+pp(5)
 954      xden = xsq+qq(4)
 955      do ii = 1,3
 956        xnum = xnum*xsq+pp(ii)
 957        xden = xden*xsq+qq(ii)
 958      end do
 959      res = xx*xnum/xden
 960    end if
 961    if (isw==-1) res = -res
 962    res = 1.0e0_dp-res
 963    done=1
 964  end if
 965 
 966 !.477 < abs(yy) < 4.0 , evaluate approximation for erfc
 967  if (xx<=4.0e0_dp .and. done==0 ) then
 968    xsq = xx*xx
 969    xnum = p1(8)*xx+p1(9)
 970    xden = xx+q1(8)
 971    do ii=1,7
 972      xnum = xnum*xx+p1(ii)
 973      xden = xden*xx+q1(ii)
 974    end do
 975    res = xnum/xden
 976    res = res* exp(-xsq)
 977    if (isw.eq.-1) res = 2.0e0_dp-res
 978    done=1
 979  end if
 980 
 981 !y > 13.3e0_dp
 982  if (isw > 0 .and. xx > xbig .and. done==0 ) then
 983    res = 0.0e0_dp
 984    done=1
 985  end if
 986 
 987 !4.0 < yy < 13.3e0_dp  .or. -6.375e0_dp < yy < -4.0
 988 !evaluate minimax approximation for erfc
 989  if ( ( isw > 0 .or. xx < xlarge ) .and. done==0 ) then
 990    xsq = xx*xx
 991    xi = 1.0e0_dp/xsq
 992    xnum= p2(5)*xi+p2(6)
 993    xden = xi+q2(5)
 994    do ii = 1,4
 995      xnum = xnum*xi+p2(ii)
 996      xden = xden*xi+q2(ii)
 997    end do
 998    res = (sqrpi+xi*xnum/xden)/xx
 999    res = res* exp(-xsq)
1000    if (isw.eq.-1) res = 2.0e0_dp-res
1001  end if
1002 
1003 !All cases have been investigated
1004  derfc_yy = res
1005 
1006 end function abi_derfc

m_special_funcs/besjm [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 besjm

FUNCTION

 Spherical bessel function of order nn. Handles nn=0,1,2,3,4, or 5 only.

INPUTS

  arg= scaling to be applied to xx(nx)
  nn=order of spherical bessel function (only 0 through 5 allowed)
  cosx(1:nx)=cosines of arg*xx(1:nx)
  xx(1:nx)=set of dimensionless arguments of function
  nx=number of arguments
  sinx(1:nx)=sines of arg*xx(1:nx)

OUTPUT

  besjx(1:nx)=returned values

NOTES

 besj(nn,y)=$ j_{nn}(y) =(\frac{\pi}{2y})^{\frac{1}{2}}J(nn+\frac{1}{2},y)$
 where J=Bessel function of the first kind.
 besjm compute multiple values, and relies on precomputed values of sin and cos of y.
 The argument y is arg*xx(ix), for ix from 1 to nx
 The values of xx must be positive, and ordered by increasing order
 At small arg, the higher orders have so much cancellation that the
 analytic expression is very poor computationally.  In that case we
 use a rational polynomial approximation.

PARENTS

      m_special_funcs,mlwfovlp_radial,psp1nl

CHILDREN

      splint

SOURCE

1196 subroutine besjm(arg,besjx,cosx,nn,nx,sinx,xx)
1197 
1198 
1199 !This section has been created automatically by the script Abilint (TD).
1200 !Do not modify the following lines by hand.
1201 #undef ABI_FUNC
1202 #define ABI_FUNC 'besjm'
1203 !End of the abilint section
1204 
1205  implicit none
1206 
1207 !This section has been created automatically by the script Abilint (TD).
1208 !Do not modify the following lines by hand.
1209 #undef ABI_FUNC
1210 #define ABI_FUNC 'besjm'
1211 !End of the abilint section
1212 
1213 !Arguments ------------------------------------
1214 !scalars
1215  integer,intent(in) :: nn,nx
1216  real(dp),intent(in) :: arg
1217 !arrays
1218  real(dp),intent(in) :: cosx(nx),sinx(nx),xx(nx)
1219  real(dp),intent(out) :: besjx(nx)
1220 
1221 !Local variables-------------------------------
1222 !scalars
1223  integer :: ix,switchx
1224 !Series or rational polynomial coefficients
1225  real(dp),parameter :: b01=1.d0/6.d0,b02=1.d0/120.d0,b03=1.d0/5040.d0
1226  real(dp),parameter :: b04=1.d0/362880.d0,b11=0.8331251468724171d-1
1227  real(dp),parameter :: b12=0.2036961284395412d-2,b13=0.1932970379901801d-4
1228  real(dp),parameter :: b14=0.6526053169009489d-7,b21=0.5867824627555163d-1
1229  real(dp),parameter :: b22=0.1152501878595934d-2,b23=0.1011071389414764d-4
1230  real(dp),parameter :: b24=0.4172322111421287d-7,b25=0.6790616688656543d-10
1231  real(dp),parameter :: b31=0.439131885807176d-1,b32=0.6813139609887099d-3
1232  real(dp),parameter :: b33=0.4899103784264755d-5,b34=0.17025590795625d-7
1233  real(dp),parameter :: b35=0.2382642910613347d-10,b41=0.3587477991030971d-1
1234  real(dp),parameter :: b42=0.4833719855268907d-3,b43=0.3238388977796242d-5
1235  real(dp),parameter :: b44=0.1171802513125112d-7,b45=0.223261650431992d-10
1236  real(dp),parameter :: b46=.1800045587335951d-13,b51=0.295232406376567d-1
1237  real(dp),parameter :: b52=0.3359864457080573d-3,b53=0.19394750603618d-5
1238  real(dp),parameter :: b54=0.6143166228216219d-8,b55=0.10378501636108d-10
1239  real(dp),parameter :: b56=.749975122872713d-14
1240  real(dp),parameter :: c11=0.1668748531275829d-1,c12=0.1342812442426702d-3
1241  real(dp),parameter :: c13=0.6378249315355233d-6,c14=0.1573564527360138d-8
1242  real(dp),parameter :: c21=0.127503251530198d-1,c22=0.7911240539893565d-4
1243  real(dp),parameter :: c23=0.3044380758068054d-6,c24=0.7439837832363479d-9
1244  real(dp),parameter :: c25=0.9515065658793124d-12,c31=0.1164236697483795d-1
1245  real(dp),parameter :: c32=0.654858636312224d-4,c33=0.2265576367562734d-6
1246  real(dp),parameter :: c34=0.4929905563217352d-9,c35=0.555120465710914d-12
1247  real(dp),parameter :: c41=0.9579765544235745d-2,c42=0.4468999977536864d-4
1248  real(dp),parameter :: c43=0.1315634305905896d-6,c44=0.2615492488301639d-9
1249  real(dp),parameter :: c45=0.3387473312408129d-12,c46=.2280866204624012d-15
1250  real(dp),parameter :: c51=0.8938297823881763d-2,c52=0.3874149021633025d-4
1251  real(dp),parameter :: c53=0.1054692715135225d-6,c54=0.192879620987602d-9
1252  real(dp),parameter :: c55=0.2284469423833734d-12,c56=0.139729234332572d-15
1253  real(dp),parameter :: ffnth=1.d0/15.d0,o10395=1.d0/10395d0,oo105=1.d0/105.d0
1254  real(dp),parameter :: oo945=1.d0/945.d0
1255  real(dp) :: bot,rr,rsq,top
1256  character(len=500) :: message
1257 
1258 ! *************************************************************************
1259 
1260  if (nn==0) then
1261 
1262    switchx=nx+1
1263    do ix=1,nx
1264      rr=arg*xx(ix)
1265      if (rr<=1.d-1) then
1266        rsq=rr*rr
1267        besjx(ix)=1.d0-rsq*(b01-rsq*(b02-rsq*(b03-rsq*b04)))
1268      else
1269        switchx=ix
1270        exit
1271      end if
1272    end do
1273 
1274    do ix=switchx,nx
1275      rr=arg*xx(ix)
1276      besjx(ix)=sinx(ix)/rr
1277    end do
1278 
1279  else if (nn==1) then
1280 
1281    switchx=nx+1
1282    do ix=1,nx
1283      rr=arg*xx(ix)
1284      if (rr<=1.d0) then
1285        rsq=rr*rr
1286        top=1.d0-rsq*(b11-rsq*(b12-rsq*(b13-rsq*b14)))
1287        bot=1.d0+rsq*(c11+rsq*(c12+rsq*(c13+rsq*c14)))
1288        besjx(ix)=third*rr*top/bot
1289      else
1290        switchx=ix
1291        exit
1292      end if
1293    end do
1294 
1295    do ix=switchx,nx
1296      rr=arg*xx(ix)
1297      rsq=rr*rr
1298      besjx(ix)=(sinx(ix)-rr*cosx(ix))/rsq
1299    end do
1300 
1301  else if (nn==2) then
1302 
1303    switchx=nx+1
1304    do ix=1,nx
1305      rr=arg*xx(ix)
1306      if (rr<=2.d0) then
1307        rsq=rr*rr
1308        top=1.d0-rsq*(b21-rsq*(b22-rsq*(b23-rsq*(b24-rsq*b25))))
1309        bot=1.d0+rsq*(c21+rsq*(c22+rsq*(c23+rsq*(c24+rsq*c25))))
1310        besjx(ix)=ffnth*rsq*top/bot
1311      else
1312        switchx=ix
1313        exit
1314      end if
1315    end do
1316 
1317    do ix=switchx,nx
1318      rr=arg*xx(ix)
1319      rsq=rr*rr
1320      besjx(ix)=((3.d0-rsq)*sinx(ix)-3.d0*rr*cosx(ix))/(rr*rsq)
1321    end do
1322 
1323  else if (nn==3) then
1324 
1325    switchx=nx+1
1326    do ix=1,nx
1327      rr=arg*xx(ix)
1328      if (rr<=2.d0) then
1329        rsq=rr*rr
1330        top=1.d0-rsq*(b31-rsq*(b32-rsq*(b33-rsq*(b34-rsq*b35))))
1331        bot=1.d0+rsq*(c31+rsq*(c32+rsq*(c33+rsq*(c34+rsq*c35))))
1332        besjx(ix)=rr*rsq*oo105*top/bot
1333      else
1334        switchx=ix
1335        exit
1336      end if
1337    end do
1338 
1339    do ix=switchx,nx
1340      rr=arg*xx(ix)
1341      rsq=rr*rr
1342      besjx(ix)=( (15.d0-6.d0*rsq)*sinx(ix)&
1343 &     + rr*(rsq-15.d0)  *cosx(ix) ) /(rsq*rsq)
1344    end do
1345 
1346  else if (nn==4) then
1347 
1348    switchx=nx+1
1349    do ix=1,nx
1350      rr=arg*xx(ix)
1351      if (rr<=4.d0) then
1352        rsq=rr*rr
1353        top=1.d0-rsq*(b41-rsq*(b42-rsq*(b43-rsq*(b44-rsq*(b45-rsq*b46)))))
1354        bot=1.d0+rsq*(c41+rsq*(c42+rsq*(c43+rsq*(c44+rsq*(c45+rsq*c46)))))
1355        besjx(ix)=rsq*rsq*oo945*top/bot
1356      else
1357        switchx=ix
1358        exit
1359      end if
1360    end do
1361 
1362    do ix=switchx,nx
1363      rr=arg*xx(ix)
1364      rsq=rr*rr
1365      besjx(ix)=( (105.d0-rsq*(45.d0-rsq)) *sinx(ix)&
1366 &     + rr * (10.d0*rsq-105.d0)  *cosx(ix) ) /(rsq*rsq*rr)
1367    end do
1368 
1369  else if (nn==5) then
1370 
1371    switchx=nx+1
1372    do ix=1,nx
1373      rr=arg*xx(ix)
1374      if (rr<=4.d0) then
1375        rsq=rr*rr
1376        top=1.d0-rsq*(b51-rsq*(b52-rsq*(b53-rsq*(b54-rsq*(b55-rsq*b56)))))
1377        bot=1.d0+rsq*(c51+rsq*(c52+rsq*(c53+rsq*(c54+rsq*(c55+rsq*c56)))))
1378        besjx(ix)=rsq*rsq*rr*o10395*top/bot
1379      else
1380        switchx=ix
1381        exit
1382      end if
1383    end do
1384 
1385    do ix=switchx,nx
1386      rr=arg*xx(ix)
1387      rsq=rr*rr
1388      besjx(ix)=( (945.d0-rsq*(420.d0-rsq*15.d0)) *sinx(ix)&
1389 &     + rr * (945.d0-rsq*(105.d0-rsq))  *cosx(ix) ) /(rsq*rsq*rr)
1390    end do
1391 
1392  else
1393    write(message, '(a,i0,a)' )' besjm only defined for nn in [0,5]; input was nn=',nn,'.'
1394    MSG_BUG(message)
1395  end if
1396 
1397 end subroutine besjm

m_special_funcs/binomcoeff [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 factorial

FUNCTION

 Calculates n!/( k!* (n-k)!). Returns a real (dp)

INPUTS

   nn=number to use

OUTPUT

   binomcoeff= n!/( k!* (n-k)!)  (real dp)

PARENTS

CHILDREN

SOURCE

329 elemental function binomcoeff(n,k)
330 
331 
332 !This section has been created automatically by the script Abilint (TD).
333 !Do not modify the following lines by hand.
334 #undef ABI_FUNC
335 #define ABI_FUNC 'binomcoeff'
336 !End of the abilint section
337 
338  implicit none
339 
340 !Arguments ---------------------------------------------
341 !scalars
342  integer,intent(in) :: n,k
343  real(dp) :: binomcoeff
344 
345 ! *********************************************************************
346 
347  binomcoeff=factorial(n)/(factorial(k)*factorial(n-k))
348 
349 end function binomcoeff

m_special_funcs/bose_einstein [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

  bose_einstein

FUNCTION

  Returns the Bose Einstein distribution for T and energy
  presumes everything is in Hartree!!!! Not Kelvin for T

INPUTS

  energy = electron energy level
  temperature = T

PARENTS

SOURCE

1563 function bose_einstein(energy, temperature)
1564 
1565 
1566 !This section has been created automatically by the script Abilint (TD).
1567 !Do not modify the following lines by hand.
1568 #undef ABI_FUNC
1569 #define ABI_FUNC 'bose_einstein'
1570 !End of the abilint section
1571 
1572  implicit none
1573 
1574 !Arguments ------------------------------------
1575 !scalars
1576  real(dp),intent(in) :: energy, temperature
1577  real(dp) :: bose_einstein
1578 
1579 !Local variables-------------------------------
1580 !scalars
1581  real(dp) :: arg
1582  character(len=500) :: message
1583 
1584 ! *************************************************************************
1585 
1586  bose_einstein = zero
1587  if (temperature > tol12) then
1588    arg = energy/temperature
1589    if(arg > tol12 .and. arg < 600._dp)then
1590      bose_einstein = one / (exp(arg)  - one)
1591    else if (arg < tol12) then
1592      write(message,'(a)') 'No Bose Einstein for negative energies'
1593      MSG_WARNING(message)
1594    end if
1595  else if (arg < tol12) then
1596    write(message,'(a)') 'No Bose Einstein for negative or 0 T'
1597    MSG_WARNING(message)
1598  end if
1599 
1600 
1601 end function bose_einstein

m_special_funcs/clp [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 clp

FUNCTION

 clp(x)= x-1, if x>1/2
         x+1, if x<-1/2

INPUTS

  x= input variable

OUTPUT

  clp= resulting function

PARENTS

      nhatgrid

SOURCE

166 pure function clp(x)
167 
168 
169 !This section has been created automatically by the script Abilint (TD).
170 !Do not modify the following lines by hand.
171 #undef ABI_FUNC
172 #define ABI_FUNC 'clp'
173 !End of the abilint section
174 
175  implicit none
176 
177 !Arguments ------------------------------------
178 !scalars
179  real(dp) :: clp
180  real(dp),intent(in) :: x
181 
182 ! **********************************************************************
183 
184  if(x > half) then
185    clp=x-one
186  elseif(x < -half) then
187    clp=x+one
188  else
189    clp=x
190  end if
191 
192 end function clp

m_special_funcs/dirac_delta [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 dirac_delta

FUNCTION

  Approximate Dirac delta with normal distribution.
    delta(x,sigma) = 1/(sigma sqrt(pi)) e^{-x**2/(sigma**2)}

INPUTS

   arg=Argument of the approximated Delta.
   sigma=Broadening factor.

PARENTS

CHILDREN

SOURCE

621 elemental function dirac_delta(arg,sigma)
622 
623 
624 !This section has been created automatically by the script Abilint (TD).
625 !Do not modify the following lines by hand.
626 #undef ABI_FUNC
627 #define ABI_FUNC 'dirac_delta'
628 !End of the abilint section
629 
630  implicit none
631 
632 !Arguments ---------------------------------------------
633 !scalars
634  real(dp),intent(in) :: arg,sigma
635  real(dp) :: dirac_delta
636 
637 !Local variables ---------------------------------------
638 !scalars
639  real(dp) :: xx
640 
641 ! *********************************************************************
642 
643  xx=arg/sigma
644  dirac_delta = exp(-xx*xx) / (sigma*sqrt(pi))
645 
646 end function dirac_delta

m_special_funcs/factorial [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 factorial

FUNCTION

 Calculates N!. Returns a (dp) real.

INPUTS

   nn=number to use

OUTPUT

   factorial= n! (real)

PARENTS

CHILDREN

SOURCE

214 elemental function factorial(nn)
215 
216 
217 !This section has been created automatically by the script Abilint (TD).
218 !Do not modify the following lines by hand.
219 #undef ABI_FUNC
220 #define ABI_FUNC 'factorial'
221 !End of the abilint section
222 
223  implicit none
224 
225 !Arguments ---------------------------------------------
226 !scalars
227  integer,intent(in) :: nn
228  real(dp) :: factorial
229 
230 !Local variables ---------------------------------------
231 !scalars
232  integer :: ii
233  real(dp) :: ff
234 
235 ! *********************************************************************
236 
237  ff=one
238  do ii=2,nn
239    ff=ff*ii
240  end do
241 
242  factorial=ff
243 
244 end function factorial

m_special_funcs/fermi_dirac [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

  fermi_dirac

FUNCTION

  Returns the Fermi Dirac distribution for T and energy wrt Fermi level
  presumes everything is in Hartree!!!! Not Kelvin for T

INPUTS

  energy = electron energy level
  mu = chemical potential
  temperature = T

PARENTS

SOURCE

1508 function fermi_dirac(energy, mu, temperature)
1509 
1510 
1511 !This section has been created automatically by the script Abilint (TD).
1512 !Do not modify the following lines by hand.
1513 #undef ABI_FUNC
1514 #define ABI_FUNC 'fermi_dirac'
1515 !End of the abilint section
1516 
1517  implicit none
1518 
1519 !Arguments ------------------------------------
1520 !scalars
1521  real(dp),intent(in) :: energy, mu, temperature
1522  real(dp) :: fermi_dirac
1523 
1524 !Local variables-------------------------------
1525 !scalars
1526  real(dp) :: arg
1527 
1528 ! *************************************************************************
1529 
1530  fermi_dirac = zero
1531  if (temperature > tol12) then
1532    arg = (energy-mu)/temperature
1533    if(arg < -600._dp)then ! far below Ef
1534      fermi_dirac = one
1535    else if (arg < 600._dp)then ! around Ef
1536      fermi_dirac = one / (exp(arg)  + one)
1537    end if
1538  else  ! T is too small - just step function
1539    if (mu-energy > tol12) fermi_dirac = one
1540  end if
1541 
1542 end function fermi_dirac

m_special_funcs/gaussian [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 gaussian

FUNCTION

  Return the values of the normalized Gaussian distribution
    Gauss(arg,sigma) =  1/(sigma SQRT(2*pi)) e^{-arg**2/(2*sigma**2)}

INPUTS

   arg=Argument of the Gaussian.
   sigma=Standard deviation

PARENTS

CHILDREN

SOURCE

669 elemental function gaussian(arg,sigma)
670 
671 
672 !This section has been created automatically by the script Abilint (TD).
673 !Do not modify the following lines by hand.
674 #undef ABI_FUNC
675 #define ABI_FUNC 'gaussian'
676 !End of the abilint section
677 
678  implicit none
679 
680 !Arguments ---------------------------------------------
681 !scalars
682  real(dp),intent(in) :: arg,sigma
683  real(dp) :: gaussian
684 
685 !Local variables ---------------------------------------
686 !scalars
687  real(dp) :: xx
688 
689 ! *********************************************************************
690 
691  xx=arg/(sqrt2*sigma)
692  gaussian = exp(-xx*xx) / (sigma*sqrt(two_pi))
693 
694 end function gaussian

m_special_funcs/gspline_eval [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

  gspline_eval

FUNCTION

  Evaluate the gaussian approximant and its primitive at (xmesh - x0)

INPUTS

  self<gspline_t>=Object used to spline the gaussian approximant
  x0=Shift to be given to xmesh
  nx=Number of points in input mesh.
  xmesh(nx)=Frequency points (not necessarly linear).

OUTPUT

  weights(nx,2)=First slice contains the gaussian approximant on xmesh.
   The second slice stores the primitive.

PARENTS

CHILDREN

SOURCE

1977 pure subroutine gspline_eval(self, x0, nx, xmesh, weights)
1978 
1979 
1980 !This section has been created automatically by the script Abilint (TD).
1981 !Do not modify the following lines by hand.
1982 #undef ABI_FUNC
1983 #define ABI_FUNC 'gspline_eval'
1984 !End of the abilint section
1985 
1986  implicit none
1987 
1988 !Arguments ------------------------------------
1989 !scalars
1990  integer,intent(in) :: nx
1991  real(dp),intent(in) :: x0
1992  type(gspline_t),intent(in) :: self
1993 !arrays
1994  real(dp),intent(in) :: xmesh(nx)
1995  real(dp),intent(out) :: weights(nx,2)
1996 
1997 !Local variables ------------------------------
1998 !scalars
1999  integer :: ix,jspl
2000  real(dp) :: xx,absx,aa,bb,cc,dd
2001  logical :: isneg
2002  !real(dp) :: int_values(nx)
2003 
2004 ! *************************************************************************
2005 
2006  do ix=1,nx
2007    xx = xmesh(ix) - x0; absx = abs(xx); isneg = xx < zero
2008    if (absx >= self%xmax) then
2009      ! Region in which gauss(x) is negligible.
2010      weights(ix,1) = zero
2011      if (isneg) then
2012        weights(ix,2) = zero
2013      else
2014        weights(ix,2) = one
2015      end if
2016    else
2017      ! Spline functions at |x| and recover the value at x:
2018      ! g(x) = g(-x); G(-x) = 1 - G(x)
2019      jspl = 1 + int((absx - self%xmin) * self%stepm1); dd = absx - self%xvals(jspl)
2020      bb = dd * self%stepm1
2021      aa = one - bb
2022      cc = aa*(aa**2-one) * self%step2div6
2023      dd = bb*(bb**2-one) * self%step2div6
2024 
2025      weights(ix,1) = aa*self%svals(jspl,1) + bb*self%svals(jspl+1,1) + cc*self%svals(jspl,2) + dd*self%svals(jspl+1,2)
2026      weights(ix,2) = aa*self%svals(jspl,3) + bb*self%svals(jspl+1,3) + cc*self%svals(jspl,4) + dd*self%svals(jspl+1,4)
2027      if (isneg) weights(ix,2) = one - weights(ix,2)
2028    end if
2029  end do
2030 
2031  !call simpson_int(nx,xmesh(2) - xmesh(1),weights(:,1),int_values)
2032  !do ix=1,nx
2033  !  write(99,*)xmesh(ix), weights(ix,1), dirac_delta(xx, self%sigma), weights(ix,2), int_values(ix)
2034  !end do
2035 
2036 end subroutine gspline_eval

m_special_funcs/gspline_free [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

  gspline_free

FUNCTION

  Free dynamic memory

INPUTS

  self<gspline_t>=Object used to spline the gaussian approximant

PARENTS

CHILDREN

SOURCE

2055 subroutine gspline_free(self)
2056 
2057 
2058 !This section has been created automatically by the script Abilint (TD).
2059 !Do not modify the following lines by hand.
2060 #undef ABI_FUNC
2061 #define ABI_FUNC 'gspline_free'
2062 !End of the abilint section
2063 
2064  implicit none
2065 
2066 !Arguments ------------------------------------
2067 !scalars
2068  type(gspline_t),intent(inout) :: self
2069 
2070 ! *************************************************************************
2071 
2072  if (allocated(self%xvals)) then
2073    ABI_FREE(self%xvals)
2074  end if
2075  if (allocated(self%svals)) then
2076    ABI_FREE(self%svals)
2077  end if
2078 
2079 end subroutine gspline_free

m_special_funcs/gspline_new [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

  gspline_new

FUNCTION

  Build object to spline the gaussian approximant and its primitive.

INPUTS

  sigma=Broadening parameter.

PARENTS

CHILDREN

SOURCE

1905 type (gspline_t) function gspline_new(sigma) result(new)
1906 
1907 
1908 !This section has been created automatically by the script Abilint (TD).
1909 !Do not modify the following lines by hand.
1910 #undef ABI_FUNC
1911 #define ABI_FUNC 'gspline_new'
1912 !End of the abilint section
1913 
1914  implicit none
1915 
1916 !Arguments ------------------------------------
1917 !scalars
1918  real(dp),intent(in) :: sigma
1919 
1920 !Local variables ------------------------------
1921 !scalars
1922  integer :: ii
1923  real(dp) :: ybcbeg, ybcend
1924 ! *************************************************************************
1925 
1926  new%nspline = 5 * 1024; new%sigma = sigma
1927  ABI_CHECK(sigma > zero, sjoin("invalid sigma:", ftoa(sigma)))
1928  new%xmin = zero
1929  new%xmax = sigma * sqrt(-log(sigma * sqrt(pi) * tol12)) ! gauss(xmax) = tol12
1930  new%step = (new%xmax - new%xmin) / (new%nspline - 1)
1931  new%stepm1 = one / new%step; new%step2div6 = new%step**2 / six
1932 
1933  ABI_MALLOC(new%xvals, (new%nspline))
1934  do ii=1,new%nspline
1935    new%xvals(ii) = new%xmin + (ii-1) * new%step
1936  end do
1937  new%xmax = new%xvals(new%nspline)
1938 
1939  ! Spline the gaussian approximant.
1940  ABI_MALLOC(new%svals, (new%nspline, 4))
1941  new%svals(:, 1) = dirac_delta(new%xvals, sigma)
1942  ybcbeg = - (two * new%xmin / sigma**2) * new%svals(1,1)
1943  ybcend = - (two * new%xmax / sigma**2) * new%svals(new%nspline,1)
1944  call spline(new%xvals, new%svals(:,1), new%nspline, ybcbeg, ybcend, new%svals(:,2))
1945 
1946  ! Spline the primitive: 1/2 [1 + erf(x/sigma)]
1947  new%svals(:, 3) = half * (one + abi_derf(new%xvals / new%sigma))
1948  call spline(new%xvals, new%svals(:,3), new%nspline, new%svals(1,1), new%svals(new%nspline, 1), new%svals(:,4))
1949  !do ii=1,new%nspline; write(98,*)new%xvals(ii),new%svals(ii,3),new%svals(ii,4); end do
1950 
1951 end function gspline_new

m_special_funcs/gspline_t [ Types ]

[ Top ] [ m_special_funcs ] [ Types ]

NAME

 gspline_t

FUNCTION

  Object used to interpolate the gaussian approximant and its primitive with cubic spline.
  Particularly useful if we are computing DOSes with many k-points/bands
  because one can significantly decrease the number of calls to exponential functions.

SOURCE

112  type,public :: gspline_t
113 
114    integer :: nspline
115     ! Number of points used in spline table.
116 
117    real(dp) :: sigma
118     ! Broadening parameter.
119 
120    real(dp) :: xmin,xmax
121     ! Min and max x in spline mesh. Only positive xs are stored in memory
122     ! The values at -x are reconstructed by symmetry.
123     ! xmin is usually zero, xmax is the point where the gaussian == tol16.
124     ! g(x) is set to zero if x > xmin.
125 
126    real(dp) :: step, stepm1, step2div6
127     ! Step of the linear mesh used in spline and associated coeffients.
128 
129    real(dp),allocatable :: xvals(:)
130     ! xvals(nspline)
131     ! The xvalues used in the spline
132 
133    real(dp),allocatable :: svals(:,:)
134     ! svals(nspline,4)
135     ! Internal tables with spline data.
136 
137  end type gspline_t

m_special_funcs/IRadFnH [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 IRadFnH

FUNCTION

  IRadFnH(a,b,n,l,Z): Integral of radial function of atomic wavefunction between a and b.
  recursive programming using simpson's rule
  iteration depth of m=8 corresponds to relative error of 10^(-12).

INPUTS

   a lower limit for integration
   b upper limit for integration
   n principal quantum number
   l quantum number

OUTPUT

  IRadFnH(a,b,n,l,Z) (dp)

PARENTS

CHILDREN

  Laguerre
  factorial
  RadFnH

SOURCE

535 recursive function IRadFnH(a,b,n,l,Z,m) result(x)
536 
537 
538 !This section has been created automatically by the script Abilint (TD).
539 !Do not modify the following lines by hand.
540 #undef ABI_FUNC
541 #define ABI_FUNC 'IRadFnH'
542 !End of the abilint section
543 
544  implicit none
545 
546 !Arguments ---------------------------------------------
547 !scalars
548  integer,intent(in),optional  :: n,l,m
549  real(dp),intent(in):: a
550  real(dp),intent(in),optional :: b,Z
551 
552 !Local variables ---------------------------------------
553 !scalars
554  integer   :: nn,ll,mm
555  real(dp)  :: h,bb,ZZ,x
556 
557 ! *********************************************************************
558 
559  if (present(n)) then
560    nn=n
561  else
562    nn=3
563  end if
564 
565  if (present(l)) then
566    ll=l
567  else
568    ll=2
569  end if
570 
571  if (present(Z)) then
572    ZZ=Z
573  else
574    ZZ=28
575  end if
576 
577  if (present(b)) then
578    bb=b
579  else
580    bb=100.0_dp
581  end if
582 
583  if (present(m)) then
584    mm=m
585  else
586    mm=0
587  end if
588 
589  h=(bb-a)/2.0_dp
590  if (mm<8) then
591   !h=2*h/exp(1.0_dp)
592   x=IRadFnH(a,a+h,nn,ll,ZZ,mm+1)+IRadFnH(a+h,bb,nn,ll,ZZ,mm+1)
593  else
594   x=RadFnH(a,nn,ll,ZZ)**2*a**2+4.0_dp*RadFnH(a+h,nn,ll,ZZ)**2*(a+h)**2
595   x=h/3.0_dp*(x+RadFnH(bb,nn,ll,ZZ)**2*bb**2)
596  end if
597 
598 end function IRadFnH

m_special_funcs/jlspline_free [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 jlspline_free

FUNCTION

  deallocate memory

PARENTS

      m_cut3d,partial_dos_fractions

CHILDREN

      splint

SOURCE

1802 subroutine jlspline_free(jlspl)
1803 
1804 
1805 !This section has been created automatically by the script Abilint (TD).
1806 !Do not modify the following lines by hand.
1807 #undef ABI_FUNC
1808 #define ABI_FUNC 'jlspline_free'
1809 !End of the abilint section
1810 
1811  implicit none
1812 
1813 !Arguments ------------------------------------
1814  type(jlspline_t),intent(inout) :: jlspl
1815 
1816 ! *********************************************************************
1817 
1818  if (allocated(jlspl%xx)) then
1819    ABI_FREE(jlspl%xx)
1820  end if
1821  if (allocated(jlspl%bess_spl)) then
1822    ABI_FREE(jlspl%bess_spl)
1823  end if
1824  if (allocated(jlspl%bess_spl_der)) then
1825    ABI_FREE(jlspl%bess_spl_der)
1826  end if
1827 
1828 end subroutine jlspline_free

m_special_funcs/jlspline_integral [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 jlspline_integral

INPUTS

OUTPUT

FUNCTION

PARENTS

SOURCE

1847 real(dp) function jlspline_integral(jlspl, il, qq, powr, nr, rcut)  result(res)
1848 
1849 
1850 !This section has been created automatically by the script Abilint (TD).
1851 !Do not modify the following lines by hand.
1852 #undef ABI_FUNC
1853 #define ABI_FUNC 'jlspline_integral'
1854 !End of the abilint section
1855 
1856  implicit none
1857 
1858 !Arguments ------------------------------------
1859  integer,intent(in) :: il,nr,powr
1860  real(dp),intent(in) :: qq, rcut
1861  type(jlspline_t),intent(in) :: jlspl
1862 
1863 !Local variables ---------------------------------------
1864  integer :: ierr
1865  real(dp) :: step
1866 !arrays
1867  real(dp):: xfit(nr),yfit(nr),rr(nr)
1868 ! *********************************************************************
1869 
1870  step = rcut / (nr - 1)
1871  rr = arth(zero, step, nr)
1872  xfit = qq * rr
1873  call splint(jlspl%nx, jlspl%xx, jlspl%bess_spl(:,il), jlspl%bess_spl_der(:,il), nr, xfit, yfit, ierr=ierr)
1874 
1875  if (ierr /= 0) then
1876    write(std_out,*)"qq, rcut, qq*rcut, maxarg", qq, rcut, qq*rcut, jlspl%maxarg
1877    write(std_out,*)"x[0], x[-1]",jlspl%xx(1),jlspl%xx(jlspl%nx)
1878    write(std_out,*)"minval xfit: ",minval(xfit)
1879    write(std_out,*)"maxval xfit: ",maxval(xfit)
1880    MSG_ERROR("splint returned ierr != 0")
1881  end if
1882 
1883  if (powr /= 1) yfit = yfit * (rr ** powr)
1884  res = simpson(step, yfit)
1885 
1886 end function jlspline_integral

m_special_funcs/jlspline_new [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 jlspline_new

FUNCTION

 Pre-calculate the j_v(y) for recip_ylm on regular grid
     NOTE: spherical Bessel function small j!

INPUTS

  nx = max number of points on grid for integral
  delta = space between integral arguments
  mlang= max angular momentum

OUTPUT

  bess_spl=array of integrals
  bess_spl_der=array of derivatives of integrals
  xx=coordinates of points belonging to the grid

PARENTS

      m_cut3d,partial_dos_fractions

CHILDREN

      besjm,spline

SOURCE

1716 type(jlspline_t) function jlspline_new(nx, delta, mlang) result(new)
1717 
1718 
1719 !This section has been created automatically by the script Abilint (TD).
1720 !Do not modify the following lines by hand.
1721 #undef ABI_FUNC
1722 #define ABI_FUNC 'jlspline_new'
1723 !End of the abilint section
1724 
1725  implicit none
1726 
1727 !Arguments ------------------------------------
1728 !scalars
1729  integer,intent(in) :: nx,mlang
1730  real(dp),intent(in) :: delta
1731 
1732 !Local variables -------------------------
1733 !scalars
1734  integer :: ix,ll
1735  real(dp) :: yp1,ypn
1736 !arrays
1737  real(dp),allocatable :: cosbessx(:),sinbessx(:)
1738 
1739 ! *********************************************************************
1740 
1741  if (nx < 2) then
1742    MSG_ERROR('need more than one point for the interpolation routines')
1743  end if
1744 
1745  new%nx = nx; new%mlang = mlang; new%delta = delta; new%maxarg = (nx-1) * delta
1746  ABI_MALLOC(new%xx, (nx))
1747  ABI_MALLOC(new%bess_spl, (nx, mlang))
1748  ABI_MALLOC(new%bess_spl_der, (nx, mlang))
1749 
1750  !-----------------------------------------------------------------
1751  !Bessel function into array
1752  !-----------------------------------------------------------------
1753  ! integration grid is nfiner times finer than the interpolation grid
1754  ABI_MALLOC(sinbessx, (nx))
1755  ABI_MALLOC(cosbessx, (nx))
1756 
1757  ! could be done by chain rule for cos sin (is it worth it?) but
1758  ! precision problems as numerical errors are propagated.
1759  do ix=1,nx
1760    new%xx(ix) = (ix-1) * delta
1761    sinbessx(ix) = sin(new%xx(ix))
1762    cosbessx(ix) = cos(new%xx(ix))
1763  end do
1764 
1765  ! fill bess_spl array
1766  do ll=0,mlang-1
1767    call besjm(one,new%bess_spl(:,ll+1),cosbessx,ll,nx,sinbessx,new%xx)
1768 
1769    ! call spline to get 2nd derivative (reuse in splint later)
1770    yp1 = zero; ypn = zero
1771    call spline(new%xx, new%bess_spl(:,ll+1), nx, yp1, ypn, new%bess_spl_der(:,ll+1))
1772  end do
1773 
1774 !write(std_out,*) ' bess funct  0   1   2   3   4'
1775 !do ix=1,nx
1776 !write(std_out,*) xx(ix), (new%bess_spl(ix,ll),ll=1,mlang)
1777 !end do
1778 
1779  ABI_FREE(sinbessx)
1780  ABI_FREE(cosbessx)
1781 
1782 end function jlspline_new

m_special_funcs/jlspline_t [ Types ]

[ Top ] [ m_special_funcs ] [ Types ]

NAME

 jlspline_t

FUNCTION

  Object used to interpolate Bessel functions

SOURCE

65  public :: fermi_dirac        ! Fermi Dirac distribution
66  public :: bose_einstein      ! Bose Einstein distribution
67 
68  type,public :: jlspline_t
69 
70    integer :: nx
71    ! number of points on linear mesh used in spline.
72 
73    integer :: mlang
74    ! mlang= max angular momentum + 1
75 
76    real(dp) :: delta
77    ! Step of linear mesh.
78 
79    real(dp) :: maxarg
80    ! max arg value.
81 
82    real(dp),allocatable :: xx(:)
83    ! xx(nx)
84    ! coordinates of points belonging to the grid
85 
86    real(dp),allocatable :: bess_spl(:,:)
87    ! bess_spl(nx,mlang)
88    ! bessel functions computed on the linear mesh
89 
90    real(dp),allocatable :: bess_spl_der(:,:)
91    ! bess_spl_der(nx,mlang)
92    ! the second derivatives of the cubic spline.
93 
94  end type jlspline_t

m_special_funcs/k_fermi [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

  k_fermi

FUNCTION

  Returns the Fermi wave vector corresponding to the local value of the real space density rhor.

INPUTS

  rhor=Local density in real space.

PARENTS

SOURCE

1620 elemental function k_fermi(rhor)
1621 
1622 
1623 !This section has been created automatically by the script Abilint (TD).
1624 !Do not modify the following lines by hand.
1625 #undef ABI_FUNC
1626 #define ABI_FUNC 'k_fermi'
1627 !End of the abilint section
1628 
1629  implicit none
1630 
1631 !Arguments ------------------------------------
1632 !scalars
1633  real(dp),intent(in) :: rhor
1634  real(dp) :: k_fermi
1635 
1636 !Local variables-------------------------------
1637 !scalars
1638  real(dp),parameter :: pisq=pi**2
1639 
1640 ! *************************************************************************
1641 
1642  k_fermi = (three*pisq*rhor)**third
1643 
1644 end function k_fermi

m_special_funcs/k_thfermi [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

  k_thfermi

FUNCTION

  Returns the Thomas-Fermi wave vector corresponding to the local value of the real space density rhor.

INPUTS

  rhor=Local density in real space.

PARENTS

SOURCE

1663 elemental function k_thfermi(rhor)
1664 
1665 
1666 !This section has been created automatically by the script Abilint (TD).
1667 !Do not modify the following lines by hand.
1668 #undef ABI_FUNC
1669 #define ABI_FUNC 'k_thfermi'
1670 !End of the abilint section
1671 
1672  implicit none
1673 
1674 !Arguments ------------------------------------
1675 !scalars
1676  real(dp),intent(in) :: rhor
1677  real(dp) :: k_thfermi
1678 
1679 !Local variables-------------------------------
1680 !scalars
1681  real(dp),parameter :: pisq=pi**2
1682 
1683 ! *************************************************************************
1684 
1685  k_thfermi = SQRT(four*k_fermi(rhor)*piinv)
1686 
1687 end function k_thfermi

m_special_funcs/laguerre [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 laguerre

FUNCTION

 Laguerre(x,n,a). Returns a (dp) real.

INPUTS

   x position
   n order of laguerre polynomial
   a

OUTPUT

   Laguerre(x,n,a) (dp)

PARENTS

CHILDREN

   factorial

SOURCE

377 function laguerre(x,n,a)
378 
379 
380 !This section has been created automatically by the script Abilint (TD).
381 !Do not modify the following lines by hand.
382 #undef ABI_FUNC
383 #define ABI_FUNC 'laguerre'
384 !End of the abilint section
385 
386  implicit none
387 
388 !Arguments ---------------------------------------------
389 !scalars
390  integer,intent(in),optional :: n,a
391  real(dp)                    :: laguerre
392  real(dp),intent(in)         :: x
393 
394 !Local variables ---------------------------------------
395 !scalars
396  integer :: ii, nn, aa
397 
398 !arrays
399  real(dp),allocatable :: ff(:)
400 
401 ! *********************************************************************
402 
403  if (present(n)) then
404    nn=n
405  else
406    nn=1
407  end if
408 
409  if (present(a)) then
410    aa=a
411  else
412    aa=0
413  end if
414  ABI_ALLOCATE(ff,(nn+1))
415  ff=0.0_dp
416  ff=(/ (binomcoeff(nn+aa,nn-ii)*((-1.0_dp)*x)**ii/factorial(ii) ,ii=0,nn) /)
417  laguerre=sum(ff)
418 
419  ABI_DEALLOCATE(ff)
420 
421 end function laguerre

m_special_funcs/permutations [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 permutations

FUNCTION

 Returns N!/(N-k)!  if N>=0 and N-k>0
                    otherwise 0 is returned
 Output is real

INPUTS

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

OUTPUT

   permutations= n!/(n-k)! (real)

PARENTS

      green_atomic_hubbard

CHILDREN

SOURCE

270 pure function permutations(nn,kk)
271 
272 
273 !This section has been created automatically by the script Abilint (TD).
274 !Do not modify the following lines by hand.
275 #undef ABI_FUNC
276 #define ABI_FUNC 'permutations'
277 !End of the abilint section
278 
279  implicit none
280 
281 !Arguments ---------------------------------------------
282 !scalars
283  integer,intent(in) :: kk,nn
284  real(dp) :: permutations
285 
286 !Local variables ---------------------------------------
287 !scalars
288  integer :: ii
289  real(dp) :: pp
290 
291 ! *********************************************************************
292 
293  if ((nn>=0).and.((nn-kk)>=0)) then
294    pp=one
295    do ii=nn-kk+1,nn
296      pp=pp*ii
297    end do
298  else
299    pp=zero
300  end if
301 
302  permutations=pp
303 
304 end function permutations

m_special_funcs/RadFnH [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 RadFnH

FUNCTION

  RadFnH(r,n,l,Z) radial function of atomic wavefunction with nuclear charge Z.
  for quantum number n, and l.
  Default: Fe 3d function. Returns a (dp) real.

INPUTS

   r radius
   n principal quantum number
   l quantum number

OUTPUT

  RadFnH(r,n,l,Z) (dp)

PARENTS

CHILDREN

  Laguerre
  factorial

SOURCE

453 function RadFnH(r,n,l,Z)
454 
455 
456 !This section has been created automatically by the script Abilint (TD).
457 !Do not modify the following lines by hand.
458 #undef ABI_FUNC
459 #define ABI_FUNC 'RadFnH'
460 !End of the abilint section
461 
462  implicit none
463 
464 !Arguments ---------------------------------------------
465 !scalars
466  integer,intent(in),optional :: n,l
467  real(dp) :: RadFnH
468  real(dp),intent(in) :: r
469  real(dp),intent(in),optional :: Z
470 
471 !Local variables ---------------------------------------
472 !scalars
473  integer   :: nn,ll
474  real(dp)  :: ff,rr,ZZ
475 
476 ! *********************************************************************
477 
478  if (present(n)) then
479    nn=n
480  else
481    nn=3
482  end if
483 
484  if (present(l)) then
485    ll=l
486  else
487    ll=2
488  end if
489 
490  if (present(Z)) then
491    ZZ=Z
492  else
493    ZZ=28.0_dp
494  end if
495 
496  rr=ZZ*r/nn
497  ff=exp(log(ZZ*1.0_dp)*(3.0_dp/2.0_dp))*2/nn**2
498  ff=ff*sqrt(factorial(nn-ll-1)/factorial(nn+ll))*(2*rr)**ll
499  RadFnH=ff*exp(-1*rr)*laguerre(2*rr,nn-ll-1,2*ll+1)
500 
501 end function RadFnH

m_special_funcs/sbf8 [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 sbf8

FUNCTION

 Computes set of spherical bessel functions using accurate algorithm
 based on downward recursion in order and normalization sum.
 Power series used at small arguments.

INPUTS

  nm=maximum angular momentum wanted + one
  xx=argument of sbf

OUTPUT

  sb_out(nm)=values of spherical bessel functions for l=0,nm-1

PARENTS

      posdoppler,psp8nl,qijb_kk,qmc_prep_ctqmc,smatrix_pawinit

CHILDREN

SOURCE

1423 subroutine sbf8(nm,xx,sb_out)
1424 
1425 
1426 !This section has been created automatically by the script Abilint (TD).
1427 !Do not modify the following lines by hand.
1428 #undef ABI_FUNC
1429 #define ABI_FUNC 'sbf8'
1430 !End of the abilint section
1431 
1432  implicit none
1433 
1434 !Arguments----------------------------------------------------------
1435 !scalars
1436  integer,intent(in) :: nm
1437  real(dp),intent(in) :: xx
1438 !arrays
1439  real(dp),intent(out) :: sb_out(nm)
1440 
1441 !Local variables-------------------------------
1442 !scalars
1443  integer :: nlim,nn
1444  real(dp) :: fn,sn,xi,xn,xs
1445 !arrays
1446  real(dp),allocatable :: sb(:)
1447 
1448 ! *************************************************************************
1449 
1450  if(xx<= 1.0e-36_dp) then
1451 !  zero argument section
1452    sb_out(:)=zero
1453    sb_out(1)=one
1454  else if(xx<1.e-3_dp) then
1455 !  small argument section
1456    xn=one
1457    xs=half*xx**2
1458    do nn=1,nm
1459      sb_out(nn)=xn*(one - xs*(one - xs/(4*nn+6))/(2*nn+1))
1460      xn=xx*xn/(2*nn+1)
1461    end do
1462  else
1463 !  recursion method
1464    if(xx<one) then
1465      nlim=nm+int(15.0e0_dp*xx)+1
1466    else
1467      nlim=nm+int(1.36e0_dp*xx)+15
1468    end if
1469    ABI_ALLOCATE(sb,(nlim+1))
1470    nn=nlim
1471    xi=one/xx
1472    sb(nn+1)=zero
1473    sb(nn)=1.e-18_dp
1474    sn=dble(2*nn-1)*1.e-36_dp
1475    do nn=nlim-1,1,-1
1476      sb(nn)=dble(2*nn+1)*xi*sb(nn+1) - sb(nn+2)
1477    end do
1478    do nn=1,nlim-1
1479      sn=sn + dble(2*nn-1)*sb(nn)*sb(nn)
1480    end do
1481    fn=1.d0/sqrt(sn)
1482    sb_out(:)=fn*sb(1:nm)
1483    ABI_DEALLOCATE(sb)
1484  end if
1485 
1486 end subroutine sbf8