TABLE OF CONTENTS


ABINIT/GAMMA_FUNCTION [ Functions ]

[ Top ] [ Functions ]

NAME

  GAMMA_FUNCTION

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

SOURCE

872 subroutine GAMMA_FUNCTION(X,GA)
873 
874 !       ====================================================
875 !       Purpose: This program computes the gamma function
876 !                Gamma(x) using subroutine GAMMA
877 !       Examples:
878 !                   x          Gamma(x)
879 !                ----------------------------
880 !                  1/3       2.678938534708
881 !                  0.5       1.772453850906
882 !                 -0.5      -3.544907701811
883 !                 -1.5       2.363271801207
884 !                  5.0      24.000000000000
885 !       ====================================================
886 !
887 !  This routine was downloaded from UIUC:
888 !  http://jin.ece.uiuc.edu/routines/routines.html
889 !
890 !  The programs appear to accompany a book "Computation of Special
891 !  Functions" (1996) John Wiley and Sons, but are distributed online
892 !  by the authors. Exact copyright should be checked.
893 !
894 !  Authors / copyright:
895 !     Shanjie Zhang and Jianming Jin
896 !     Proposed contact is:  j-jin1@uiuc.edu
897 !
898 !  20 October 2008:
899 !     Incorporated into ABINIT by M. Verstraete
900 !
901 !
902 !
903 !       ==================================================
904 !       Purpose: Compute the gamma function Gamma(x)
905 !       Input :  x  --- Argument of Gamma(x)
906 !                       ( x is not equal to 0,-1,-2, etc )
907 !       Output:  GA --- Gamma(x)
908 !       ==================================================
909 !
910 
911   ! arguments
912 
913   real(dp),intent(in) :: x
914   real(dp),intent(out) :: ga
915 
916   ! local variables
917   integer :: k,m
918   real(dp) :: m1,z,r,gr
919   real(dp) :: G(26)
920 
921   ! source code:
922 
923   ! initialization of reference data
924   G=(/1.0D0,0.5772156649015329D0, &
925      &  -0.6558780715202538D0, -0.420026350340952D-1, &
926      &   0.1665386113822915D0,-.421977345555443D-1, &
927      &  -.96219715278770D-2, .72189432466630D-2, &
928      &  -.11651675918591D-2, -.2152416741149D-3, &
929      &   .1280502823882D-3, -.201348547807D-4, &
930      &  -.12504934821D-5, .11330272320D-5, &
931      &  -.2056338417D-6, .61160950D-8, &
932      &   .50020075D-8, -.11812746D-8, &
933      &   .1043427D-9, .77823D-11, &
934      &  -.36968D-11, .51D-12, &
935      &  -.206D-13, -.54D-14, .14D-14, .1D-15/)
936 
937 
938   ! for the integer case, do explicit factorial
939   if (X==int(X)) then
940     if (X > 0.0D0) then
941       GA=1.0D0
942       M1=X-1
943       do K=2,int(M1)
944         GA=GA*K
945       end do
946     else
947       GA=1.0D+300
948     end if
949   ! for the integer case, do explicit factorial
950   else
951     if (abs(X) > 1.0D0) then
952       Z=abs(X)
953       M=int(Z)
954       R=1.0D0
955       do K=1,M
956         R=R*(Z-K)
957       end do
958       Z=Z-M
959     else
960       Z=X
961     end if
962     GR=G(26)
963     do K=25,1,-1
964       GR=GR*Z+G(K)
965     end do
966     GA=1.0D0/(GR*Z)
967     if (abs(X) > 1.0D0) then
968       GA=GA*R
969       if (X < 0.0D0) GA=-PI/(X*GA*SIN(PI*X))
970     end if
971   end if
972   return
973 
974 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-2022 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

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

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

SOURCE

591 elemental function abi_derf(yy) result(derf_yy)
592 
593 !Arguments ------------------------------------
594 !scalars
595  real(dp),intent(in) :: yy
596  real(dp) :: derf_yy
597 
598 !Local variables-------------------------------
599  integer          ::  done,ii,isw
600 ! coefficients for 0.0 <= yy < .477
601  real(dp), parameter :: &
602 &  pp(5)=(/ 113.8641541510502e0_dp, 377.4852376853020e0_dp,  &
603 &           3209.377589138469e0_dp, .1857777061846032e0_dp,  &
604 &           3.161123743870566e0_dp /)
605  real(dp), parameter :: &
606 &  qq(4)=(/ 244.0246379344442e0_dp, 1282.616526077372e0_dp,  &
607 &           2844.236833439171e0_dp, 23.60129095234412e0_dp/)
608 ! coefficients for .477 <= yy <= 4.0
609  real(dp), parameter :: &
610 &  p1(9)=(/ 8.883149794388376e0_dp, 66.11919063714163e0_dp,  &
611 &           298.6351381974001e0_dp, 881.9522212417691e0_dp,  &
612 &           1712.047612634071e0_dp, 2051.078377826071e0_dp,  &
613 &           1230.339354797997e0_dp, 2.153115354744038e-8_dp, &
614 &           .5641884969886701e0_dp /)
615  real(dp), parameter :: &
616 &  q1(8)=(/ 117.6939508913125e0_dp, 537.1811018620099e0_dp,  &
617 &           1621.389574566690e0_dp, 3290.799235733460e0_dp,  &
618 &           4362.619090143247e0_dp, 3439.367674143722e0_dp,  &
619 &           1230.339354803749e0_dp, 15.74492611070983e0_dp/)
620   ! coefficients for 4.0 < y,
621  real(dp), parameter :: &
622 &  p2(6)=(/ -3.603448999498044e-01_dp, -1.257817261112292e-01_dp,   &
623 &           -1.608378514874228e-02_dp, -6.587491615298378e-04_dp,   &
624 &           -1.631538713730210e-02_dp, -3.053266349612323e-01_dp/)
625  real(dp), parameter :: &
626 &  q2(5)=(/ 1.872952849923460e0_dp   , 5.279051029514284e-01_dp,    &
627 &           6.051834131244132e-02_dp , 2.335204976268692e-03_dp,    &
628 &           2.568520192289822e0_dp /)
629  real(dp), parameter :: &
630 &  sqrpi=.5641895835477563e0_dp, xbig=13.3e0_dp, xlarge=6.375e0_dp, xmin=1.0e-10_dp
631  real(dp) ::  res,xden,xi,xnum,xsq,xx
632 
633 ! ******************************************************************
634 
635  xx = yy
636  isw = 1
637 !Here change the sign of xx, and keep track of it thanks to isw
638  if (xx<0.0e0_dp) then
639    isw = -1
640    xx = -xx
641  end if
642 
643  done=0
644 
645 !Residual value, if yy < -6.375e0_dp
646  res=-1.0e0_dp
647 
648 !abs(yy) < .477, evaluate approximation for erfc
649  if (xx<0.477e0_dp) then
650 !  xmin is a very small number
651    if (xx<xmin) then
652      res = xx*pp(3)/qq(3)
653    else
654      xsq = xx*xx
655      xnum = pp(4)*xsq+pp(5)
656      xden = xsq+qq(4)
657      do ii = 1,3
658        xnum = xnum*xsq+pp(ii)
659        xden = xden*xsq+qq(ii)
660      end do
661      res = xx*xnum/xden
662    end if
663    if (isw==-1) res = -res
664    done=1
665  end if
666 
667 !.477 < abs(yy) < 4.0 , evaluate approximation for erfc
668  if (xx<=4.0e0_dp .and. done==0 ) then
669    xsq = xx*xx
670    xnum = p1(8)*xx+p1(9)
671    xden = xx+q1(8)
672    do ii=1,7
673      xnum = xnum*xx+p1(ii)
674      xden = xden*xx+q1(ii)
675    end do
676    res = xnum/xden
677    res = res* exp(-xsq)
678    if (isw.eq.-1) then
679      res = res-1.0e0_dp
680    else
681      res=1.0e0_dp-res
682    end if
683    done=1
684  end if
685 
686 !y > 13.3e0_dp
687  if (isw > 0 .and. xx > xbig .and. done==0 ) then
688    res = 1.0e0_dp
689    done=1
690  end if
691 
692 !4.0 < yy < 13.3e0_dp  .or. -6.375e0_dp < yy < -4.0
693 !evaluate minimax approximation for erfc
694  if ( ( isw > 0 .or. xx < xlarge ) .and. done==0 ) then
695    xsq = xx*xx
696    xi = 1.0e0_dp/xsq
697    xnum= p2(5)*xi+p2(6)
698    xden = xi+q2(5)
699    do ii = 1,4
700      xnum = xnum*xi+p2(ii)
701      xden = xden*xi+q2(ii)
702    end do
703    res = (sqrpi+xi*xnum/xden)/xx
704    res = res* exp(-xsq)
705    if (isw.eq.-1) then
706      res = res-1.0e0_dp
707    else
708      res=1.0e0_dp-res
709    end if
710  end if
711 
712 !All cases have been investigated
713  derf_yy = res
714 
715 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

