TABLE OF CONTENTS
- ABINIT/m_paw_numeric
- m_paw_numeric/paw_derfc
- m_paw_numeric/paw_jbessel
- m_paw_numeric/paw_smooth
- m_paw_numeric/paw_solvbes
- m_paw_numeric/paw_sort_dp
- m_paw_numeric/paw_spline
- m_paw_numeric/paw_splint
- m_paw_numeric/paw_splint_der
- m_paw_numeric/paw_uniform_splfit
- m_special_funcs/paw_jbessel_4spline
ABINIT/m_paw_numeric [ Modules ]
NAME
m_paw_numeric
FUNCTION
Wrappers for various numeric operations (spline, sort, ...)
COPYRIGHT
Copyright (C) 2012-2024 ABINIT group (MT,TR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
FOR DEVELOPPERS: in order to preserve the portability of libPAW library, please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
SOURCE
20 #include "libpaw.h" 21 22 module m_paw_numeric 23 24 USE_DEFS 25 USE_MSG_HANDLING 26 USE_MEMORY_PROFILING 27 28 implicit none 29 30 private 31 32 !public procedures 33 public:: paw_spline 34 public:: paw_splint 35 public:: paw_splint_der 36 public:: paw_uniform_splfit 37 public:: paw_smooth 38 public:: paw_sort_dp 39 public:: paw_jbessel 40 public:: paw_solvbes 41 public:: paw_jbessel_4spline 42 public:: paw_derfc
m_paw_numeric/paw_derfc [ Functions ]
[ Top ] [ m_paw_numeric ] [ Functions ]
NAME
paw_derfc
FUNCTION
Evaluates the complementary error function in real(dp).
INPUTS
yy
OUTPUT
derfc_yy=complementary error function of yy
SOURCE
997 elemental function paw_derfc(yy) result(derfc_yy) 998 999 !Arguments ------------------------------------ 1000 !scalars 1001 real(dp),intent(in) :: yy 1002 real(dp) :: derfc_yy 1003 1004 !Local variables------------------------------- 1005 integer :: done,ii,isw 1006 ! coefficients for 0.0 <= yy < .477 1007 real(dp), parameter :: & 1008 & pp(5)=(/ 113.8641541510502e0_dp, 377.4852376853020e0_dp, & 1009 & 3209.377589138469e0_dp, .1857777061846032e0_dp, & 1010 & 3.161123743870566e0_dp /) 1011 real(dp), parameter :: & 1012 & qq(4)=(/ 244.0246379344442e0_dp, 1282.616526077372e0_dp, & 1013 & 2844.236833439171e0_dp, 23.60129095234412e0_dp/) 1014 ! coefficients for .477 <= yy <= 4.0 1015 real(dp), parameter :: & 1016 & p1(9)=(/ 8.883149794388376e0_dp, 66.11919063714163e0_dp, & 1017 & 298.6351381974001e0_dp, 881.9522212417691e0_dp, & 1018 & 1712.047612634071e0_dp, 2051.078377826071e0_dp, & 1019 & 1230.339354797997e0_dp, 2.153115354744038e-8_dp, & 1020 & .5641884969886701e0_dp /) 1021 real(dp), parameter :: & 1022 & q1(8)=(/ 117.6939508913125e0_dp, 537.1811018620099e0_dp, & 1023 & 1621.389574566690e0_dp, 3290.799235733460e0_dp, & 1024 & 4362.619090143247e0_dp, 3439.367674143722e0_dp, & 1025 & 1230.339354803749e0_dp, 15.74492611070983e0_dp/) 1026 ! coefficients for 4.0 < y, 1027 real(dp), parameter :: & 1028 & p2(6)=(/ -3.603448999498044e-01_dp, -1.257817261112292e-01_dp, & 1029 & -1.608378514874228e-02_dp, -6.587491615298378e-04_dp, & 1030 & -1.631538713730210e-02_dp, -3.053266349612323e-01_dp/) 1031 real(dp), parameter :: & 1032 & q2(5)=(/ 1.872952849923460e0_dp , 5.279051029514284e-01_dp, & 1033 & 6.051834131244132e-02_dp , 2.335204976268692e-03_dp, & 1034 & 2.568520192289822e0_dp /) 1035 real(dp), parameter :: & 1036 & sqrpi=.5641895835477563e0_dp, xbig=13.3e0_dp, xlarge=6.375e0_dp, xmin=1.0e-10_dp 1037 real(dp) :: res,xden,xi,xnum,xsq,xx 1038 1039 !****************************************************************** 1040 1041 xx = yy 1042 isw = 1 1043 !Here change the sign of xx, and keep track of it thanks to isw 1044 if (xx<0.0e0_dp) then 1045 isw = -1 1046 xx = -xx 1047 end if 1048 1049 done=0 1050 1051 !Residual value, if yy < -6.375e0_dp 1052 res=2.0e0_dp 1053 1054 !abs(yy) < .477, evaluate approximation for erfc 1055 if (xx<0.477e0_dp) then 1056 ! xmin is a very small number 1057 if (xx<xmin) then 1058 res = xx*pp(3)/qq(3) 1059 else 1060 xsq = xx*xx 1061 xnum = pp(4)*xsq+pp(5) 1062 xden = xsq+qq(4) 1063 do ii = 1,3 1064 xnum = xnum*xsq+pp(ii) 1065 xden = xden*xsq+qq(ii) 1066 end do 1067 res = xx*xnum/xden 1068 end if 1069 if (isw==-1) res = -res 1070 res = 1.0e0_dp-res 1071 done=1 1072 end if 1073 1074 !.477 < abs(yy) < 4.0 , evaluate approximation for erfc 1075 if (xx<=4.0e0_dp .and. done==0 ) then 1076 xsq = xx*xx 1077 xnum = p1(8)*xx+p1(9) 1078 xden = xx+q1(8) 1079 do ii=1,7 1080 xnum = xnum*xx+p1(ii) 1081 xden = xden*xx+q1(ii) 1082 end do 1083 res = xnum/xden 1084 res = res* exp(-xsq) 1085 if (isw.eq.-1) res = 2.0e0_dp-res 1086 done=1 1087 end if 1088 1089 !y > 13.3e0_dp 1090 if (isw > 0 .and. xx > xbig .and. done==0 ) then 1091 res = 0.0e0_dp 1092 done=1 1093 end if 1094 1095 !4.0 < yy < 13.3e0_dp .or. -6.375e0_dp < yy < -4.0 1096 !evaluate minimax approximation for erfc 1097 if ( ( isw > 0 .or. xx < xlarge ) .and. done==0 ) then 1098 xsq = xx*xx 1099 xi = 1.0e0_dp/xsq 1100 xnum= p2(5)*xi+p2(6) 1101 xden = xi+q2(5) 1102 do ii = 1,4 1103 xnum = xnum*xi+p2(ii) 1104 xden = xden*xi+q2(ii) 1105 end do 1106 res = (sqrpi+xi*xnum/xden)/xx 1107 res = res* exp(-xsq) 1108 if (isw.eq.-1) res = 2.0e0_dp-res 1109 end if 1110 1111 !All cases have been investigated 1112 derfc_yy = res 1113 1114 end function paw_derfc
m_paw_numeric/paw_jbessel [ Functions ]
[ Top ] [ m_paw_numeric ] [ Functions ]
NAME
paw_jbessel
FUNCTION
Compute spherical Bessel function j_l(x) and derivative(s)
INPUTS
ll=l-order of the Bessel function order=1 if first derivative is requested 2 if first and second derivatives are requested xx=where to compute j_l
OUTPUT
bes= Bessel function j_l at xx besp= first derivative of j_l at xx (only if order>=1) bespp= second derivative of j_l at xx (only if order=2)
SOURCE
686 subroutine paw_jbessel(bes,besp,bespp,ll,order,xx) 687 688 !Arguments --------------------------------------------- 689 !scalars 690 integer,intent(in) :: ll,order 691 real(dp),intent(in) :: xx 692 real(dp),intent(out) :: bes,besp,bespp 693 694 !Local variables --------------------------------------- 695 !scalars 696 integer,parameter :: imax=40 697 integer :: ii,il 698 real(dp),parameter :: prec=1.d-15 699 real(dp) :: besp1,fact,factp,factpp,jn,jnp,jnpp,jr,xx2,xxinv 700 character(len=200) :: msg 701 702 ! ********************************************************************* 703 704 if (order>2) then 705 msg='Wrong order in paw_jbessel!' 706 LIBPAW_ERROR(msg) 707 end if 708 709 if (abs(xx)<prec) then 710 bes=zero;if (ll==0) bes=one 711 if (order>=1) then 712 besp=zero;if (ll==1) besp=third 713 end if 714 if (order==2) then 715 bespp=zero 716 if (ll==0) bespp=-third 717 if (ll==2) bespp=2._dp/15._dp 718 end if 719 return 720 end if 721 722 xxinv=one/xx 723 if (order==0) then 724 factp=zero 725 factpp=zero 726 jnp=zero 727 jnpp=zero 728 end if 729 730 if (xx<one) then 731 xx2=0.5_dp*xx*xx 732 fact=one 733 do il=1,ll 734 fact=fact*xx/dble(2*il+1) 735 end do 736 jn=one;jr=one;ii=0 737 do while(abs(jr)>=prec.and.ii<imax) 738 ii=ii+1;jr=-jr*xx2/dble(ii*(2*(ll+ii)+1)) 739 jn=jn+jr 740 end do 741 bes=jn*fact 742 if (abs(jr)>prec) then 743 msg='Bessel function did not converge!' 744 LIBPAW_ERROR(msg) 745 end if 746 if (order>=1) then 747 factp=fact*xx/dble(2*ll+3) 748 jnp=one;jr=one;ii=0 749 do while(abs(jr)>=prec.AND.ii<imax) 750 ii=ii+1;jr=-jr*xx2/dble(ii*(2*(ll+ii)+3)) 751 jnp=jnp+jr 752 end do 753 besp=-jnp*factp+jn*fact*xxinv*dble(ll) 754 if (abs(jr)>prec) then 755 msg='1st der. of Bessel function did not converge!' 756 LIBPAW_ERROR(msg) 757 end if 758 end if 759 if (order==2) then 760 factpp=factp*xx/dble(2*ll+5) 761 jnpp=one;jr=one;ii=0 762 do while(abs(jr)>=prec.AND.ii<imax) 763 ii=ii+1;jr=-jr*xx2/dble(ii*(2*(ll+ii)+5)) 764 jnpp=jnpp+jr 765 end do 766 besp1=-jnpp*factpp+jnp*factp*xxinv*dble(ll+1) 767 if (abs(jr)>prec) then 768 msg='2nd der. of Bessel function did not converge !' 769 LIBPAW_ERROR(msg) 770 end if 771 end if 772 else 773 jn =sin(xx)*xxinv 774 jnp=(-cos(xx)+jn)*xxinv 775 do il=2,ll+1 776 jr=-jn+dble(2*il-1)*jnp*xxinv 777 jn=jnp;jnp=jr 778 end do 779 bes=jn 780 if (order>=1) besp =-jnp+jn *xxinv*dble(ll) 781 if (order==2) besp1= jn -jnp*xxinv*dble(ll+2) 782 end if 783 784 if (order==2) bespp=-besp1+besp*ll*xxinv-bes*ll*xxinv*xxinv 785 786 end subroutine paw_jbessel
m_paw_numeric/paw_smooth [ Functions ]
[ Top ] [ m_paw_numeric ] [ Functions ]
NAME
paw_smooth
FUNCTION
Smooth an array of given ordinates (y's) that are in order of increasing abscissas (x's), but without using the abscissas themselves supposed to be equally spaced.
INPUTS
it=number of abscissas to treat mesh=size of the array (number of abscissas)
OUTPUT
SIDE EFFECTS
a(mesh)=array to be smoothed
SOURCE
526 subroutine paw_smooth(a,mesh,it) 527 528 !Arguments ------------------------------------ 529 !scalars 530 integer, intent(in) :: it,mesh 531 !arrays 532 real(dp), intent(inout) :: a(mesh) 533 534 !Local variables------------------------------- 535 !scalars 536 integer :: i,k 537 !arrays 538 real(dp) :: asm(mesh) 539 540 ! ************************************************************************* 541 542 asm(1:4) = zero ! ?? Correct me ... 543 do k=1,it 544 asm(5)=0.2_dp*(a(3)+a(4)+a(5)+a(6)+a(7)) 545 asm(mesh-4)=0.2_dp*(a(mesh-2)+a(mesh-3)+a(mesh-4)+& 546 & a(mesh-5)+a(mesh-6)) 547 asm(mesh-3)=0.2_dp*(a(mesh-1)+a(mesh-2)+a(mesh-3)+& 548 & a(mesh-4)+a(mesh-5)) 549 asm(mesh-2)=0.2_dp*(a(mesh)+a(mesh-1)+a(mesh-2)+& 550 & a(mesh-3)+a(mesh-4)) 551 asm(mesh-1)=0.25_dp*(a(mesh)+a(mesh-1)+a(mesh-2)+a(mesh-3)) 552 asm(mesh)=1.0_dp/3.0_dp*(a(mesh)+a(mesh-1)+a(mesh-2)) 553 do i=6,mesh-5 554 asm(i)=0.1_dp *a(i)+0.1_dp*(a(i+1)+a(i-1))+& 555 & 0.1_dp *(a(i+2)+a(i-2))+& 556 & 0.1_dp *(a(i+3)+a(i-3))+& 557 & 0.1_dp *(a(i+4)+a(i-4))+& 558 & 0.05_dp*(a(i+5)+a(i-5)) 559 end do 560 do i=1,mesh 561 a(i)=asm(i) 562 end do 563 end do 564 565 end subroutine paw_smooth
m_paw_numeric/paw_solvbes [ Functions ]
[ Top ] [ m_paw_numeric ] [ Functions ]
NAME
paw_solvbes
FUNCTION
Find nq first roots of instrinsic equation: alpha.jl(Q) + beta.Q.djl/dr(Q) = 0
INPUTS
alpha,beta= factors in intrinsic equation ll= l quantum number nq= number of roots to find
OUTPUT
root(nq)= roots of instrinsic equation
SOURCE
809 subroutine paw_solvbes(root,alpha,beta,ll,nq) 810 811 !Arguments ------------------------------------ 812 !scalars 813 integer :: ll,nq 814 real(dp) :: alpha,beta 815 !arrays 816 real(dp) :: root(nq) 817 818 !Local variables------------------------------- 819 !scalars 820 integer :: nroot 821 real(dp),parameter :: dh=0.1_dp,tol=tol14 822 real(dp) :: dum,hh,jbes,jbesp,qq,qx,y1,y2 823 824 ! ************************************************************************* 825 826 qq=dh;nroot=0 827 828 do while (nroot<nq) 829 call paw_jbessel(jbes,jbesp,dum,ll,1,qq) 830 y1=alpha*jbes+beta*qq*jbesp 831 qq=qq+dh 832 call paw_jbessel(jbes,jbesp,dum,ll,1,qq) 833 y2=alpha*jbes+beta*qq*jbesp 834 835 do while (y1*y2>=zero) 836 qq=qq+dh 837 call paw_jbessel(jbes,jbesp,dum,ll,1,qq) 838 y2=alpha*jbes+beta*qq*jbesp 839 end do 840 841 hh=dh;qx=qq 842 do while (hh>tol) 843 hh=half*hh 844 if (y1*y2<zero) then 845 qx=qx-hh 846 else 847 qx=qx+hh 848 end if 849 call paw_jbessel(jbes,jbesp,dum,ll,1,qx) 850 y2=alpha*jbes+beta*qx*jbesp 851 end do 852 nroot=nroot+1 853 root(nroot)=qx 854 855 end do 856 857 end subroutine paw_solvbes
m_paw_numeric/paw_sort_dp [ Functions ]
[ Top ] [ m_paw_numeric ] [ Functions ]
NAME
paw_sort_dp
FUNCTION
Sort real(dp) array list(n) into ascending numerical order using Heapsort algorithm, while making corresponding rearrangement of the integer array iperm. Consider that two real(dp) numbers within tolerance tol are equal.
INPUTS
n intent(in) dimension of the list tol intent(in) numbers within tolerance are equal list(n) intent(inout) list of real(dp) numbers to be sorted iperm(n) intent(inout) iperm(i)=i (very important)
OUTPUT
list(n) sorted list iperm(n) index of permutation given the right ascending order
SOURCE
592 subroutine paw_sort_dp(n,list,iperm,tol) 593 594 !Arguments ------------------------------------ 595 !scalars 596 integer, intent(in) :: n 597 real(dp), intent(in) :: tol 598 !arrays 599 integer, intent(inout) :: iperm(n) 600 real(dp), intent(inout) :: list(n) 601 602 !Local variables------------------------------- 603 !scalars 604 integer :: l,ir,iap,i,j 605 real(dp) :: ap 606 character(len=500) :: msg 607 !arrays 608 609 ! ************************************************************************* 610 611 !Accomodate case of array of length 1: already sorted! 612 if (n==1) return 613 614 !Should not call with n<1 615 if (n<1) then 616 write(msg,'(a,i12,2a)') & 617 & 'paw_sort_dp has been called with array length n=',n, ch10, & 618 & ' having a value less than 1. This is not allowed.' 619 LIBPAW_ERROR(msg) 620 end if 621 622 !Conduct the usual sort 623 l=n/2+1 ; ir=n 624 do ! Infinite do-loop 625 if (l>1) then 626 l=l-1 627 ap=list(l) 628 iap=iperm(l) 629 else ! l<=1 630 ap=list(ir) 631 iap=iperm(ir) 632 list(ir)=list(1) 633 iperm(ir)=iperm(1) 634 ir=ir-1 635 if (ir==1) then 636 list(1)=ap 637 iperm(1)=iap 638 exit ! This is the end of this algorithm 639 end if 640 end if ! l>1 641 i=l 642 j=l+l 643 do while (j<=ir) 644 if (j<ir) then 645 if ( list(j)<list(j+1)-tol .or. & 646 & (list(j)<list(j+1)+tol.and.iperm(j)<iperm(j+1))) j=j+1 647 endif 648 if (ap<list(j)-tol.or.(ap<list(j)+tol.and.iap<iperm(j))) then 649 list(i)=list(j) 650 iperm(i)=iperm(j) 651 i=j 652 j=j+j 653 else 654 j=ir+1 655 end if 656 end do 657 list(i)=ap 658 iperm(i)=iap 659 end do ! End infinite do-loop 660 661 end subroutine paw_sort_dp
m_paw_numeric/paw_spline [ Functions ]
[ Top ] [ m_paw_numeric ] [ Functions ]
NAME
paw_spline
FUNCTION
Computes the second derivatives of a cubic spline
INPUTS
* Input, integer N, the number of data points; N must be at least 2. In the special case where N = 2 and IBCBEG = IBCEND = 0, the spline will actually be linear. * Input, real(dp) T(N), the knot values, that is, the points where data is specified. The knot values should be distinct, and increasing. * Input, real(dp) Y(N), the data values to be interpolated. * Input, real(dp) YBCBEG, YBCEND, the values to be used in the boundary conditions if IBCBEG or IBCEND is equal to 1 or 2.
OUTPUT
Output, real(dp) YPP(N), the second derivatives of the cubic spline. Work space, real(dp) DIAG(N) - should be removed ...
SOURCE
74 subroutine paw_spline(t,y,n,ybcbeg,ybcend,ypp) 75 76 !******************************************************************************* 77 ! Discussion: 78 ! For data interpolation, the user must call SPLINE_CUBIC_SET to 79 ! determine the second derivative data, passing in the data to be 80 ! interpolated, and the desired boundary conditions. 81 ! The data to be interpolated, plus the SPLINE_CUBIC_SET output, 82 ! defines the spline. The user may then call SPLINE_CUBIC_VAL to 83 ! evaluate the spline at any point. 84 ! The cubic spline is a piecewise cubic polynomial. The intervals 85 ! are determined by the "knots" or abscissas of the data to be 86 ! interpolated. The cubic spline has continous first and second 87 ! derivatives over the entire interval of interpolation. 88 ! For any point T in the interval T(IVAL), T(IVAL+1), the form of 89 ! the spline is 90 ! SPL(T) = A(IVAL) 91 ! + B(IVAL) * ( T - T(IVAL) ) 92 ! + C(IVAL) * ( T - T(IVAL) )**2 93 ! + D(IVAL) * ( T - T(IVAL) )**3 94 ! If we assume that we know the values Y(*) and YPP(*), which represent 95 ! the values and second derivatives of the spline at each knot, then 96 ! the coefficients can be computed as: 97 ! A(IVAL) = Y(IVAL) 98 ! B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) ) 99 ! - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6 100 ! C(IVAL) = YPP(IVAL) / 2 101 ! D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) ) 102 ! Since the first derivative of the spline is 103 ! SPL'(T) = B(IVAL) 104 ! + 2 * C(IVAL) * ( T - T(IVAL) ) 105 ! + 3 * D(IVAL) * ( T - T(IVAL) )**2, 106 ! the requirement that the first derivative be continuous at interior 107 ! knot I results in a total of N-2 equations, of the form: 108 ! B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1)) 109 ! + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))**2 = B(IVAL) 110 ! or, setting H(IVAL) = T(IVAL+1) - T(IVAL) 111 ! ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1) 112 ! - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6 113 ! + YPP(IVAL-1) * H(IVAL-1) 114 ! + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2 115 ! = 116 ! ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) 117 ! - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6 118 ! or 119 ! YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) ) 120 ! + YPP(IVAL) * H(IVAL) 121 ! = 122 ! 6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL) 123 ! - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1) 124 ! Boundary conditions must be applied at the first and last knots. 125 ! The resulting tridiagonal system can be solved for the YPP values. 126 ! 127 ! Author: 128 ! John Burkardt, modified by Xavier Gonze 129 130 !Arguments ------------------------------------ 131 !scalars 132 integer,intent(in) :: n 133 real(dp),intent(in) :: ybcbeg,ybcend 134 !arrays 135 real(dp),intent(in) :: t(n),y(n) 136 real(dp),intent(out) :: ypp(n) 137 138 !Local variables------------------------------- 139 !scalars 140 integer :: ibcbeg,ibcend,i,k 141 real(dp) :: ratio,pinv 142 character(len=500) :: msg 143 !arrays 144 real(dp),allocatable :: tmp(:) 145 146 ! ************************************************************************* 147 148 !Check 149 if (n<=1) then 150 write(msg,'(6a,i8)') ch10, & 151 & 'SPLINE_CUBIC_SET - Fatal error!',ch10, & 152 & ' The number of knots must be at least 2.',ch10, & 153 & ' The input value of N = ', n 154 LIBPAW_ERROR(msg) 155 end if 156 157 LIBPAW_ALLOCATE(tmp,(n)) 158 159 do i=1,n-1 160 if (t(i)>=t(i+1)) then 161 write(msg,'(6a,i8,a,es19.12,2a,i8,a,es19.12)') ch10, & 162 & 'SPLINE_CUBIC_SET - Fatal error!',ch10, & 163 & ' The knots must be strictly increasing, but',ch10, & 164 & ' T(', i,') = ', t(i), ch10, & 165 & ' T(',i+1,') = ', t(i+1) 166 LIBPAW_ERROR(msg) 167 end if 168 end do 169 170 ibcbeg=1;if(ybcbeg>1.0d+30)ibcbeg=0 171 ibcend=1;if(ybcend>1.0d+30)ibcend=0 172 173 !Set the first and last equations 174 if (ibcbeg==0) then 175 ypp(1) = 0._dp 176 tmp(1) = 0._dp 177 else if ( ibcbeg == 1 ) then 178 ypp(1) = -0.5_dp 179 tmp(1) = (3._dp/(t(2)-t(1)))*((y(2)-y(1))/(t(2)-t(1))-ybcbeg) 180 end if 181 if (ibcend==0) then 182 ypp(n) = 0._dp 183 tmp(n) = 0._dp 184 else if ( ibcend == 1 ) then 185 ypp(n) = 0.5_dp 186 tmp(n) = (3._dp/(t(n)-t(n-1)))*(ybcend-(y(n)-y(n-1))/(t(n)-t(n-1))) 187 end if 188 189 !Set the intermediate equations 190 do i=2,n-1 191 ratio=(t(i)-t(i-1))/(t(i+1)-t(i-1)) 192 pinv = 1.0_dp/(ratio*ypp(i-1) + 2.0_dp) 193 ypp(i) = (ratio-1.0_dp)*pinv 194 tmp(i)=(6.0_dp*((y(i+1)-y(i))/(t(i+1)-t(i))-(y(i)-y(i-1)) & 195 & /(t(i)-t(i-1)))/(t(i+1)-t(i-1))-ratio*tmp(i-1))*pinv 196 if (abs(tmp(i))<1.d5*tiny(0._dp)) tmp(i)=0._dp 197 end do 198 199 !Solve the equations 200 ypp(n) = (tmp(n)-ypp(n)*tmp(n-1))/(ypp(n)*ypp(n-1)+1.0_dp) 201 do k=n-1,1,-1 202 ypp(k)=ypp(k)*ypp(k+1)+tmp(k) 203 end do 204 205 LIBPAW_DEALLOCATE(tmp) 206 207 end subroutine paw_spline
m_paw_numeric/paw_splint [ Functions ]
[ Top ] [ m_paw_numeric ] [ Functions ]
NAME
paw_splint
FUNCTION
Compute spline interpolation of a tabulated function. There is no hypothesis about the spacing of the input grid points.
INPUTS
nspline: number of grid points of input mesh xspline(nspline): input mesh yspline(nspline): function on input mesh ysplin2(nspline): second derivative of yspline on input mesh nfit: number of points of output mesh xfit(nfit): output mesh
OUTPUT
yfit(nfit): function on output mesh [ierr]=A non-zero value is used to signal that some points in xfit exceed xspline(nspline). The input value is incremented by the number of such points.
SOURCE
235 subroutine paw_splint(nspline,xspline,yspline,ysplin2,nfit,xfit,yfit,ierr) 236 237 !Arguments ------------------------------------ 238 !scalars 239 integer,intent(in) :: nfit, nspline 240 integer,optional,intent(out) :: ierr 241 !arrays 242 real(dp),intent(in) :: xspline(nspline),yspline(nspline) 243 real(dp),intent(in) :: ysplin2(nspline),xfit(nfit) 244 real(dp),intent(out) :: yfit(nfit) 245 246 !Local variables------------------------------- 247 !scalars 248 integer :: left,i,k,right,my_err 249 real(dp) :: delarg,invdelarg,aa,bb 250 character(len=50) :: msg 251 !arrays 252 253 ! ************************************************************************* 254 255 my_err=0 256 left=1 257 do i=1,nfit 258 yfit(i)=0._dp ! Initialize for the unlikely event that rmax exceed r(mesh) 259 do k=left+1, nspline 260 if(xspline(k) >= xfit(i)) then 261 if(xspline(k-1) <= xfit(i)) then 262 right = k 263 left = k-1 264 else 265 if (k-1.eq.1 .and. i.eq.1) then 266 msg='xfit(1) < xspline(1)' 267 else 268 msg='xfit not properly ordered' 269 end if 270 LIBPAW_ERROR(msg) 271 end if 272 delarg= xspline(right) - xspline(left) 273 invdelarg= 1.0_dp/delarg 274 aa= (xspline(right)-xfit(i))*invdelarg 275 bb= (xfit(i)-xspline(left))*invdelarg 276 yfit(i) = aa*yspline(left)+bb*yspline(right) & 277 & +( (aa*aa*aa-aa)*ysplin2(left) + & 278 & (bb*bb*bb-bb)*ysplin2(right) ) *delarg*delarg/6.0_dp 279 exit 280 end if 281 end do ! k 282 if (k==nspline+1) my_err=my_err+1 ! xfit not found 283 end do ! i 284 if (present(ierr)) ierr=my_err 285 286 end subroutine paw_splint
m_paw_numeric/paw_splint_der [ Functions ]
[ Top ] [ m_paw_numeric ] [ Functions ]
NAME
paw_splint_der
FUNCTION
Compute spline interpolation of the derivative of a tabulated function. There is no hypothesis about the spacing of the input grid points.
INPUTS
nspline: number of grid points of input mesh xspline(nspline): input mesh yspline(nspline): function on input mesh ysplin2(nspline): second derivative of yspline on input mesh nfit: number of points of output mesh xfit(nfit): output mesh
OUTPUT
dydxfit(nfit): 1st-derivative of function on output mesh [ierr]=A non-zero value is used to signal that some points in xfit exceed xspline(nspline). The input value is incremented by the number of such points.
SOURCE
314 subroutine paw_splint_der(nspline,xspline,yspline,ysplin2,nfit,xfit,dydxfit,ierr) 315 316 !Arguments ------------------------------------ 317 !scalars 318 integer,intent(in) :: nfit, nspline 319 integer,optional,intent(out) :: ierr 320 !arrays 321 real(dp),intent(in) :: xspline(nspline),yspline(nspline) 322 real(dp),intent(in) :: ysplin2(nspline),xfit(nfit) 323 real(dp),intent(out) :: dydxfit(nfit) 324 325 !Local variables------------------------------- 326 !scalars 327 integer :: left,i,k,right,my_err 328 real(dp) :: delarg,invdelarg,aa,bb 329 character(len=50) :: msg 330 !arrays 331 332 ! ************************************************************************* 333 334 my_err=0 335 left=1 336 do i=1,nfit 337 dydxfit(i)=0._dp ! Initialize for the unlikely event that rmax exceed r(mesh) 338 do k=left+1, nspline 339 if(xspline(k) >= xfit(i)) then 340 if(xspline(k-1) <= xfit(i)) then 341 right = k 342 left = k-1 343 else 344 if (k-1.eq.1 .and. i.eq.1) then 345 msg='xfit(1) < xspline(1)' 346 else 347 msg='xfit not properly ordered' 348 end if 349 LIBPAW_ERROR(msg) 350 end if 351 delarg= xspline(right) - xspline(left) 352 invdelarg= 1.0_dp/delarg 353 aa= (xspline(right)-xfit(i))*invdelarg 354 bb= (xfit(i)-xspline(left))*invdelarg 355 dydxfit(i) = (yspline(right)-yspline(left))*invdelarg & 356 & -( (3.0_dp*(aa*aa)-1.0_dp) *ysplin2(left) & 357 & -(3.0_dp*(bb*bb)-1.0_dp) *ysplin2(right) ) *delarg/6.0_dp 358 exit 359 end if 360 end do ! k 361 if (k==nspline+1) my_err=my_err+1 ! xfit not found 362 end do ! i 363 if (present(ierr)) ierr=my_err 364 365 end subroutine paw_splint_der
m_paw_numeric/paw_uniform_splfit [ Functions ]
[ Top ] [ m_paw_numeric ] [ Functions ]
NAME
paw_uniform_splfit
FUNCTION
Evaluate cubic spline fit to get function values on input set of ORDERED, UNIFORMLY SPACED points. Optionally gives derivatives (first and second) at those points too. If point lies outside the range of arg, assign the extremal point values to these points, and zero derivative.
INPUTS
arg(numarg)=equally spaced arguments (spacing delarg) for data to which spline was fit. fun(numarg,2)=function values to which spline was fit and spline fit to second derivatives (from Numerical Recipes spline). ider= see above newarg(numnew)=new values of arguments at which function is desired. numarg=number of arguments at which spline was fit. numnew=number of arguments at which function values are desired.
OUTPUT
derfun(numnew)=(optional) values of first or second derivative of function. This is only computed for ider=1 or 2; otherwise derfun not used. newfun(numnew)=values of function at newarg(numnew). This is only computed for ider=0 or 1.
NOTES
if ider=0, compute only the function (contained in fun) if ider=1, compute the function (contained in fun) and its first derivative (in derfun) if ider=2, compute only the second derivative of the function (in derfun)
SOURCE
404 subroutine paw_uniform_splfit(arg,derfun,fun,ider,newarg,newfun,numarg,numnew) 405 406 !Arguments ------------------------------------ 407 !scalars 408 integer, intent(in) :: ider,numarg,numnew 409 !arrays 410 real(dp), intent(in) :: arg(numarg),fun(numarg,2),newarg(numnew) 411 real(dp), intent(out) :: derfun(numnew) 412 real(dp), intent(inout) :: newfun(numnew) 413 414 !Local variables------------------------------- 415 !scalars 416 integer :: i,jspl 417 real(dp) :: argmin,delarg,d,aa,bb,cc,dd 418 character(len=500) :: msg 419 !arrays 420 421 ! ************************************************************************* 422 423 !argmin is smallest x value in spline fit; delarg is uniform spacing of spline argument 424 argmin=arg(1) 425 delarg=(arg(numarg)-argmin)/dble(numarg-1) 426 427 if(delarg<tol12)then 428 write(msg,'(a,es16.8)') 'delarg should be strictly positive, while delarg= ',delarg 429 LIBPAW_ERROR(msg) 430 endif 431 432 jspl=-1 433 434 !Do one loop for no grads, other for grads: 435 if (ider==0) then 436 437 ! Spline index loop for no grads: 438 do i=1,numnew 439 if (newarg(i).ge.arg(numarg)) then 440 newfun(i)=fun(numarg,1) 441 else if (newarg(i).le.arg(1)) then 442 newfun(i)=fun(1,1) 443 else 444 jspl=1+int((newarg(i)-argmin)/delarg) 445 d=newarg(i)-arg(jspl) 446 bb = d/delarg 447 aa = 1.0d0-bb 448 cc = aa*(aa**2-1.0d0)*(delarg**2/6.0d0) 449 dd = bb*(bb**2-1.0d0)*(delarg**2/6.0d0) 450 newfun(i)=aa*fun(jspl,1)+bb*fun(jspl+1,1)+cc*fun(jspl,2)+dd*fun(jspl+1,2) 451 end if 452 end do 453 454 else if(ider==1)then 455 456 ! Spline index loop includes grads: 457 do i=1,numnew 458 if (newarg(i).ge.arg(numarg)) then 459 newfun(i)=fun(numarg,1) 460 derfun(i)=0.0d0 461 else if (newarg(i).le.arg(1)) then 462 newfun(i)=fun(1,1) 463 derfun(i)=0.0d0 464 else 465 ! cubic spline interpolation: 466 jspl=1+int((newarg(i)-arg(1))/delarg) 467 d=newarg(i)-arg(jspl) 468 bb = d/delarg 469 aa = 1.0d0-bb 470 cc = aa*(aa**2-1.0d0)*(delarg**2/6.0d0) 471 dd = bb*(bb**2-1.0d0)*(delarg**2/6.0d0) 472 newfun(i)=aa*fun(jspl,1)+bb*fun(jspl+1,1)+cc*fun(jspl,2)+dd*fun(jspl+1,2) 473 ! spline fit to first derivative: 474 ! note correction of Numerical Recipes sign error 475 derfun(i) = (fun(jspl+1,1)-fun(jspl,1))/delarg + & 476 & (-(3.d0*aa**2-1.d0)*fun(jspl,2)+ & 477 & (3.d0*bb**2-1.d0)*fun(jspl+1,2)) * delarg/6.0d0 478 end if 479 end do 480 481 else if (ider==2) then 482 483 do i=1,numnew 484 if (newarg(i).ge.arg(numarg)) then 485 derfun(i)=0.0d0 486 else if (newarg(i).le.arg(1)) then 487 derfun(i)=0.0d0 488 else 489 ! cubic spline interpolation: 490 jspl=1+int((newarg(i)-argmin)/delarg) 491 d=newarg(i)-arg(jspl) 492 bb = d/delarg 493 aa = 1.0d0-bb 494 ! second derivative of spline (piecewise linear function) 495 derfun(i) = aa*fun(jspl,2)+bb*fun(jspl+1,2) 496 end if 497 end do 498 499 end if 500 501 end subroutine paw_uniform_splfit
m_special_funcs/paw_jbessel_4spline [ Functions ]
[ Top ] [ m_special_funcs ] [ Functions ]
NAME
paw_jbessel_4spline
FUNCTION
Compute spherical Bessel functions and derivatives. A polynomial approximation is employed for q-->0.
INPUTS
ll=l-order of the Bessel function tol=tolerance below which a Polynomial approximation is employed both for jl and its derivative (if required) order=1 if only first derivative is requested 2 if first and second derivatives are requested xx=where to compute j_l
OUTPUT
bes=Spherical Bessel function j_l at xx besp= first derivative of j_l at xx (only if order>=1)
TODO
Remove inline definitions, they are obsolete in F2003
SOURCE
887 subroutine paw_jbessel_4spline(bes,besp,ll,order,xx,tol) 888 889 !Arguments --------------------------------------------- 890 !scalars 891 integer,intent(in) :: ll,order 892 real(dp),intent(in) :: xx,tol 893 real(dp),intent(out) :: bes,besp 894 895 !Local variables --------------------------------------- 896 !scalars 897 real(dp) :: bespp 898 !real(dp) :: arg,bes0a,bes0ap,bes0b,bes0bp,bes1a,bes1ap,bes1b,bes1bp 899 !real(dp) :: bes2a,bes2ap,bes2b,bes2bp,bes3a,bes3ap,bes3b,bes3bp 900 character(len=100) :: msg 901 902 ! ********************************************************************* 903 904 ! === l=0,1,2 and 3 spherical Bessel functions (and derivatives) === 905 ! Statement functions are obsolete. Sorry ... 906 !bes0a(arg)=1.0_dp-arg**2/6.0_dp*(1.0_dp-arg**2/20.0_dp) 907 !bes0b(arg)=sin(arg)/arg 908 !bes1a(arg)=(10.0_dp-arg*arg)*arg/30.0_dp 909 !bes1b(arg)=(sin(arg)-arg*cos(arg))/arg**2 910 !bes2a(arg)=arg*arg/15.0_dp-arg**4/210.0_dp 911 !bes2b(arg)=((3.0_dp-arg**2)*sin(arg)-3.0_dp*arg*cos(arg))/arg**3 912 !bes3a(arg)=arg*arg*arg/105.0_dp-arg**5/1890.0_dp+arg**7/83160.0_dp 913 !bes3b(arg)=(15.0_dp*sin(arg)-15.0_dp*arg*cos(arg)-6.0_dp*arg**2*sin(arg)+arg**3*cos(arg))/arg**4 914 !bes0ap(arg)=(-10.0_dp+arg*arg)*arg/30.0_dp 915 !bes0bp(arg)=-(sin(arg)-arg*cos(arg))/arg**2 916 !bes1ap(arg)=(10.0_dp-3.0_dp*arg*arg)/30.0_dp 917 !bes1bp(arg)=((arg*arg-2.0_dp)*sin(arg)+2.0_dp*arg*cos(arg))/arg**3 918 !bes2ap(arg)=(1.0_dp-arg*arg/7.0_dp)*2.0_dp*arg/15.0_dp 919 !bes2bp(arg)=((4.0_dp*arg*arg-9.0_dp)*sin(arg)+(9.0_dp-arg*arg)*arg*cos(arg))/arg**4 920 !bes3ap(arg)=(1.0_dp/35-arg*arg/378.0_dp+arg**4/11880.0_dp)*arg*arg 921 !bes3bp(arg)=((-60.0_dp+27.0_dp*arg*arg-arg**4)*sin(arg)+(60.0_dp*arg-7.0_dp*arg**3)*cos(arg))/arg**5 922 923 ! This is to test paw_jbessel calculation without polynomial approximation for q-->0. 924 ! call paw_jbessel(bes,besp,bespp,ll,order,xx) 925 ! RETURN 926 927 if (order>2) then 928 msg='Wrong order in paw_jbessel_4spline' 929 LIBPAW_ERROR(msg) 930 end if 931 932 select case (ll) 933 case (0) 934 if (xx<TOL) then 935 bes=1.0_dp-xx**2/6.0_dp*(1.0_dp-xx**2/20.0_dp) 936 if (order>=1) besp=(-10.0_dp+xx*xx)*xx/30.0_dp 937 else 938 bes=sin(xx)/xx 939 if (order>=1) besp=-(sin(xx)-xx*cos(xx))/xx**2 940 end if 941 942 case (1) 943 if (xx<TOL) then 944 bes=(10.0_dp-xx*xx)*xx/30.0_dp 945 if (order>=1) besp=(10.0_dp-3.0_dp*xx*xx)/30.0_dp 946 else 947 bes=(sin(xx)-xx*cos(xx))/xx**2 948 if (order>=1) besp=((xx*xx-2.0_dp)*sin(xx)+2.0_dp*xx*cos(xx))/xx**3 949 end if 950 951 case (2) 952 if (xx<TOL) then 953 bes=xx*xx/15.0_dp-xx**4/210.0_dp 954 if (order>=1) besp=(1.0_dp-xx*xx/7.0_dp)*2.0_dp*xx/15.0_dp 955 else 956 bes=((3.0_dp-xx**2)*sin(xx)-3.0_dp*xx*cos(xx))/xx**3 957 if (order>=1) besp=((4.0_dp*xx*xx-9.0_dp)*sin(xx)+(9.0_dp-xx*xx)*xx*cos(xx))/xx**4 958 end if 959 960 case (3) 961 if (xx<TOL) then 962 bes=xx*xx*xx/105.0_dp-xx**5/1890.0_dp+xx**7/83160.0_dp 963 if (order>=1) besp=(1.0_dp/35-xx*xx/378.0_dp+xx**4/11880.0_dp)*xx*xx 964 else 965 bes=(15.0_dp*sin(xx)-15.0_dp*xx*cos(xx)-6.0_dp*xx**2*sin(xx)+xx**3*cos(xx))/xx**4 966 if (order>=1) besp=((-60.0_dp+27.0_dp*xx*xx-xx**4)*sin(xx)+(60.0_dp*xx-7.0_dp*xx**3)*cos(xx))/xx**5 967 end if 968 969 case (4:) 970 call paw_jbessel(bes,besp,bespp,ll,order,xx) 971 972 case default 973 write(msg,'(a,i4)')' wrong value for ll = ',ll 974 LIBPAW_BUG(msg) 975 end select 976 977 end subroutine paw_jbessel_4spline