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-2018 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

21 #include "libpaw.h"
22 
23 module m_paw_numeric
24 
25  USE_DEFS
26  USE_MSG_HANDLING
27  USE_MEMORY_PROFILING
28 
29  implicit none
30 
31  private
32 
33 !public procedures
34  public:: paw_spline
35  public:: paw_splint
36  public:: paw_splint_der
37  public:: paw_uniform_splfit
38  public:: paw_smooth
39  public:: paw_sort_dp
40  public:: paw_jbessel
41  public:: paw_solvbes
42  public:: paw_jbessel_4spline
43  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

PARENTS

     screened_coul_kernel

CHILDREN

SOURCE

1143 elemental function paw_derfc(yy) result(derfc_yy)
1144 
1145 
1146 !This section has been created automatically by the script Abilint (TD).
1147 !Do not modify the following lines by hand.
1148 #undef ABI_FUNC
1149 #define ABI_FUNC 'paw_derfc'
1150 !End of the abilint section
1151 
1152  implicit none
1153 
1154 !Arguments ------------------------------------
1155 !scalars
1156  real(dp),intent(in) :: yy
1157  real(dp) :: derfc_yy
1158 
1159 !Local variables-------------------------------
1160  integer          ::  done,ii,isw
1161 ! coefficients for 0.0 <= yy < .477
1162  real(dp), parameter :: &
1163 &  pp(5)=(/ 113.8641541510502e0_dp, 377.4852376853020e0_dp,  &
1164 &           3209.377589138469e0_dp, .1857777061846032e0_dp,  &
1165 &           3.161123743870566e0_dp /)
1166  real(dp), parameter :: &
1167 &  qq(4)=(/ 244.0246379344442e0_dp, 1282.616526077372e0_dp,  &
1168 &           2844.236833439171e0_dp, 23.60129095234412e0_dp/)
1169 ! coefficients for .477 <= yy <= 4.0
1170  real(dp), parameter :: &
1171 &  p1(9)=(/ 8.883149794388376e0_dp, 66.11919063714163e0_dp,  &
1172 &           298.6351381974001e0_dp, 881.9522212417691e0_dp,  &
1173 &           1712.047612634071e0_dp, 2051.078377826071e0_dp,  &
1174 &           1230.339354797997e0_dp, 2.153115354744038e-8_dp, &
1175 &           .5641884969886701e0_dp /)
1176  real(dp), parameter :: &
1177 &  q1(8)=(/ 117.6939508913125e0_dp, 537.1811018620099e0_dp,  &
1178 &           1621.389574566690e0_dp, 3290.799235733460e0_dp,  &
1179 &           4362.619090143247e0_dp, 3439.367674143722e0_dp,  &
1180 &           1230.339354803749e0_dp, 15.74492611070983e0_dp/)
1181  ! coefficients for 4.0 < y,
1182  real(dp), parameter :: &
1183 &  p2(6)=(/ -3.603448999498044e-01_dp, -1.257817261112292e-01_dp,   &
1184 &           -1.608378514874228e-02_dp, -6.587491615298378e-04_dp,   &
1185 &           -1.631538713730210e-02_dp, -3.053266349612323e-01_dp/)
1186  real(dp), parameter :: &
1187 &  q2(5)=(/ 1.872952849923460e0_dp   , 5.279051029514284e-01_dp,    &
1188 &           6.051834131244132e-02_dp , 2.335204976268692e-03_dp,    &
1189 &           2.568520192289822e0_dp /)
1190  real(dp), parameter :: &
1191 &  sqrpi=.5641895835477563e0_dp, xbig=13.3e0_dp, xlarge=6.375e0_dp, xmin=1.0e-10_dp
1192  real(dp) ::  res,xden,xi,xnum,xsq,xx
1193 
1194 !******************************************************************
1195 
1196  xx = yy
1197  isw = 1
1198 !Here change the sign of xx, and keep track of it thanks to isw
1199  if (xx<0.0e0_dp) then
1200    isw = -1
1201    xx = -xx
1202  end if
1203 
1204  done=0
1205 
1206 !Residual value, if yy < -6.375e0_dp
1207  res=2.0e0_dp
1208 
1209 !abs(yy) < .477, evaluate approximation for erfc
1210  if (xx<0.477e0_dp) then
1211 !  xmin is a very small number
1212    if (xx<xmin) then
1213      res = xx*pp(3)/qq(3)
1214    else
1215      xsq = xx*xx
1216      xnum = pp(4)*xsq+pp(5)
1217      xden = xsq+qq(4)
1218      do ii = 1,3
1219        xnum = xnum*xsq+pp(ii)
1220        xden = xden*xsq+qq(ii)
1221      end do
1222      res = xx*xnum/xden
1223    end if
1224    if (isw==-1) res = -res
1225    res = 1.0e0_dp-res
1226    done=1
1227  end if
1228 
1229 !.477 < abs(yy) < 4.0 , evaluate approximation for erfc
1230  if (xx<=4.0e0_dp .and. done==0 ) then
1231    xsq = xx*xx
1232    xnum = p1(8)*xx+p1(9)
1233    xden = xx+q1(8)
1234    do ii=1,7
1235      xnum = xnum*xx+p1(ii)
1236      xden = xden*xx+q1(ii)
1237    end do
1238    res = xnum/xden
1239    res = res* exp(-xsq)
1240    if (isw.eq.-1) res = 2.0e0_dp-res
1241    done=1
1242  end if
1243 
1244 !y > 13.3e0_dp
1245  if (isw > 0 .and. xx > xbig .and. done==0 ) then
1246    res = 0.0e0_dp
1247    done=1
1248  end if
1249 
1250 !4.0 < yy < 13.3e0_dp  .or. -6.375e0_dp < yy < -4.0
1251 !evaluate minimax approximation for erfc
1252  if ( ( isw > 0 .or. xx < xlarge ) .and. done==0 ) then
1253    xsq = xx*xx
1254    xi = 1.0e0_dp/xsq
1255    xnum= p2(5)*xi+p2(6)
1256    xden = xi+q2(5)
1257    do ii = 1,4
1258      xnum = xnum*xi+p2(ii)
1259      xden = xden*xi+q2(ii)
1260    end do
1261    res = (sqrpi+xi*xnum/xden)/xx
1262    res = res* exp(-xsq)
1263    if (isw.eq.-1) res = 2.0e0_dp-res
1264  end if
1265 
1266 !All cases have been investigated
1267  derfc_yy = res
1268 
1269 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)