SOURCE

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

SOURCE

1007 subroutine besjm(arg,besjx,cosx,nn,nx,sinx,xx)
1008 
1009 !Arguments ------------------------------------
1010 !scalars
1011  integer,intent(in) :: nn,nx
1012  real(dp),intent(in) :: arg
1013 !arrays
1014  real(dp),intent(in) :: cosx(nx),sinx(nx),xx(nx)
1015  real(dp),intent(out) :: besjx(nx)
1016 
1017 !Local variables-------------------------------
1018 !scalars
1019  integer :: ix,switchx
1020 !Series or rational polynomial coefficients
1021  real(dp),parameter :: b01=1.d0/6.d0,b02=1.d0/120.d0,b03=1.d0/5040.d0
1022  real(dp),parameter :: b04=1.d0/362880.d0,b11=0.8331251468724171d-1
1023  real(dp),parameter :: b12=0.2036961284395412d-2,b13=0.1932970379901801d-4
1024  real(dp),parameter :: b14=0.6526053169009489d-7,b21=0.5867824627555163d-1
1025  real(dp),parameter :: b22=0.1152501878595934d-2,b23=0.1011071389414764d-4
1026  real(dp),parameter :: b24=0.4172322111421287d-7,b25=0.6790616688656543d-10
1027  real(dp),parameter :: b31=0.439131885807176d-1,b32=0.6813139609887099d-3
1028  real(dp),parameter :: b33=0.4899103784264755d-5,b34=0.17025590795625d-7
1029  real(dp),parameter :: b35=0.2382642910613347d-10,b41=0.3587477991030971d-1
1030  real(dp),parameter :: b42=0.4833719855268907d-3,b43=0.3238388977796242d-5
1031  real(dp),parameter :: b44=0.1171802513125112d-7,b45=0.223261650431992d-10
1032  real(dp),parameter :: b46=.1800045587335951d-13,b51=0.295232406376567d-1
1033  real(dp),parameter :: b52=0.3359864457080573d-3,b53=0.19394750603618d-5
1034  real(dp),parameter :: b54=0.6143166228216219d-8,b55=0.10378501636108d-10
1035  real(dp),parameter :: b56=.749975122872713d-14
1036  real(dp),parameter :: c11=0.1668748531275829d-1,c12=0.1342812442426702d-3
1037  real(dp),parameter :: c13=0.6378249315355233d-6,c14=0.1573564527360138d-8
1038  real(dp),parameter :: c21=0.127503251530198d-1,c22=0.7911240539893565d-4
1039  real(dp),parameter :: c23=0.3044380758068054d-6,c24=0.7439837832363479d-9
1040  real(dp),parameter :: c25=0.9515065658793124d-12,c31=0.1164236697483795d-1
1041  real(dp),parameter :: c32=0.654858636312224d-4,c33=0.2265576367562734d-6
1042  real(dp),parameter :: c34=0.4929905563217352d-9,c35=0.555120465710914d-12
1043  real(dp),parameter :: c41=0.9579765544235745d-2,c42=0.4468999977536864d-4
1044  real(dp),parameter :: c43=0.1315634305905896d-6,c44=0.2615492488301639d-9
1045  real(dp),parameter :: c45=0.3387473312408129d-12,c46=.2280866204624012d-15
1046  real(dp),parameter :: c51=0.8938297823881763d-2,c52=0.3874149021633025d-4
1047  real(dp),parameter :: c53=0.1054692715135225d-6,c54=0.192879620987602d-9
1048  real(dp),parameter :: c55=0.2284469423833734d-12,c56=0.139729234332572d-15
1049  real(dp),parameter :: ffnth=1.d0/15.d0,o10395=1.d0/10395d0,oo105=1.d0/105.d0
1050  real(dp),parameter :: oo945=1.d0/945.d0
1051  real(dp) :: bot,rr,rsq,top
1052  character(len=500) :: message
1053 
1054 ! *************************************************************************
1055 
1056  if (nn==0) then
1057 
1058    switchx=nx+1
1059    do ix=1,nx
1060      rr=arg*xx(ix)
1061      if (rr<=1.d-1) then
1062        rsq=rr*rr
1063        besjx(ix)=1.d0-rsq*(b01-rsq*(b02-rsq*(b03-rsq*b04)))
1064      else
1065        switchx=ix
1066        exit
1067      end if
1068    end do
1069 
1070    do ix=switchx,nx
1071      rr=arg*xx(ix)
1072      besjx(ix)=sinx(ix)/rr
1073    end do
1074 
1075  else if (nn==1) then
1076 
1077    switchx=nx+1
1078    do ix=1,nx
1079      rr=arg*xx(ix)
1080      if (rr<=1.d0) then
1081        rsq=rr*rr
1082        top=1.d0-rsq*(b11-rsq*(b12-rsq*(b13-rsq*b14)))
1083        bot=1.d0+rsq*(c11+rsq*(c12+rsq*(c13+rsq*c14)))
1084        besjx(ix)=third*rr*top/bot
1085      else
1086        switchx=ix
1087        exit
1088      end if
1089    end do
1090 
1091    do ix=switchx,nx
1092      rr=arg*xx(ix)
1093      rsq=rr*rr
1094      besjx(ix)=(sinx(ix)-rr*cosx(ix))/rsq
1095    end do
1096 
1097  else if (nn==2) then
1098 
1099    switchx=nx+1
1100    do ix=1,nx
1101      rr=arg*xx(ix)
1102      if (rr<=2.d0) then
1103        rsq=rr*rr
1104        top=1.d0-rsq*(b21-rsq*(b22-rsq*(b23-rsq*(b24-rsq*b25))))
1105        bot=1.d0+rsq*(c21+rsq*(c22+rsq*(c23+rsq*(c24+rsq*c25))))
1106        besjx(ix)=ffnth*rsq*top/bot
1107      else
1108        switchx=ix
1109        exit
1110      end if
1111    end do
1112 
1113    do ix=switchx,nx
1114      rr=arg*xx(ix)
1115      rsq=rr*rr
1116      besjx(ix)=((3.d0-rsq)*sinx(ix)-3.d0*rr*cosx(ix))/(rr*rsq)
1117    end do
1118 
1119  else if (nn==3) then
1120 
1121    switchx=nx+1
1122    do ix=1,nx
1123      rr=arg*xx(ix)
1124      if (rr<=2.d0) then
1125        rsq=rr*rr
1126        top=1.d0-rsq*(b31-rsq*(b32-rsq*(b33-rsq*(b34-rsq*b35))))
1127        bot=1.d0+rsq*(c31+rsq*(c32+rsq*(c33+rsq*(c34+rsq*c35))))
1128        besjx(ix)=rr*rsq*oo105*top/bot
1129      else
1130        switchx=ix
1131        exit
1132      end if
1133    end do
1134 
1135    do ix=switchx,nx
1136      rr=arg*xx(ix)
1137      rsq=rr*rr
1138      besjx(ix)=( (15.d0-6.d0*rsq)*sinx(ix)&
1139 &     + rr*(rsq-15.d0)  *cosx(ix) ) /(rsq*rsq)
1140    end do
1141 
1142  else if (nn==4) then
1143 
1144    switchx=nx+1
1145    do ix=1,nx
1146      rr=arg*xx(ix)
1147      if (rr<=4.d0) then
1148        rsq=rr*rr
1149        top=1.d0-rsq*(b41-rsq*(b42-rsq*(b43-rsq*(b44-rsq*(b45-rsq*b46)))))
1150        bot=1.d0+rsq*(c41+rsq*(c42+rsq*(c43+rsq*(c44+rsq*(c45+rsq*c46)))))
1151        besjx(ix)=rsq*rsq*oo945*top/bot
1152      else
1153        switchx=ix
1154        exit
1155      end if
1156    end do
1157 
1158    do ix=switchx,nx
1159      rr=arg*xx(ix)
1160      rsq=rr*rr
1161      besjx(ix)=( (105.d0-rsq*(45.d0-rsq)) *sinx(ix)&
1162 &     + rr * (10.d0*rsq-105.d0)  *cosx(ix) ) /(rsq*rsq*rr)
1163    end do
1164 
1165  else if (nn==5) then
1166 
1167    switchx=nx+1
1168    do ix=1,nx
1169      rr=arg*xx(ix)
1170      if (rr<=4.d0) then
1171        rsq=rr*rr
1172        top=1.d0-rsq*(b51-rsq*(b52-rsq*(b53-rsq*(b54-rsq*(b55-rsq*b56)))))
1173        bot=1.d0+rsq*(c51+rsq*(c52+rsq*(c53+rsq*(c54+rsq*(c55+rsq*c56)))))
1174        besjx(ix)=rsq*rsq*rr*o10395*top/bot
1175      else
1176        switchx=ix
1177        exit
1178      end if
1179    end do
1180 
1181    do ix=switchx,nx
1182      rr=arg*xx(ix)
1183      rsq=rr*rr
1184      besjx(ix)=( (945.d0-rsq*(420.d0-rsq*15.d0)) *sinx(ix)&
1185 &     + rr * (945.d0-rsq*(105.d0-rsq))  *cosx(ix) ) /(rsq*rsq*rr)
1186    end do
1187 
1188  else
1189    write(message, '(a,i0,a)' )' besjm only defined for nn in [0,5]; input was nn=',nn,'.'
1190    ABI_BUG(message)
1191  end if
1192 
1193 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)

