TABLE OF CONTENTS
- ABINIT/GAMMA_FUNCTION
- ABINIT/m_special_funcs
- m_special_funcs/abi_derf
- m_special_funcs/abi_derfc
- m_special_funcs/besjm
- m_special_funcs/binomcoeff
- m_special_funcs/bose_einstein
- m_special_funcs/clp
- m_special_funcs/dip12
- m_special_funcs/dip32
- m_special_funcs/djp12
- m_special_funcs/djp32
- m_special_funcs/factorial
- m_special_funcs/fermi_dirac
- m_special_funcs/gaussian
- m_special_funcs/gspline_eval
- m_special_funcs/gspline_free
- m_special_funcs/gspline_new
- m_special_funcs/gspline_t
- m_special_funcs/IRadFnH
- m_special_funcs/jlspline_free
- m_special_funcs/jlspline_integral
- m_special_funcs/jlspline_new
- m_special_funcs/jlspline_t
- m_special_funcs/k_fermi
- m_special_funcs/k_thfermi
- m_special_funcs/laguerre
- m_special_funcs/levi_civita_3
- m_special_funcs/lorentzian
- m_special_funcs/permutations
- m_special_funcs/RadFnH
- m_special_funcs/sbf8
ABINIT/GAMMA_FUNCTION [ 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 ]
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