PARENTS

      m_paw_atom,m_paw_finegrid,m_paw_numeric,m_vcoul

CHILDREN

      paw_jbessel

SOURCE

782 subroutine paw_jbessel(bes,besp,bespp,ll,order,xx)
783 
784 
785 !This section has been created automatically by the script Abilint (TD).
786 !Do not modify the following lines by hand.
787 #undef ABI_FUNC
788 #define ABI_FUNC 'paw_jbessel'
789 !End of the abilint section
790 
791  implicit none
792 
793 !Arguments ---------------------------------------------
794 !scalars
795  integer,intent(in) :: ll,order
796  real(dp),intent(in) :: xx
797  real(dp),intent(out) :: bes,besp,bespp
798 
799 !Local variables ---------------------------------------
800 !scalars
801  integer,parameter :: imax=40
802  integer :: ii,il
803  real(dp),parameter :: prec=1.d-15
804  real(dp) :: besp1,fact,factp,factpp,jn,jnp,jnpp,jr,xx2,xxinv
805  character(len=200) :: msg
806 
807 ! *********************************************************************
808 
809  if (order>2) then
810    msg='Wrong order in paw_jbessel!'
811    MSG_ERROR(msg)
812  end if
813 
814  if (abs(xx)<prec) then
815    bes=zero;if (ll==0) bes=one
816    if (order>=1) then
817      besp=zero;if (ll==1) besp=third
818    end if
819    if (order==2) then
820      bespp=zero
821      if (ll==0) bespp=-third
822      if (ll==2) bespp=2._dp/15._dp
823    end if
824    return
825  end if
826 
827  xxinv=one/xx
828  if (order==0) then
829    factp=zero
830    factpp=zero
831    jnp=zero
832    jnpp=zero
833  end if
834 
835  if (xx<one) then
836    xx2=0.5_dp*xx*xx
837    fact=one
838    do il=1,ll
839      fact=fact*xx/dble(2*il+1)
840    end do
841    jn=one;jr=one;ii=0
842    do while(abs(jr)>=prec.and.ii<imax)
843      ii=ii+1;jr=-jr*xx2/dble(ii*(2*(ll+ii)+1))
844      jn=jn+jr
845    end do
846    bes=jn*fact
847    if (abs(jr)>prec) then
848      msg='Bessel function did not converge!'
849      MSG_ERROR(msg)
850    end if
851    if (order>=1) then
852      factp=fact*xx/dble(2*ll+3)
853      jnp=one;jr=one;ii=0
854      do while(abs(jr)>=prec.AND.ii<imax)
855        ii=ii+1;jr=-jr*xx2/dble(ii*(2*(ll+ii)+3))
856        jnp=jnp+jr
857      end do
858      besp=-jnp*factp+jn*fact*xxinv*dble(ll)
859      if (abs(jr)>prec) then
860        msg='1st der. of Bessel function did not converge!'
861        MSG_ERROR(msg)
862      end if
863    end if
864    if (order==2) then
865      factpp=factp*xx/dble(2*ll+5)
866      jnpp=one;jr=one;ii=0
867      do while(abs(jr)>=prec.AND.ii<imax)
868        ii=ii+1;jr=-jr*xx2/dble(ii*(2*(ll+ii)+5))
869        jnpp=jnpp+jr
870      end do
871      besp1=-jnpp*factpp+jnp*factp*xxinv*dble(ll+1)
872      if (abs(jr)>prec) then
873        msg='2nd der. of Bessel function did not converge !'
874        MSG_ERROR(msg)
875      end if
876    end if
877  else
878    jn =sin(xx)*xxinv
879    jnp=(-cos(xx)+jn)*xxinv
880    do il=2,ll+1
881      jr=-jn+dble(2*il-1)*jnp*xxinv
882      jn=jnp;jnp=jr
883    end do
884    bes=jn
885    if (order>=1) besp =-jnp+jn *xxinv*dble(ll)
886    if (order==2) besp1= jn -jnp*xxinv*dble(ll+2)
887  end if
888 
889  if (order==2) bespp=-besp1+besp*ll*xxinv-bes*ll*xxinv*xxinv
890 
891 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