SOURCE

289 elemental function binomcoeff(n,k)
290 
291 !Arguments ---------------------------------------------
292 !scalars
293  integer,intent(in) :: n,k
294  real(dp) :: binomcoeff
295 
296 ! *********************************************************************
297 
298  binomcoeff=factorial(n)/(factorial(k)*factorial(n-k))
299 
300 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

SOURCE

1332 function bose_einstein(energy, temperature)
1333 
1334 !Arguments ------------------------------------
1335 !scalars
1336  real(dp),intent(in) :: energy, temperature
1337  real(dp) :: bose_einstein
1338 
1339 !Local variables-------------------------------
1340 !scalars
1341  real(dp) :: arg
1342  character(len=500) :: message
1343 
1344 ! *************************************************************************
1345 
1346  bose_einstein = zero
1347  if (temperature > tol12) then
1348    arg = energy/temperature
1349    if(arg > tol12 .and. arg < 600._dp)then
1350      bose_einstein = one / (exp(arg)  - one)
1351    else if (arg < tol12) then
1352      write(message,'(a)') 'No Bose Einstein for negative energies'
1353      ABI_WARNING(message)
1354    end if
1355  else
1356    write(message,'(a)') 'No Bose Einstein for negative or 0 T'
1357    ABI_WARNING(message)
1358  end if
1359 
1360 
1361 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

SOURCE

166 pure function clp(x)
167 
168 !Arguments ------------------------------------
169 !scalars
170  real(dp) :: clp
171  real(dp),intent(in) :: x
172 
173 ! **********************************************************************
174 
175  if(x > half) then
176    clp=x-one
177  elseif(x < -half) then
178    clp=x+one
179  else
180    clp=x
181  end if
182 
183 end function clp

