TABLE OF CONTENTS


ABINIT/m_paw_numeric [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_numeric

FUNCTION

  Wrappers for various numeric operations (spline, sort, ...)

COPYRIGHT

  Copyright (C) 2012-2022 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