PARENTS

      m_pawpsp

CHILDREN

      paw_jbessel

SOURCE

592 subroutine paw_smooth(a,mesh,it)
593 
594 
595 !This section has been created automatically by the script Abilint (TD).
596 !Do not modify the following lines by hand.
597 #undef ABI_FUNC
598 #define ABI_FUNC 'paw_smooth'
599 !End of the abilint section
600 
601  implicit none
602 
603 !Arguments ------------------------------------
604 !scalars
605  integer, intent(in) :: it,mesh
606 !arrays
607  real(dp), intent(inout) :: a(mesh)
608 
609 !Local variables-------------------------------
610 !scalars
611  integer :: i,k
612 !arrays
613  real(dp) :: asm(mesh)
614 
615 ! *************************************************************************
616 
617  asm(1:4) = zero ! ?? Correct me ...
618  do k=1,it
619    asm(5)=0.2_dp*(a(3)+a(4)+a(5)+a(6)+a(7))
620    asm(mesh-4)=0.2_dp*(a(mesh-2)+a(mesh-3)+a(mesh-4)+&
621 &                     a(mesh-5)+a(mesh-6))
622    asm(mesh-3)=0.2_dp*(a(mesh-1)+a(mesh-2)+a(mesh-3)+&
623 &                     a(mesh-4)+a(mesh-5))
624    asm(mesh-2)=0.2_dp*(a(mesh)+a(mesh-1)+a(mesh-2)+&
625 &                     a(mesh-3)+a(mesh-4))
626    asm(mesh-1)=0.25_dp*(a(mesh)+a(mesh-1)+a(mesh-2)+a(mesh-3))
627    asm(mesh)=1.0_dp/3.0_dp*(a(mesh)+a(mesh-1)+a(mesh-2))
628    do i=6,mesh-5
629      asm(i)=0.1_dp *a(i)+0.1_dp*(a(i+1)+a(i-1))+&
630 &           0.1_dp *(a(i+2)+a(i-2))+&
631 &           0.1_dp *(a(i+3)+a(i-3))+&
632 &           0.1_dp *(a(i+4)+a(i-4))+&
633 &           0.05_dp*(a(i+5)+a(i-5))
634    end do
635    do i=1,mesh
636      a(i)=asm(i)
637    end do
638  end do
639 
640 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

PARENTS

      m_paw_atom

CHILDREN

      paw_jbessel

SOURCE

920  subroutine paw_solvbes(root,alpha,beta,ll,nq)
921 
922 
923 !This section has been created automatically by the script Abilint (TD).
924 !Do not modify the following lines by hand.
925 #undef ABI_FUNC
926 #define ABI_FUNC 'paw_solvbes'
927 !End of the abilint section
928 
929  implicit none
930 
931 !Arguments ------------------------------------
932 !scalars
933  integer :: ll,nq
934  real(dp) :: alpha,beta
935 !arrays
936  real(dp) :: root(nq)
937 
938 !Local variables-------------------------------
939 !scalars
940  integer :: nroot
941  real(dp),parameter :: dh=0.1_dp,tol=tol14
942  real(dp) :: dum,hh,jbes,jbesp,qq,qx,y1,y2
943 
944 ! *************************************************************************
945 
946  qq=dh;nroot=0
947 
948  do while (nroot<nq)
949    call paw_jbessel(jbes,jbesp,dum,ll,1,qq)
950    y1=alpha*jbes+beta*qq*jbesp
951    qq=qq+dh
952    call paw_jbessel(jbes,jbesp,dum,ll,1,qq)
953    y2=alpha*jbes+beta*qq*jbesp
954 
955    do while (y1*y2>=zero)
956      qq=qq+dh
957      call paw_jbessel(jbes,jbesp,dum,ll,1,qq)
958      y2=alpha*jbes+beta*qq*jbesp
959    end do
960 
961    hh=dh;qx=qq
962    do while (hh>tol)
963      hh=half*hh
964      if (y1*y2<zero) then
965        qx=qx-hh
966      else
967        qx=qx+hh
968      end if
969      call paw_jbessel(jbes,jbesp,dum,ll,1,qx)
970      y2=alpha*jbes+beta*qx*jbesp
971    end do
972    nroot=nroot+1
973    root(nroot)=qx
974 
975  end do
976 
977 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

PARENTS

      m_paw_finegrid

CHILDREN

      paw_jbessel

SOURCE