m_special_funcs/dip12 [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

  dip12

FUNCTION

  Returns the complete Fermi integral of order 1/2.
  Based on an analytical approximation.

INPUTS

  gamma=complete Fermi integral argument

OUTPUT

  dip12=resulting function

SOURCE

1381 function dip12(gamma)
1382 ! Arguments -------------------------------
1383 ! Scalars
1384  real(dp),intent(in) :: gamma
1385 
1386 ! Local variables -------------------------
1387 ! Scalars
1388  real(dp) :: d,dip12,dy
1389 
1390 ! *********************************************************************
1391 
1392  if (gamma.lt.3.) then
1393    dy=exp(gamma)
1394    if (gamma+1.9375.LE.0) then
1395      dip12=dy*&
1396      & (1.-dy*(0.35355283-dy*(0.19242767-dy*(0.12456909-dy*&
1397      & (8.5114507E-02-dy*4.551794E-02)))))
1398    else
1399      d=gamma-0.5
1400      dip12=dy*(0.677695804-d*(0.187773135+d*(2.16197521E-02-d*&
1401      & (9.23703807E-03+d*&
1402      & (1.71735167E-03-d*(6.07913775E-04+d*&
1403      & (1.1448629E-04-d*&
1404      & (4.544432E-05+d*(6.4719368E-06-d*(3.794983E-06+d*&
1405      & (1.7338029E-07-d*&
1406      & (3.5546516E-07-d*(3.7329191E-08+d*&
1407      & (3.3097822E-08-d*&
1408      & (8.3190193E-09+d*(2.2752769E-09-d*(7.836005E-10+d*&
1409      & (7.519551E-11-d*2.960006E-11))))))))))))))))))
1410    end if
1411  else if (gamma.lt.20.) then
1412    if (gamma.lt.10.) then
1413      d=gamma-6.5
1414      dip12=12.839811+d*&
1415      & (2.844774+d*(0.114920926-d*(3.43733039E-03-d*&
1416      & (2.3980356E-04-d*&
1417      & (2.0201888E-05-d*(1.5219883E-06-d*&
1418      & (6.2770524E-08+d*&
1419      & (4.8830336E-09-d*(2.1031164E-09-d*(5.785753E-10-d*&
1420      & (7.233066E-11-d*1.230727E-12)))))))))))
1421    else
1422      d=gamma-14.5
1423      dip12=41.7799227+d*&
1424      & (4.2881461+d*(7.45407825E-02-d*(8.79243296E-04-d*&
1425      & (2.38288861E-05-d*&
1426      & (8.82474867E-07-d*(3.82865217E-08-d*&
1427      & (1.9274292E-09-d*(1.42248669E-10-d*8.17019813E-12))))))))
1428    end if
1429  else
1430    d=1./gamma
1431    dy=gamma*sqrt(gamma)/1.329340388
1432    dip12=dy*(1.-d*(9.354E-07-d*(1.2338391-d*(6.77931E-03-d*&
1433    & 1.17871643))))
1434  end if
1435  dip12=dip12*0.88622692
1436 end function dip12

m_special_funcs/dip32 [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

  dip32

FUNCTION

  Returns the complete Fermi integral of order 3/2.
  Based on an analytical approximation.

INPUTS

  gamma=complete Fermi integral argument

OUTPUT

  dip32=resulting function

SOURCE

1456 function dip32(gamma)
1457 ! Arguments -------------------------------
1458 ! Scalars
1459  real(dp),intent(in) :: gamma
1460 
1461 ! Local variables -------------------------
1462 ! Scalars
1463  real(dp) :: d,dip32,dval
1464 
1465 ! *********************************************************************
1466 
1467  if (gamma.GT.1.75) then
1468    dval=gamma*gamma*SQRT(gamma)
1469    if (gamma.LT.4.5) then
1470      d=gamma-3.125
1471      dip32=(1.27623+0.596065*gamma+0.3*dval)*&
1472      & (1.0055558385-d*(5.23889494E-03+d*&
1473      & (3.13523144E-03-d*(3.06124286E-03-d*&
1474      & (1.3644667E-03-d*&
1475      & (4.1528384E-04-d*(8.901188E-05-d*(1.079979E-05+d*&
1476      & (2.29058E-06-d*(2.58985E-06-d*7.30909E-07))))))))))
1477    else if (gamma.LT.12.) then
1478      if (gamma.LT.8.) then
1479        d=gamma-6.25
1480        dip32=(2.01508+0.425775*gamma+0.3*dval)*&
1481        & (1.000387131-d*(3.93626295E-04+d*&
1482        & (2.55710115E-04-d*&
1483        & (1.57383494E-04-d*(5.0286036E-05-d*&
1484        & (1.2073559865E-05-d*&
1485        & (2.4909523213E-06-d*(5.244328548E-07-d*&
1486        & 8.0884033896E-08))))))))
1487      else
1488        d=gamma-10.
1489        dip32=0.3*dval*&
1490        & (1.064687247-d*(1.22972303E-02-d*(1.8362121E-03-d*&
1491        & (2.433558E-04-d*(3.018186E-05-d*(3.5694E-06-d*&
1492        & (4.11212E-07-d*(5.2151E-08-d*5.8424E-09))))))))
1493      end if
1494    else
1495      d=1./gamma
1496      dip32=0.30090111127*dval*&
1497      & (1.-d*(2.863E-06-d*(6.168876549-d*&
1498      & (1.740553E-02+d*(1.425257+d*2.95887)))))
1499    end if
1500  else if (gamma+0.75.LE.0) then
1501    d=EXP(gamma)
1502    dip32=d*&
1503    & (1.-d*(1.76775246E-01-d*(6.4124584E-02-d*&
1504    & (3.1027055E-02-d*(1.6797637E-02-d*&
1505    & (8.212636E-03-d*(2.384106E-03)))))))
1506  else
1507    d=gamma-0.5
1508    dip32=EXP(gamma)*(0.846691-0.128948*gamma)*&
1509    & (1.034064158+d*(2.778947E-03-d*&
1510    & (3.572502805E-02+d*(3.0411645E-03-d*&
1511    & (1.7380548E-03+d*(2.7756776E-04-d*&
1512    & (8.08302E-05+d*(1.59606E-05-d*&
1513    & (3.8144E-06+d*7.4446E-07)))))))))
1514  end if
1515  dip32=dip32*1.32934038
1516 end function dip32

m_special_funcs/djp12 [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

  djp12

FUNCTION

  Returns the incomplete Fermi integral of order 1/2.
  Based on an analytical approximation.

INPUTS

  xcut=lower bound of the incomplete Fermi integral
  gamma=incomplete Fermi integral argument

OUTPUT

  djp12=resulting function

SOURCE

1537 function djp12(xcut,gamma)
1538 ! Arguments -------------------------------
1539 ! Scalars
1540  real(dp),intent(in) :: xcut,gamma
1541  real(dp) :: djp12
1542 
1543 ! Local variables -------------------------
1544 ! Scalars
1545  real(dp) :: d2h,db,dc,dd,de,df1,df2,df3,dh
1546  real(dp) :: ds,dt,dv,dw,dxm,dxp
1547  integer :: i,ind,iq,k,nm,np,nq
1548 ! Arrays
1549  real(dp) :: dq(5),df(101),dy(101)
1550 
1551 ! *********************************************************************
1552 
1553  dh=0.2D+0
1554  d2h=0.4D+0
1555  nm=101
1556  ind=0
1557  dq=(/1.D+0,2.828427124D+0,5.196152423D+0,8.D+0,1.118033989D+1/)
1558 
1559  djp12=0.D+0
1560  dxm=gamma-1.5D+1
1561  if (xcut.gt.dxm) then
1562    if (ind.eq.0) then
1563      do i=1,nm
1564        dy(i)=-1.5D+1+(i-1)*dh
1565        df(i)=1.D+0+dexp(dy(i))
1566      end do
1567      ind=1
1568    end if
1569    dxp=gamma+5.D+0
1570    if (xcut.lt.dxp) then
1571      dc=dxp
1572    else
1573      dc=xcut
1574    end if
1575    db=dexp(gamma-dc)
1576    dt=db
1577    do iq=1,5
1578      dd=iq*dc
1579      ds=dsqrt(dd)
1580      dw=1.+.3275911*ds
1581      dw=1.D+0/dw
1582      dv=dw*(.2258368458D+0+&
1583      & dw*(-.2521286676D+0+dw*(1.2596951294D+0+&
1584      & dw*(-1.2878224530D+0+dw*(.9406460699D+0)))))
1585      dv=dv+ds
1586      de=dt*dv/dq(iq)
1587      djp12=djp12+de
1588      if (dabs(de).lt.(1.D-07*djp12)) exit
1589      dt=-dt*db
1590    end do
1591    if (xcut.ge.dxp) return
1592    np=(dxp-xcut)/dh
1593    np=2*(np/2)
1594    np=nm-np
1595    nq=(15.-gamma)/dh
1596    nq=1+2*(nq/2)
1597    if (np.lt.nq) np=nq
1598    if (np.le.nm) then
1599      df3=0.D+0
1600      dt=dy(np)+gamma
1601      dv=(dt-xcut)/2.D+0
1602      df3=0.D0
1603      if (dt.ge.1.D-13) df3=dsqrt(dt)
1604      if (dabs(dv).ge.1.D-13) then
1605        df1=dsqrt(xcut)
1606        dt=df1+df3
1607        dw=(dv+dv)/(dt*dt)
1608        df2=dw*dw
1609        df2=df2+df2
1610        db=df2*(df2+7.D+0)+7.D+1
1611        dc=7.D+0*(1.D+1-df2)
1612        dc=dc*dw
1613        dd=-df2*(df2-2.8D+1)+1.4D+2
1614        dd=dd+dd
1615        ds=dt*((db-dc)/(1.D+0+dexp(xcut-gamma))+&
1616        & dd/(1.D+0+dexp(xcut+dv-gamma))+(db+dc)/df(np))
1617        ds=ds*dv/4.2D+2
1618        djp12=djp12+ds
1619      end if
1620      if (np.ne.nm) then
1621        ds=0.D+0
1622        np=np+2
1623        do k=np,nm,2
1624          df1=df3
1625          df3=dsqrt(dy(k)+gamma)
1626          dt=df1+df3
1627          dw=d2h/(dt*dt)
1628          df2=dw*dw
1629          df2=df2+df2
1630          db=df2*(df2+7.D+0)+7.D+1
1631          dc=7.D+0*(1.D+1-df2)
1632          dc=dc*dw
1633          dd=-df2*(df2-2.8D+1)+1.4D+2
1634          dd=dd+dd
1635          ds=ds+dt*((db-dc)/df(k-2)+dd/df(k-1)+(db+dc)/df(k))
1636        end do
1637        ds=ds*dh/4.2D+2
1638        djp12=djp12+ds
1639      end if
1640      if (xcut.ge.dxm) return
1641    end if
1642  end if
1643  djp12=dip12(gamma)-xcut*dsqrt(xcut)/1.5D+0
1644 end function djp12

m_special_funcs/djp32 [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

  djp32

FUNCTION

  Returns the incomplete Fermi integral of order 3/2.
  Based on an analytical approximation.

INPUTS

  xcut=lower bound of the incomplete Fermi integral
  gamma=incomplete Fermi integral argument

OUTPUT

  djp32=resulting function

SOURCE

1665 function djp32(xcut,gamma)
1666 ! Arguments -------------------------------
1667 ! Scalars
1668  real(dp),intent(in) :: xcut,gamma
1669  real(dp) :: djp32
1670 
1671 ! Local variables -------------------------
1672 ! Scalars
1673  real(dp) :: d2h,db,dc,dd,de,df1,df2,df3,dh
1674  real(dp) :: ds,dt,dv,dw,dx1,dx2
1675  real(dp) :: dx3,dxm,dxp
1676  integer :: i,ind,iq,k,nm,np,nq
1677 ! Arrays
1678  real(dp) :: dq(5),df(101),dy(101)
1679 
1680 ! *********************************************************************
1681 
1682  dh=0.2D+0
1683  d2h=0.4D+0
1684  nm=101
1685  ind=0
1686  dq=(/1.D+0,5.656854228D+0,1.558845727D+1,3.2D+1,5.590169945D+0/)
1687 
1688  djp32=0.D+0
1689  dxm=gamma-1.5D+1
1690  if (xcut.GT.dxm) then
1691    if (ind.EQ.0) then
1692      do i=1,nm
1693        dy(i)=-1.5D+1+(i-1)*dh
1694        df(i)=1.D+0+DEXP(dy(i))
1695      end do
1696      ind=1
1697    end if
1698    dxp=gamma+5.D+0
1699    if (xcut.LT.dxp) then
1700      dc=dxp
1701    else
1702      dc=xcut
1703    end if
1704    db=DEXP(gamma-dc)
1705    dt=db
1706    do iq=1,5
1707      dd=iq*dc
1708      ds=DSQRT(dd)
1709      dw=1.+.3275911*ds
1710      dw=1.D+0/dw
1711      dv=dw*(.2258368458D+0+&
1712      & dw*(-.2521286676D+0+dw*(1.2596951294D+0+&
1713      & dw*(-1.2878224530D+0+dw*(.9406460699D+0)))))
1714      dv=dv+ds
1715      dv=1.5D+0*dv+ds*dd
1716      de=dt*dv/dq(iq)
1717      djp32=djp32+de
1718      if (DABS(de).LT.(1.D-07*djp32)) exit
1719      dt=-dt*db
1720    end do
1721    if (xcut.GE.dxp) return
1722    np=(dxp-xcut)/dh
1723    np=2*(np/2)
1724    np=nm-np
1725    nq=(15.-gamma)/dh
1726    nq=1+2*(nq/2)
1727    if (np.LT.nq) np=nq
1728    if (np.LE.nm) then
1729      df3=0.D+0
1730      dt=dy(np)+gamma
1731      dv=(dt-xcut)/2.D+0
1732      df3=DSQRT(dt)
1733      dx3=dt
1734      if (DABS(dv).GE.1.D-13) then
1735        df1=DSQRT(xcut)
1736        dt=df1+df3
1737        dw=(dv+dv)/(dt*dt)
1738        df2=dw*dw
1739        df2=df2+df2
1740        db=df2*(df2+7.D+0)+7.D+1
1741        dc=7.D+0*(1.D+1-df2)
1742        dc=dc*dw
1743        dd=-df2*(df2-2.8D+1)+1.4D+2
1744        dd=dd+dd
1745        ds=dt*((db-dc)*xcut/(1.D+0+DEXP(xcut-gamma))+dd*(xcut+dv)&
1746        & /(1.D+0+DEXP(xcut+dv-gamma))+(db+dc)*(dy(np)+gamma)/df(np))
1747        ds=ds*dv/4.2D+2
1748        djp32=djp32+ds
1749      end if
1750      if (np.NE.nm) then
1751        ds=0.D+0
1752        np=np+2
1753        do k=np,nm,2
1754          dx1=dx3
1755          df1=df3
1756          dx2=dy(k-1)+gamma
1757          dx3=dy(k)+gamma
1758          df3=DSQRT(dx3)
1759          dt=df1+df3
1760          dw=d2h/(dt*dt)
1761          df2=dw*dw
1762          df2=df2+df2
1763          db=df2*(df2+7.D+0)+7.D+1
1764          dc=7.D+0*(1.D+1-df2)
1765          dc=dc*dw
1766          dd=-df2*(df2-2.8D+1)+1.4D+2
1767          dd=dd+dd
1768          ds=ds+dt*((db-dc)*dx1/df(k-2)+dd*dx2/df(k-1)&
1769          & +(db+dc)*dx3/df(k))
1770        end do
1771        ds=ds*dh/4.2D+2
1772        djp32=djp32+ds
1773      end if
1774      if (xcut.GE.dxm) return
1775    end if
1776  end if
1777  djp32=dip32(gamma)-xcut*xcut*DSQRT(xcut)/2.5D+0
1778 end function djp32

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)

SOURCE

201 elemental function factorial(nn)
202 
203 !Arguments ---------------------------------------------
204 !scalars
205  integer,intent(in) :: nn
206  real(dp) :: factorial
207 
208 !Local variables ---------------------------------------
209 !scalars
210  integer :: ii
211  real(dp) :: ff
212 
213 ! *********************************************************************
214 
215  ff=one
216  do ii=2,nn
217    ff=ff*ii
218  end do
219 
220  factorial=ff
221 
222 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

SOURCE

1288 function fermi_dirac(energy, mu, temperature)
1289 
1290 !Arguments ------------------------------------
1291 !scalars
1292  real(dp),intent(in) :: energy, mu, temperature
1293  real(dp) :: fermi_dirac
1294 
1295 !Local variables-------------------------------
1296 !scalars
1297  real(dp) :: arg
1298 
1299 ! *************************************************************************
1300 
1301  fermi_dirac = zero
1302  if (temperature > tol12) then
1303    arg = (energy-mu)/temperature
1304    if(arg < -600._dp)then ! far below Ef
1305      fermi_dirac = one
1306    else if (arg < 600._dp)then ! around Ef
1307      fermi_dirac = one / (exp(arg)  + one)
1308    end if
1309  else  ! T is too small - just step function
1310    if (mu-energy > tol12) fermi_dirac = one
1311  end if
1312 
1313 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

SOURCE

524 elemental function gaussian(arg, sigma)
525 
526 !Arguments ---------------------------------------------
527 !scalars
528  real(dp),intent(in) :: arg,sigma
529  real(dp) :: gaussian
530 
531 !Local variables ---------------------------------------
532 !scalars
533  real(dp) :: xx
534 
535 ! *********************************************************************
536 
537  xx = arg / (sqrt2 * sigma)
538  gaussian = exp(-xx*xx) / (sigma * sqrt(two_pi))
539 
540 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.

SOURCE

2104 pure subroutine gspline_eval(self, x0, nx, xmesh, weights)
2105 
2106 !Arguments ------------------------------------
2107 !scalars
2108  integer,intent(in) :: nx
2109  real(dp),intent(in) :: x0
2110  type(gspline_t),intent(in) :: self
2111 !arrays
2112  real(dp),intent(in) :: xmesh(nx)
2113  real(dp),intent(out) :: weights(nx,2)
2114 
2115 !Local variables ------------------------------
2116 !scalars
2117  integer :: ix,jspl
2118  real(dp) :: xx,absx,aa,bb,cc,dd
2119  logical :: isneg
2120  !real(dp) :: int_values(nx)
2121 
2122 ! *************************************************************************
2123 
2124  do ix=1,nx
2125    xx = xmesh(ix) - x0; absx = abs(xx); isneg = xx < zero
2126    if (absx >= self%xmax) then
2127      ! Region in which gauss(x) is negligible.
2128      weights(ix,1) = zero
2129      if (isneg) then
2130        weights(ix,2) = zero
2131      else
2132        weights(ix,2) = one
2133      end if
2134    else
2135      ! Spline functions at |x| and recover the value at x:
2136      ! g(x) = g(-x); G(-x) = 1 - G(x)
2137      jspl = 1 + int((absx - self%xmin) * self%stepm1); dd = absx - self%xvals(jspl)
2138      bb = dd * self%stepm1
2139      aa = one - bb
2140      cc = aa*(aa**2-one) * self%step2div6
2141      dd = bb*(bb**2-one) * self%step2div6
2142 
2143      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)
2144      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)
2145      if (isneg) weights(ix,2) = one - weights(ix,2)
2146    end if
2147  end do
2148 
2149  !call simpson_int(nx,xmesh(2) - xmesh(1),weights(:,1),int_values)
2150  !do ix=1,nx
2151  !  write(99,*)xmesh(ix), weights(ix,1), gaussian(xx, self%sigma), weights(ix,2), int_values(ix)
2152  !end do
2153 
2154 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

SOURCE

2169 subroutine gspline_free(self)
2170 
2171 !Arguments ------------------------------------
2172 !scalars
2173  type(gspline_t),intent(inout) :: self
2174 
2175 ! *************************************************************************
2176 
2177  if (allocated(self%xvals)) then
2178    ABI_FREE(self%xvals)
2179  end if
2180  if (allocated(self%svals)) then
2181    ABI_FREE(self%svals)
2182  end if
2183 
2184 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.

SOURCE

2045 type (gspline_t) function gspline_new(sigma) result(new)
2046 
2047 !Arguments ------------------------------------
2048 !scalars
2049  real(dp),intent(in) :: sigma
2050 
2051 !Local variables ------------------------------
2052 !scalars
2053  integer :: ii
2054  real(dp) :: ybcbeg, ybcend
2055 ! *************************************************************************
2056 
2057  new%nspline = 5 * 1024; new%sigma = sigma
2058  ABI_CHECK(sigma > zero, sjoin("invalid sigma:", ftoa(sigma)))
2059  new%xmin = zero
2060  new%xmax = sigma * sqrt(-log(sigma * sqrt(pi) * tol12)) ! gauss(xmax) = tol12
2061  new%step = (new%xmax - new%xmin) / (new%nspline - 1)
2062  new%stepm1 = one / new%step; new%step2div6 = new%step**2 / six
2063 
2064  ABI_MALLOC(new%xvals, (new%nspline))
2065  do ii=1,new%nspline
2066    new%xvals(ii) = new%xmin + (ii-1) * new%step
2067  end do
2068  new%xmax = new%xvals(new%nspline)
2069 
2070  ! Spline the gaussian approximant.
2071  ABI_MALLOC(new%svals, (new%nspline, 4))
2072  new%svals(:, 1) = gaussian(new%xvals, sigma)
2073  ybcbeg = - (two * new%xmin / sigma**2) * new%svals(1,1)
2074  ybcend = - (two * new%xmax / sigma**2) * new%svals(new%nspline,1)
2075  call spline(new%xvals, new%svals(:,1), new%nspline, ybcbeg, ybcend, new%svals(:,2))
2076 
2077  ! Spline the primitive: 1/2 [1 + erf(x/sigma)]
2078  new%svals(:, 3) = half * (one + abi_derf(new%xvals / new%sigma))
2079  call spline(new%xvals, new%svals(:,3), new%nspline, new%svals(1,1), new%svals(new%nspline, 1), new%svals(:,4))
2080  !do ii=1,new%nspline; write(98,*)new%xvals(ii),new%svals(ii,3),new%svals(ii,4); end do
2081 
2082 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

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

SOURCE

450 recursive function IRadFnH(a,b,n,l,Z,m) result(x)
451 
452 !Arguments ---------------------------------------------
453 !scalars
454  integer,intent(in),optional  :: n,l,m
455  real(dp),intent(in):: a
456  real(dp),intent(in),optional :: b,Z
457 
458 !Local variables ---------------------------------------
459 !scalars
460  integer   :: nn,ll,mm
461  real(dp)  :: h,bb,ZZ,x
462 
463 ! *********************************************************************
464 
465  if (present(n)) then
466    nn=n
467  else
468    nn=3
469  end if
470 
471  if (present(l)) then
472    ll=l
473  else
474    ll=2
475  end if
476 
477  if (present(Z)) then
478    ZZ=Z
479  else
480    ZZ=28
481  end if
482 
483  if (present(b)) then
484    bb=b
485  else
486    bb=100.0_dp
487  end if
488 
489  if (present(m)) then
490    mm=m
491  else
492    mm=0
493  end if
494 
495  h=(bb-a)/2.0_dp
496  if (mm<8) then
497   !h=2*h/exp(1.0_dp)
498   x=IRadFnH(a,a+h,nn,ll,ZZ,mm+1)+IRadFnH(a+h,bb,nn,ll,ZZ,mm+1)
499  else
500   x=RadFnH(a,nn,ll,ZZ)**2*a**2+4.0_dp*RadFnH(a+h,nn,ll,ZZ)**2*(a+h)**2
501   x=h/3.0_dp*(x+RadFnH(bb,nn,ll,ZZ)**2*bb**2)
502  end if
503 
504 end function IRadFnH

m_special_funcs/jlspline_free [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 jlspline_free

FUNCTION

  deallocate memory

SOURCE

1966 subroutine jlspline_free(jlspl)
1967 
1968 !Arguments ------------------------------------
1969  type(jlspline_t),intent(inout) :: jlspl
1970 
1971 ! *********************************************************************
1972 
1973  if (allocated(jlspl%xx)) then
1974    ABI_FREE(jlspl%xx)
1975  end if
1976  if (allocated(jlspl%bess_spl)) then
1977    ABI_FREE(jlspl%bess_spl)
1978  end if
1979  if (allocated(jlspl%bess_spl_der)) then
1980    ABI_FREE(jlspl%bess_spl_der)
1981  end if
1982 
1983 end subroutine jlspline_free

m_special_funcs/jlspline_integral [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 jlspline_integral

INPUTS

OUTPUT

FUNCTION

SOURCE

2000 real(dp) function jlspline_integral(jlspl, il, qq, powr, nr, rcut)  result(res)
2001 
2002 !Arguments ------------------------------------
2003  integer,intent(in) :: il,nr,powr
2004  real(dp),intent(in) :: qq, rcut
2005  type(jlspline_t),intent(in) :: jlspl
2006 
2007 !Local variables ---------------------------------------
2008  integer :: ierr
2009  real(dp) :: step
2010 !arrays
2011  real(dp):: xfit(nr),yfit(nr),rr(nr)
2012 ! *********************************************************************
2013 
2014  step = rcut / (nr - 1)
2015  rr = arth(zero, step, nr)
2016  xfit = qq * rr
2017  call splint(jlspl%nx, jlspl%xx, jlspl%bess_spl(:,il), jlspl%bess_spl_der(:,il), nr, xfit, yfit, ierr=ierr)
2018 
2019  if (ierr /= 0) then
2020    write(std_out,*)"qq, rcut, qq*rcut, maxarg", qq, rcut, qq*rcut, jlspl%maxarg
2021    write(std_out,*)"x[0], x[-1]",jlspl%xx(1),jlspl%xx(jlspl%nx)
2022    write(std_out,*)"minval xfit: ",minval(xfit)
2023    write(std_out,*)"maxval xfit: ",maxval(xfit)
2024    ABI_ERROR("splint returned ierr != 0")
2025  end if
2026 
2027  if (powr /= 1) yfit = yfit * (rr ** powr)
2028  res = simpson(step, yfit)
2029 
2030 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

SOURCE

1895 type(jlspline_t) function jlspline_new(nx, delta, mlang) result(new)
1896 
1897 !Arguments ------------------------------------
1898 !scalars
1899  integer,intent(in) :: nx,mlang
1900  real(dp),intent(in) :: delta
1901 
1902 !Local variables -------------------------
1903 !scalars
1904  integer :: ix,ll
1905  real(dp) :: yp1,ypn
1906 !arrays
1907  real(dp),allocatable :: cosbessx(:),sinbessx(:)
1908 
1909 ! *********************************************************************
1910 
1911  if (nx < 2) then
1912    ABI_ERROR('need more than one point for the interpolation routines')
1913  end if
1914 
1915  new%nx = nx; new%mlang = mlang; new%delta = delta; new%maxarg = (nx-1) * delta
1916  ABI_MALLOC(new%xx, (nx))
1917  ABI_MALLOC(new%bess_spl, (nx, mlang))
1918  ABI_MALLOC(new%bess_spl_der, (nx, mlang))
1919 
1920  !-----------------------------------------------------------------
1921  !Bessel function into array
1922  !-----------------------------------------------------------------
1923  ! integration grid is nfiner times finer than the interpolation grid
1924  ABI_MALLOC(sinbessx, (nx))
1925  ABI_MALLOC(cosbessx, (nx))
1926 
1927  ! could be done by chain rule for cos sin (is it worth it?) but
1928  ! precision problems as numerical errors are propagated.
1929  do ix=1,nx
1930    new%xx(ix) = (ix-1) * delta
1931    sinbessx(ix) = sin(new%xx(ix))
1932    cosbessx(ix) = cos(new%xx(ix))
1933  end do
1934 
1935  ! fill bess_spl array
1936  do ll=0,mlang-1
1937    call besjm(one,new%bess_spl(:,ll+1),cosbessx,ll,nx,sinbessx,new%xx)
1938 
1939    ! call spline to get 2nd derivative (reuse in splint later)
1940    yp1 = zero; ypn = zero
1941    call spline(new%xx, new%bess_spl(:,ll+1), nx, yp1, ypn, new%bess_spl_der(:,ll+1))
1942  end do
1943 
1944 !write(std_out,*) ' bess funct  0   1   2   3   4'
1945 !do ix=1,nx
1946 !write(std_out,*) xx(ix), (new%bess_spl(ix,ll),ll=1,mlang)
1947 !end do
1948 
1949  ABI_FREE(sinbessx)
1950  ABI_FREE(cosbessx)
1951 
1952 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

71  type,public :: jlspline_t
72 
73    integer :: nx
74    ! number of points on linear mesh used in spline.
75 
76    integer :: mlang
77    ! mlang= max angular momentum + 1
78 
79    real(dp) :: delta
80    ! Step of linear mesh.
81 
82    real(dp) :: maxarg
83    ! max arg value.
84 
85    real(dp),allocatable :: xx(:)
86    ! xx(nx)
87    ! coordinates of points belonging to the grid
88 
89    real(dp),allocatable :: bess_spl(:,:)
90    ! bess_spl(nx,mlang)
91    ! bessel functions computed on the linear mesh
92 
93    real(dp),allocatable :: bess_spl_der(:,:)
94    ! bess_spl_der(nx,mlang)
95    ! the second derivatives of the cubic spline.
96 
97  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.

SOURCE

1795 elemental function k_fermi(rhor)
1796 
1797 !Arguments ------------------------------------
1798 !scalars
1799  real(dp),intent(in) :: rhor
1800  real(dp) :: k_fermi
1801 
1802 !Local variables-------------------------------
1803 !scalars
1804  real(dp),parameter :: pisq=pi**2
1805 
1806 ! *************************************************************************
1807 
1808  k_fermi = (three*pisq*rhor)**third
1809 
1810 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.

SOURCE

1827 elemental function k_thfermi(rhor)
1828 
1829 !Arguments ------------------------------------
1830 !scalars
1831  real(dp),intent(in) :: rhor
1832  real(dp) :: k_thfermi
1833 
1834 !Local variables-------------------------------
1835 !scalars
1836  real(dp),parameter :: pisq=pi**2
1837 
1838 ! *************************************************************************
1839 
1840  k_thfermi = SQRT(four*k_fermi(rhor)*piinv)
1841 
1842 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)

SOURCE

323 function laguerre(x,n,a)
324 
325 !Arguments ---------------------------------------------
326 !scalars
327  integer,intent(in),optional :: n,a
328  real(dp)                    :: laguerre
329  real(dp),intent(in)         :: x
330 
331 !Local variables ---------------------------------------
332 !scalars
333  integer :: ii, nn, aa
334 
335 !arrays
336  real(dp),allocatable :: ff(:)
337 
338 ! *********************************************************************
339 
340  if (present(n)) then
341    nn=n
342  else
343    nn=1
344  end if
345 
346  if (present(a)) then
347    aa=a
348  else
349    aa=0
350  end if
351  ABI_MALLOC(ff,(nn+1))
352  ff=0.0_dp
353  ff=(/ (binomcoeff(nn+aa,nn-ii)*((-1.0_dp)*x)**ii/factorial(ii) ,ii=0,nn) /)
354  laguerre=sum(ff)
355 
356  ABI_FREE(ff)
357 
358 end function laguerre

m_special_funcs/levi_civita_3 [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

  levi_civita_3

FUNCTION

 Return Levi-Civita tensor of rank 3

SOURCE

1856 pure function levi_civita_3() result(ee)
1857 
1858 !Arguments ------------------------------------
1859  integer :: ee(3,3,3)
1860 
1861 ! *************************************************************************
1862 
1863  ee = 0
1864  ee(1,2,3) = 1
1865  ee(2,3,1) = 1
1866  ee(3,1,2) = 1
1867  !
1868  ee(3,2,1) = -1
1869  ee(1,3,2) = -1
1870  ee(2,1,3) = -1
1871 
1872 end function levi_civita_3

m_special_funcs/lorentzian [ Functions ]

[ Top ] [ m_special_funcs ] [ Functions ]

NAME

 lorentzian

FUNCTION

  Lorentzian function.

INPUTS

   arg=Argument of the lorentzian.
   sigma=Broadening factor

SOURCE

558 elemental function lorentzian(arg, sigma)
559 
560 !Arguments ---------------------------------------------
561 !scalars
562  real(dp),intent(in) :: arg, sigma
563  real(dp) :: lorentzian
564 
565 ! *********************************************************************
566 
567  lorentzian = piinv * sigma / (arg ** 2 + sigma ** 2)
568 
569 end function lorentzian

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)

SOURCE

243 pure function permutations(nn,kk)
244 
245 !Arguments ---------------------------------------------
246 !scalars
247  integer,intent(in) :: kk,nn
248  real(dp) :: permutations
249 
250 !Local variables ---------------------------------------
251 !scalars
252  integer :: ii
253  real(dp) :: pp
254 
255 ! *********************************************************************
256 
257  if ((nn>=0).and.((nn-kk)>=0)) then
258    pp=one
259    do ii=nn-kk+1,nn
260      pp=pp*ii
261    end do
262  else
263    pp=zero
264  end if
265 
266  permutations=pp
267 
268 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)

SOURCE

384 function RadFnH(r,n,l,Z)
385 
386 !Arguments ---------------------------------------------
387 !scalars
388  integer,intent(in),optional :: n,l
389  real(dp) :: RadFnH
390  real(dp),intent(in) :: r
391  real(dp),intent(in),optional :: Z
392 
393 !Local variables ---------------------------------------
394 !scalars
395  integer   :: nn,ll
396  real(dp)  :: ff,rr,ZZ
397 
398 ! *********************************************************************
399 
400  if (present(n)) then
401    nn=n
402  else
403    nn=3
404  end if
405 
406  if (present(l)) then
407    ll=l
408  else
409    ll=2
410  end if
411 
412  if (present(Z)) then
413    ZZ=Z
414  else
415    ZZ=28.0_dp
416  end if
417 
418  rr=ZZ*r/nn
419  ff=exp(log(ZZ*1.0_dp)*(3.0_dp/2.0_dp))*2/nn**2
420  ff=ff*sqrt(factorial(nn-ll-1)/factorial(nn+ll))*(2*rr)**ll
421  RadFnH=ff*exp(-1*rr)*laguerre(2*rr,nn-ll-1,2*ll+1)
422 
423 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

SOURCE

1214 subroutine sbf8(nm,xx,sb_out)
1215 
1216 !Arguments----------------------------------------------------------
1217 !scalars
1218  integer,intent(in) :: nm
1219  real(dp),intent(in) :: xx
1220 !arrays
1221  real(dp),intent(out) :: sb_out(nm)
1222 
1223 !Local variables-------------------------------
1224 !scalars
1225  integer :: nlim,nn
1226  real(dp) :: fn,sn,xi,xn,xs
1227 !arrays
1228  real(dp),allocatable :: sb(:)
1229 
1230 ! *************************************************************************
1231 
1232  if(xx<= 1.0e-36_dp) then
1233 !  zero argument section
1234    sb_out(:)=zero
1235    sb_out(1)=one
1236  else if(xx<1.e-3_dp) then
1237 !  small argument section
1238    xn=one
1239    xs=half*xx**2
1240    do nn=1,nm
1241      sb_out(nn)=xn*(one - xs*(one - xs/(4*nn+6))/(2*nn+1))
1242      xn=xx*xn/(2*nn+1)
1243    end do
1244  else
1245 !  recursion method
1246    if(xx<one) then
1247      nlim=nm+int(15.0e0_dp*xx)+1
1248    else
1249      nlim=nm+int(1.36e0_dp*xx)+15
1250    end if
1251    ABI_MALLOC(sb,(nlim+1))
1252    nn=nlim
1253    xi=one/xx
1254    sb(nn+1)=zero
1255    sb(nn)=1.e-18_dp
1256    sn=dble(2*nn-1)*1.e-36_dp
1257    do nn=nlim-1,1,-1
1258      sb(nn)=dble(2*nn+1)*xi*sb(nn+1) - sb(nn+2)
1259    end do
1260    do nn=1,nlim-1
1261      sn=sn + dble(2*nn-1)*sb(nn)*sb(nn)
1262    end do
1263    fn=1.d0/sqrt(sn)
1264    sb_out(:)=fn*sb(1:nm)
1265    ABI_FREE(sb)
1266  end if
1267 
1268 end subroutine sbf8