673 subroutine paw_sort_dp(n,list,iperm,tol)
674 
675 
676 !This section has been created automatically by the script Abilint (TD).
677 !Do not modify the following lines by hand.
678 #undef ABI_FUNC
679 #define ABI_FUNC 'paw_sort_dp'
680 !End of the abilint section
681 
682  implicit none
683 
684 !Arguments ------------------------------------
685 !scalars
686  integer, intent(in) :: n
687  real(dp), intent(in) :: tol
688 !arrays
689  integer, intent(inout) :: iperm(n)
690  real(dp), intent(inout) :: list(n)
691 
692 !Local variables-------------------------------
693 !scalars
694  integer :: l,ir,iap,i,j
695  real(dp) :: ap
696  character(len=500) :: msg
697 !arrays
698 
699 ! *************************************************************************
700 
701 !Accomodate case of array of length 1: already sorted!
702  if (n==1) return
703 
704 !Should not call with n<1
705  if (n<1) then
706    write(msg,'(a,i12,2a)') &
707 &   'paw_sort_dp has been called with array length n=',n, ch10, &
708 &   ' having a value less than 1.  This is not allowed.'
709    MSG_ERROR(msg)
710  end if
711 
712 !Conduct the usual sort
713  l=n/2+1 ; ir=n
714  do ! Infinite do-loop
715    if (l>1) then
716      l=l-1
717      ap=list(l)
718      iap=iperm(l)
719    else ! l<=1
720      ap=list(ir)
721      iap=iperm(ir)
722      list(ir)=list(1)
723      iperm(ir)=iperm(1)
724      ir=ir-1
725      if (ir==1) then
726        list(1)=ap
727        iperm(1)=iap
728        exit   ! This is the end of this algorithm
729      end if
730    end if ! l>1
731    i=l
732    j=l+l
733    do while (j<=ir)
734      if (j<ir) then
735        if ( list(j)<list(j+1)-tol .or.  &
736 &          (list(j)<list(j+1)+tol.and.iperm(j)<iperm(j+1))) j=j+1
737      endif
738      if (ap<list(j)-tol.or.(ap<list(j)+tol.and.iap<iperm(j))) then
739        list(i)=list(j)
740        iperm(i)=iperm(j)
741        i=j
742        j=j+j
743      else
744        j=ir+1
745      end if
746    end do
747    list(i)=ap
748    iperm(i)=iap
749  end do ! End infinite do-loop
750 
751 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 ...

PARENTS

      dfpt_eltfrxc,m_paw_atom,m_paw_gaussfit,m_paw_pwaves_lmn,m_pawpsp
      m_pawpwij,m_pawxmlps,m_psps

CHILDREN

      paw_jbessel

SOURCE

 82 subroutine paw_spline(t,y,n,ybcbeg,ybcend,ypp)
 83 
 84 !*******************************************************************************
 85 !  Discussion:
 86 !    For data interpolation, the user must call SPLINE_CUBIC_SET to
 87 !    determine the second derivative data, passing in the data to be
 88 !    interpolated, and the desired boundary conditions.
 89 !    The data to be interpolated, plus the SPLINE_CUBIC_SET output,
 90 !    defines the spline.  The user may then call SPLINE_CUBIC_VAL to
 91 !    evaluate the spline at any point.
 92 !    The cubic spline is a piecewise cubic polynomial.  The intervals
 93 !    are determined by the "knots" or abscissas of the data to be
 94 !    interpolated.  The cubic spline has continous first and second
 95 !    derivatives over the entire interval of interpolation.
 96 !    For any point T in the interval T(IVAL), T(IVAL+1), the form of
 97 !    the spline is
 98 !      SPL(T) = A(IVAL)
 99 !             + B(IVAL) * ( T - T(IVAL) )
100 !             + C(IVAL) * ( T - T(IVAL) )**2
101 !             + D(IVAL) * ( T - T(IVAL) )**3
102 !    If we assume that we know the values Y(*) and YPP(*), which represent
103 !    the values and second derivatives of the spline at each knot, then
104 !    the coefficients can be computed as:
105 !      A(IVAL) = Y(IVAL)
106 !      B(IVAL) = ( Y(IVAL+1) - Y(IVAL) ) / ( T(IVAL+1) - T(IVAL) )
107 !        - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * ( T(IVAL+1) - T(IVAL) ) / 6
108 !      C(IVAL) = YPP(IVAL) / 2
109 !      D(IVAL) = ( YPP(IVAL+1) - YPP(IVAL) ) / ( 6 * ( T(IVAL+1) - T(IVAL) ) )
110 !    Since the first derivative of the spline is
111 !      SPL'(T) =     B(IVAL)
112 !              + 2 * C(IVAL) * ( T - T(IVAL) )
113 !              + 3 * D(IVAL) * ( T - T(IVAL) )**2,
114 !    the requirement that the first derivative be continuous at interior
115 !    knot I results in a total of N-2 equations, of the form:
116 !      B(IVAL-1) + 2 C(IVAL-1) * (T(IVAL)-T(IVAL-1))
117 !      + 3 * D(IVAL-1) * (T(IVAL) - T(IVAL-1))**2 = B(IVAL)
118 !    or, setting H(IVAL) = T(IVAL+1) - T(IVAL)
119 !      ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
120 !      - ( YPP(IVAL) + 2 * YPP(IVAL-1) ) * H(IVAL-1) / 6
121 !      + YPP(IVAL-1) * H(IVAL-1)
122 !      + ( YPP(IVAL) - YPP(IVAL-1) ) * H(IVAL-1) / 2
123 !      =
124 !      ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
125 !      - ( YPP(IVAL+1) + 2 * YPP(IVAL) ) * H(IVAL) / 6
126 !    or
127 !      YPP(IVAL-1) * H(IVAL-1) + 2 * YPP(IVAL) * ( H(IVAL-1) + H(IVAL) )
128 !      + YPP(IVAL) * H(IVAL)
129 !      =
130 !      6 * ( Y(IVAL+1) - Y(IVAL) ) / H(IVAL)
131 !    - 6 * ( Y(IVAL) - Y(IVAL-1) ) / H(IVAL-1)
132 !    Boundary conditions must be applied at the first and last knots.
133 !    The resulting tridiagonal system can be solved for the YPP values.
134 !
135 !  Author:
136 !    John Burkardt, modified by Xavier Gonze
137 
138 !This section has been created automatically by the script Abilint (TD).
139 !Do not modify the following lines by hand.
140 #undef ABI_FUNC
141 #define ABI_FUNC 'paw_spline'
142 !End of the abilint section
143 
144 implicit none
145 
146 !Arguments ------------------------------------
147 !scalars
148  integer,intent(in) :: n
149  real(dp),intent(in) :: ybcbeg,ybcend
150 !arrays
151  real(dp),intent(in) :: t(n),y(n)
152  real(dp),intent(out) :: ypp(n)
153 
154 !Local variables-------------------------------
155 !scalars
156  integer :: ibcbeg,ibcend,i,k
157  real(dp) :: ratio,pinv
158  character(len=500) :: msg
159 !arrays
160  real(dp),allocatable :: tmp(:)
161 
162 ! *************************************************************************
163 
164 !Check
165  if (n<=1) then
166    write(msg,'(6a,i8)') ch10, &
167 &   'SPLINE_CUBIC_SET - Fatal error!',ch10, &
168 &   '  The number of knots must be at least 2.',ch10, &
169 &   '  The input value of N = ', n
170    MSG_ERROR(msg)
171  end if
172 
173  LIBPAW_ALLOCATE(tmp,(n))
174 
175  do i=1,n-1
176    if (t(i)>=t(i+1)) then
177    write(msg,'(6a,i8,a,es19.12,2a,i8,a,es19.12)') ch10, &
178 &   'SPLINE_CUBIC_SET - Fatal error!',ch10, &
179 &   '  The knots must be strictly increasing, but',ch10, &
180 &   '  T(',  i,') = ', t(i), ch10, &
181 &   '  T(',i+1,') = ', t(i+1)
182    MSG_ERROR(msg)
183    end if
184  end do
185 
186  ibcbeg=1;if(ybcbeg>1.0d+30)ibcbeg=0
187  ibcend=1;if(ybcend>1.0d+30)ibcend=0
188 
189 !Set the first and last equations
190  if (ibcbeg==0) then
191    ypp(1) = 0._dp
192    tmp(1) = 0._dp
193  else if ( ibcbeg == 1 ) then
194    ypp(1) = -0.5_dp
195    tmp(1) = (3._dp/(t(2)-t(1)))*((y(2)-y(1))/(t(2)-t(1))-ybcbeg)
196  end if
197  if (ibcend==0) then
198    ypp(n) = 0._dp
199    tmp(n) = 0._dp
200  else if ( ibcend == 1 ) then
201    ypp(n) = 0.5_dp
202    tmp(n) = (3._dp/(t(n)-t(n-1)))*(ybcend-(y(n)-y(n-1))/(t(n)-t(n-1)))
203  end if
204 
205 !Set the intermediate equations
206  do i=2,n-1
207    ratio=(t(i)-t(i-1))/(t(i+1)-t(i-1))
208    pinv = 1.0_dp/(ratio*ypp(i-1) + 2.0_dp)
209    ypp(i) = (ratio-1.0_dp)*pinv
210    tmp(i)=(6.0_dp*((y(i+1)-y(i))/(t(i+1)-t(i))-(y(i)-y(i-1)) &
211 &        /(t(i)-t(i-1)))/(t(i+1)-t(i-1))-ratio*tmp(i-1))*pinv
212    if (abs(tmp(i))<1.d5*tiny(0._dp)) tmp(i)=0._dp
213  end do
214 
215 !Solve the equations
216  ypp(n) = (tmp(n)-ypp(n)*tmp(n-1))/(ypp(n)*ypp(n-1)+1.0_dp)
217  do k=n-1,1,-1
218    ypp(k)=ypp(k)*ypp(k+1)+tmp(k)
219  end do
220 
221  LIBPAW_DEALLOCATE(tmp)
222 
223 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.

PARENTS

      m_paw_atom,m_paw_finegrid,m_paw_gaussfit,m_paw_pwaves_lmn,m_pawpsp
      m_pawxmlps,mkcore,mkcore_inner,mkcore_wvl,mklocl_realspace

CHILDREN

      paw_jbessel

SOURCE

258 subroutine paw_splint(nspline,xspline,yspline,ysplin2,nfit,xfit,yfit,ierr)
259 
260 
261 !This section has been created automatically by the script Abilint (TD).
262 !Do not modify the following lines by hand.
263 #undef ABI_FUNC
264 #define ABI_FUNC 'paw_splint'
265 !End of the abilint section
266 
267  implicit none
268 
269 !Arguments ------------------------------------
270 !scalars
271  integer,intent(in) :: nfit, nspline
272  integer,optional,intent(out) :: ierr
273 !arrays
274  real(dp),intent(in) :: xspline(nspline),yspline(nspline)
275  real(dp),intent(in) :: ysplin2(nspline),xfit(nfit)
276  real(dp),intent(out) :: yfit(nfit)
277 
278 !Local variables-------------------------------
279 !scalars
280  integer :: left,i,k,right,my_err
281  real(dp) :: delarg,invdelarg,aa,bb
282  character(len=50) :: msg
283 !arrays
284 
285 ! *************************************************************************
286 
287  my_err=0
288  left=1
289  do i=1,nfit
290    yfit(i)=0._dp  ! Initialize for the unlikely event that rmax exceed r(mesh)
291    do k=left+1, nspline
292      if(xspline(k) >= xfit(i)) then
293        if(xspline(k-1) <= xfit(i)) then
294          right = k
295          left = k-1
296        else
297          if (k-1.eq.1 .and. i.eq.1) then
298            msg='xfit(1) < xspline(1)'
299          else
300            msg='xfit not properly ordered'
301          end if
302          MSG_ERROR(msg)
303        end if
304        delarg= xspline(right) - xspline(left)
305        invdelarg= 1.0_dp/delarg
306        aa= (xspline(right)-xfit(i))*invdelarg
307        bb= (xfit(i)-xspline(left))*invdelarg
308        yfit(i) = aa*yspline(left)+bb*yspline(right)    &
309 &               +( (aa*aa*aa-aa)*ysplin2(left) +         &
310 &                  (bb*bb*bb-bb)*ysplin2(right) ) *delarg*delarg/6.0_dp
311        exit
312      end if
313    end do ! k
314    if (k==nspline+1) my_err=my_err+1 ! xfit not found
315  end do ! i
316  if (present(ierr)) ierr=my_err
317 
318 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.

PARENTS

      mkcore_inner,mklocl_realspace,mklocl_wavelets

CHILDREN

      paw_jbessel

SOURCE

352 subroutine paw_splint_der(nspline,xspline,yspline,ysplin2,nfit,xfit,dydxfit,ierr)
353 
354 
355 !This section has been created automatically by the script Abilint (TD).
356 !Do not modify the following lines by hand.
357 #undef ABI_FUNC
358 #define ABI_FUNC 'paw_splint_der'
359 !End of the abilint section
360 
361  implicit none
362 
363 !Arguments ------------------------------------
364 !scalars
365  integer,intent(in) :: nfit, nspline
366  integer,optional,intent(out) :: ierr
367 !arrays
368  real(dp),intent(in) :: xspline(nspline),yspline(nspline)
369  real(dp),intent(in) :: ysplin2(nspline),xfit(nfit)
370  real(dp),intent(out) :: dydxfit(nfit)
371 
372 !Local variables-------------------------------
373 !scalars
374  integer :: left,i,k,right,my_err
375  real(dp) :: delarg,invdelarg,aa,bb
376  character(len=50) :: msg
377 !arrays
378 
379 ! *************************************************************************
380 
381  my_err=0
382  left=1
383  do i=1,nfit
384    dydxfit(i)=0._dp  ! Initialize for the unlikely event that rmax exceed r(mesh)
385    do k=left+1, nspline
386      if(xspline(k) >= xfit(i)) then
387        if(xspline(k-1) <= xfit(i)) then
388          right = k
389          left = k-1
390        else
391          if (k-1.eq.1 .and. i.eq.1) then
392            msg='xfit(1) < xspline(1)'
393          else
394            msg='xfit not properly ordered'
395          end if
396          MSG_ERROR(msg)
397        end if
398        delarg= xspline(right) - xspline(left)
399        invdelarg= 1.0_dp/delarg
400        aa= (xspline(right)-xfit(i))*invdelarg
401        bb= (xfit(i)-xspline(left))*invdelarg
402        dydxfit(i) = (yspline(right)-yspline(left))*invdelarg &
403 &                  -( (3.0_dp*(aa*aa)-1.0_dp) *ysplin2(left) &
404 &                    -(3.0_dp*(bb*bb)-1.0_dp) *ysplin2(right) ) *delarg/6.0_dp
405        exit
406      end if
407    end do ! k
408    if (k==nspline+1) my_err=my_err+1 ! xfit not found
409  end do ! i
410  if (present(ierr)) ierr=my_err
411 
412 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)

PARENTS

CHILDREN

SOURCE

455 subroutine paw_uniform_splfit(arg,derfun,fun,ider,newarg,newfun,numarg,numnew)
456 
457 
458 !This section has been created automatically by the script Abilint (TD).
459 !Do not modify the following lines by hand.
460 #undef ABI_FUNC
461 #define ABI_FUNC 'paw_uniform_splfit'
462 !End of the abilint section
463 
464  implicit none
465 
466 !Arguments ------------------------------------
467 !scalars
468  integer, intent(in) :: ider,numarg,numnew
469 !arrays
470  real(dp), intent(in) :: arg(numarg),fun(numarg,2),newarg(numnew)
471  real(dp), intent(out) :: derfun(numnew)
472  real(dp), intent(inout) :: newfun(numnew)
473 
474 !Local variables-------------------------------
475 !scalars
476  integer :: i,jspl
477  real(dp) :: argmin,delarg,d,aa,bb,cc,dd
478  character(len=500) :: msg
479 !arrays
480 
481 ! *************************************************************************
482 
483 !argmin is smallest x value in spline fit; delarg is uniform spacing of spline argument
484  argmin=arg(1)
485  delarg=(arg(numarg)-argmin)/dble(numarg-1)
486 
487  if(delarg<tol12)then
488    write(msg,'(a,es16.8)') 'delarg should be strictly positive, while delarg= ',delarg
489    MSG_ERROR(msg)
490  endif
491 
492  jspl=-1
493 
494 !Do one loop for no grads, other for grads:
495  if (ider==0) then
496 
497 ! Spline index loop for no grads:
498   do i=1,numnew
499    if (newarg(i).ge.arg(numarg)) then
500     newfun(i)=fun(numarg,1)
501    else if (newarg(i).le.arg(1)) then
502     newfun(i)=fun(1,1)
503    else
504     jspl=1+int((newarg(i)-argmin)/delarg)
505     d=newarg(i)-arg(jspl)
506     bb = d/delarg
507     aa = 1.0d0-bb
508     cc = aa*(aa**2-1.0d0)*(delarg**2/6.0d0)
509     dd = bb*(bb**2-1.0d0)*(delarg**2/6.0d0)
510     newfun(i)=aa*fun(jspl,1)+bb*fun(jspl+1,1)+cc*fun(jspl,2)+dd*fun(jspl+1,2)
511    end if
512   end do
513 
514  else if(ider==1)then
515 
516 ! Spline index loop includes grads:
517   do i=1,numnew
518    if (newarg(i).ge.arg(numarg)) then
519     newfun(i)=fun(numarg,1)
520     derfun(i)=0.0d0
521    else if (newarg(i).le.arg(1)) then
522     newfun(i)=fun(1,1)
523     derfun(i)=0.0d0
524    else
525 !   cubic spline interpolation:
526     jspl=1+int((newarg(i)-arg(1))/delarg)
527     d=newarg(i)-arg(jspl)
528     bb = d/delarg
529     aa = 1.0d0-bb
530     cc = aa*(aa**2-1.0d0)*(delarg**2/6.0d0)
531     dd = bb*(bb**2-1.0d0)*(delarg**2/6.0d0)
532     newfun(i)=aa*fun(jspl,1)+bb*fun(jspl+1,1)+cc*fun(jspl,2)+dd*fun(jspl+1,2)
533 !   spline fit to first derivative:
534 !   note correction of Numerical Recipes sign error
535     derfun(i) = (fun(jspl+1,1)-fun(jspl,1))/delarg +    &
536 &      (-(3.d0*aa**2-1.d0)*fun(jspl,2)+                 &
537 &        (3.d0*bb**2-1.d0)*fun(jspl+1,2)) * delarg/6.0d0
538     end if
539   end do
540 
541  else if (ider==2) then
542 
543   do i=1,numnew
544    if (newarg(i).ge.arg(numarg)) then
545     derfun(i)=0.0d0
546    else if (newarg(i).le.arg(1)) then
547     derfun(i)=0.0d0
548    else
549 !   cubic spline interpolation:
550     jspl=1+int((newarg(i)-argmin)/delarg)
551     d=newarg(i)-arg(jspl)
552     bb = d/delarg
553     aa = 1.0d0-bb
554 !   second derivative of spline (piecewise linear function)
555     derfun(i) = aa*fun(jspl,2)+bb*fun(jspl+1,2)
556    end if
557   end do
558 
559  end if
560 
561 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

PARENTS

      m_pawpsp,m_pawpwij

CHILDREN

      paw_jbessel

SOURCE

1013 subroutine paw_jbessel_4spline(bes,besp,ll,order,xx,tol)
1014 
1015 
1016 !This section has been created automatically by the script Abilint (TD).
1017 !Do not modify the following lines by hand.
1018 #undef ABI_FUNC
1019 #define ABI_FUNC 'paw_jbessel_4spline'
1020 !End of the abilint section
1021 
1022  implicit none
1023 
1024 !This section has been created automatically by the script Abilint (TD).
1025 !Do not modify the following lines by hand.
1026 #undef ABI_FUNC
1027 #define ABI_FUNC 'paw_jbessel_4spline'
1028 !End of the abilint section
1029 
1030 !Arguments ---------------------------------------------
1031 !scalars
1032  integer,intent(in) :: ll,order
1033  real(dp),intent(in) :: xx,tol
1034  real(dp),intent(out) :: bes,besp
1035 
1036 !Local variables ---------------------------------------
1037 !scalars
1038  real(dp) :: bespp
1039 !real(dp) :: arg,bes0a,bes0ap,bes0b,bes0bp,bes1a,bes1ap,bes1b,bes1bp
1040 !real(dp) :: bes2a,bes2ap,bes2b,bes2bp,bes3a,bes3ap,bes3b,bes3bp
1041  character(len=100) :: msg
1042 
1043 ! *********************************************************************
1044 
1045 ! === l=0,1,2 and 3 spherical Bessel functions (and derivatives) ===
1046 ! Statement functions are obsolete. Sorry ...
1047 !bes0a(arg)=1.0_dp-arg**2/6.0_dp*(1.0_dp-arg**2/20.0_dp)
1048 !bes0b(arg)=sin(arg)/arg
1049 !bes1a(arg)=(10.0_dp-arg*arg)*arg/30.0_dp
1050 !bes1b(arg)=(sin(arg)-arg*cos(arg))/arg**2
1051 !bes2a(arg)=arg*arg/15.0_dp-arg**4/210.0_dp
1052 !bes2b(arg)=((3.0_dp-arg**2)*sin(arg)-3.0_dp*arg*cos(arg))/arg**3
1053 !bes3a(arg)=arg*arg*arg/105.0_dp-arg**5/1890.0_dp+arg**7/83160.0_dp
1054 !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
1055 !bes0ap(arg)=(-10.0_dp+arg*arg)*arg/30.0_dp
1056 !bes0bp(arg)=-(sin(arg)-arg*cos(arg))/arg**2
1057 !bes1ap(arg)=(10.0_dp-3.0_dp*arg*arg)/30.0_dp
1058 !bes1bp(arg)=((arg*arg-2.0_dp)*sin(arg)+2.0_dp*arg*cos(arg))/arg**3
1059 !bes2ap(arg)=(1.0_dp-arg*arg/7.0_dp)*2.0_dp*arg/15.0_dp
1060 !bes2bp(arg)=((4.0_dp*arg*arg-9.0_dp)*sin(arg)+(9.0_dp-arg*arg)*arg*cos(arg))/arg**4
1061 !bes3ap(arg)=(1.0_dp/35-arg*arg/378.0_dp+arg**4/11880.0_dp)*arg*arg
1062 !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
1063 
1064  ! This is to test paw_jbessel calculation without polynomial approximation for q-->0.
1065  ! call paw_jbessel(bes,besp,bespp,ll,order,xx)
1066  ! RETURN
1067 
1068  if (order>2) then
1069    msg='Wrong order in paw_jbessel_4spline'
1070    MSG_ERROR(msg)
1071  end if
1072 
1073  select case (ll)
1074  case (0)
1075    if (xx<TOL) then
1076      bes=1.0_dp-xx**2/6.0_dp*(1.0_dp-xx**2/20.0_dp)
1077      if (order>=1) besp=(-10.0_dp+xx*xx)*xx/30.0_dp
1078    else
1079      bes=sin(xx)/xx
1080      if (order>=1) besp=-(sin(xx)-xx*cos(xx))/xx**2
1081    end if
1082 
1083  case (1)
1084   if (xx<TOL) then
1085     bes=(10.0_dp-xx*xx)*xx/30.0_dp
1086     if (order>=1) besp=(10.0_dp-3.0_dp*xx*xx)/30.0_dp
1087   else
1088     bes=(sin(xx)-xx*cos(xx))/xx**2
1089     if (order>=1) besp=((xx*xx-2.0_dp)*sin(xx)+2.0_dp*xx*cos(xx))/xx**3
1090   end if
1091 
1092  case (2)
1093    if (xx<TOL) then
1094      bes=xx*xx/15.0_dp-xx**4/210.0_dp
1095      if (order>=1) besp=(1.0_dp-xx*xx/7.0_dp)*2.0_dp*xx/15.0_dp
1096    else
1097      bes=((3.0_dp-xx**2)*sin(xx)-3.0_dp*xx*cos(xx))/xx**3
1098      if (order>=1) besp=((4.0_dp*xx*xx-9.0_dp)*sin(xx)+(9.0_dp-xx*xx)*xx*cos(xx))/xx**4
1099    end if
1100 
1101  case (3)
1102    if (xx<TOL) then
1103      bes=xx*xx*xx/105.0_dp-xx**5/1890.0_dp+xx**7/83160.0_dp
1104      if (order>=1) besp=(1.0_dp/35-xx*xx/378.0_dp+xx**4/11880.0_dp)*xx*xx
1105    else
1106      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
1107      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
1108    end if
1109 
1110  case (4:)
1111    call paw_jbessel(bes,besp,bespp,ll,order,xx)
1112 
1113  case default
1114    write(msg,'(a,i4)')' wrong value for ll = ',ll
1115    MSG_BUG(msg)
1116  end select
1117 
1118 end subroutine paw_jbessel_4spline