Branch data Line data Source code
1 : : !!****m* ABINIT/m_paw_numeric
2 : : !! NAME
3 : : !! m_paw_numeric
4 : : !!
5 : : !! FUNCTION
6 : : !! Wrappers for various numeric operations (spline, sort, ...)
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2012-2026 ABINIT group (MT,TR)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! NOTES
15 : : !! FOR DEVELOPPERS: in order to preserve the portability of libPAW library,
16 : : !! please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
17 : : !!
18 : : !! SOURCE
19 : :
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
43 : : !!***
44 : :
45 : : CONTAINS
46 : : !===========================================================
47 : : !!***
48 : :
49 : : !----------------------------------------------------------------------
50 : :
51 : : !!****f* m_paw_numeric/paw_spline
52 : : !! NAME
53 : : !! paw_spline
54 : : !!
55 : : !! FUNCTION
56 : : !! Computes the second derivatives of a cubic spline
57 : : !!
58 : : !! INPUTS
59 : : !! * Input, integer N, the number of data points; N must be at least 2.
60 : : !! In the special case where N = 2 and IBCBEG = IBCEND = 0, the
61 : : !! spline will actually be linear.
62 : : !! * Input, real(dp) T(N), the knot values, that is, the points where data
63 : : !! is specified. The knot values should be distinct, and increasing.
64 : : !! * Input, real(dp) Y(N), the data values to be interpolated.
65 : : !! * Input, real(dp) YBCBEG, YBCEND, the values to be used in the boundary
66 : : !! conditions if IBCBEG or IBCEND is equal to 1 or 2.
67 : : !!
68 : : !! OUTPUT
69 : : !! Output, real(dp) YPP(N), the second derivatives of the cubic spline.
70 : : !! Work space, real(dp) DIAG(N) - should be removed ...
71 : : !!
72 : : !! SOURCE
73 : :
74 : 7127 : 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 : 7127 : real(dp),allocatable :: tmp(:)
145 : :
146 : : ! *************************************************************************
147 : :
148 : : !Check
149 [ - + ]: 7127 : if (n<=1) then
150 : 0 : write(msg,'(6a,i8)') ch10, &
151 : 0 : & 'SPLINE_CUBIC_SET - Fatal error!',ch10, &
152 : 0 : & ' The number of knots must be at least 2.',ch10, &
153 : 0 : & ' The input value of N = ', n
154 : 0 : LIBPAW_ERROR(msg)
155 : : end if
156 : :
157 [ - + + - : 21381 : LIBPAW_ALLOCATE(tmp,(n))
+ - + - -
+ - + ]
158 : :
159 [ + + ]: 21529833 : do i=1,n-1
160 [ - + ]: 21529833 : if (t(i)>=t(i+1)) then
161 : 0 : write(msg,'(6a,i8,a,es19.12,2a,i8,a,es19.12)') ch10, &
162 : 0 : & 'SPLINE_CUBIC_SET - Fatal error!',ch10, &
163 : 0 : & ' The knots must be strictly increasing, but',ch10, &
164 : 0 : & ' T(', i,') = ', t(i), ch10, &
165 : 0 : & ' T(',i+1,') = ', t(i+1)
166 : 0 : LIBPAW_ERROR(msg)
167 : : end if
168 : : end do
169 : :
170 [ - + ]: 7127 : ibcbeg=1;if(ybcbeg>1.0d+30)ibcbeg=0
171 [ - + ]: 7127 : ibcend=1;if(ybcend>1.0d+30)ibcend=0
172 : :
173 : : !Set the first and last equations
174 [ - + ]: 7127 : if (ibcbeg==0) then
175 : 0 : ypp(1) = 0._dp
176 : 0 : tmp(1) = 0._dp
177 : : else if ( ibcbeg == 1 ) then
178 : 7127 : ypp(1) = -0.5_dp
179 : 7127 : tmp(1) = (3._dp/(t(2)-t(1)))*((y(2)-y(1))/(t(2)-t(1))-ybcbeg)
180 : : end if
181 [ - + ]: 7127 : if (ibcend==0) then
182 : 0 : ypp(n) = 0._dp
183 : 0 : tmp(n) = 0._dp
184 : : else if ( ibcend == 1 ) then
185 : 7127 : ypp(n) = 0.5_dp
186 : 7127 : 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 [ + + ]: 21522706 : do i=2,n-1
191 : 21515579 : ratio=(t(i)-t(i-1))/(t(i+1)-t(i-1))
192 : 21515579 : pinv = 1.0_dp/(ratio*ypp(i-1) + 2.0_dp)
193 : 21515579 : 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 : 21515579 : & /(t(i)-t(i-1)))/(t(i+1)-t(i-1))-ratio*tmp(i-1))*pinv
196 [ + + ]: 21522706 : if (abs(tmp(i))<1.d5*tiny(0._dp)) tmp(i)=0._dp
197 : : end do
198 : :
199 : : !Solve the equations
200 : 7127 : ypp(n) = (tmp(n)-ypp(n)*tmp(n-1))/(ypp(n)*ypp(n-1)+1.0_dp)
201 [ + + ]: 21529833 : do k=n-1,1,-1
202 : 21529833 : ypp(k)=ypp(k)*ypp(k+1)+tmp(k)
203 : : end do
204 : :
205 : 7127 : LIBPAW_DEALLOCATE(tmp)
206 : :
207 : 7127 : end subroutine paw_spline
208 : : !!***
209 : :
210 : : !----------------------------------------------------------------------
211 : :
212 : : !!****f* m_paw_numeric/paw_splint
213 : : !! NAME
214 : : !! paw_splint
215 : : !!
216 : : !! FUNCTION
217 : : !! Compute spline interpolation of a tabulated function.
218 : : !! There is no hypothesis about the spacing of the input grid points.
219 : : !!
220 : : !! INPUTS
221 : : !! nspline: number of grid points of input mesh
222 : : !! xspline(nspline): input mesh
223 : : !! yspline(nspline): function on input mesh
224 : : !! ysplin2(nspline): second derivative of yspline on input mesh
225 : : !! nfit: number of points of output mesh
226 : : !! xfit(nfit): output mesh
227 : : !!
228 : : !! OUTPUT
229 : : !! yfit(nfit): function on output mesh
230 : : !! [ierr]=A non-zero value is used to signal that some points in xfit exceed xspline(nspline).
231 : : !! The input value is incremented by the number of such points.
232 : : !!
233 : : !! SOURCE
234 : :
235 : 73 : 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 : 73 : my_err=0
256 : 73 : left=1
257 [ + + ]: 122804 : do i=1,nfit
258 : 122731 : yfit(i)=0._dp ! Initialize for the unlikely event that rmax exceed r(mesh)
259 [ + - ]: 260378 : do k=left+1, nspline
260 [ + + ]: 260378 : if(xspline(k) >= xfit(i)) then
261 [ - + ]: 122731 : if(xspline(k-1) <= xfit(i)) then
262 : : right = k
263 : : left = k-1
264 : : else
265 [ # # ]: 0 : if (k-1.eq.1 .and. i.eq.1) then
266 : 0 : msg='xfit(1) < xspline(1)'
267 : : else
268 : 0 : msg='xfit not properly ordered'
269 : : end if
270 : 0 : LIBPAW_ERROR(msg)
271 : : end if
272 : 122731 : delarg= xspline(right) - xspline(left)
273 : 122731 : invdelarg= 1.0_dp/delarg
274 : 122731 : aa= (xspline(right)-xfit(i))*invdelarg
275 : 122731 : bb= (xfit(i)-xspline(left))*invdelarg
276 : : yfit(i) = aa*yspline(left)+bb*yspline(right) &
277 : : & +( (aa*aa*aa-aa)*ysplin2(left) + &
278 : 122731 : & (bb*bb*bb-bb)*ysplin2(right) ) *delarg*delarg/6.0_dp
279 : 122731 : exit
280 : : end if
281 : : end do ! k
282 [ - + ]: 122804 : if (k==nspline+1) my_err=my_err+1 ! xfit not found
283 : : end do ! i
284 [ - + ]: 73 : if (present(ierr)) ierr=my_err
285 : :
286 : 73 : end subroutine paw_splint
287 : : !!***
288 : :
289 : : !----------------------------------------------------------------------
290 : :
291 : : !!****f* m_paw_numeric/paw_splint_der
292 : : !! NAME
293 : : !! paw_splint_der
294 : : !!
295 : : !! FUNCTION
296 : : !! Compute spline interpolation of the derivative of a tabulated function.
297 : : !! There is no hypothesis about the spacing of the input grid points.
298 : : !!
299 : : !! INPUTS
300 : : !! nspline: number of grid points of input mesh
301 : : !! xspline(nspline): input mesh
302 : : !! yspline(nspline): function on input mesh
303 : : !! ysplin2(nspline): second derivative of yspline on input mesh
304 : : !! nfit: number of points of output mesh
305 : : !! xfit(nfit): output mesh
306 : : !!
307 : : !! OUTPUT
308 : : !! dydxfit(nfit): 1st-derivative of function on output mesh
309 : : !! [ierr]=A non-zero value is used to signal that some points in xfit exceed xspline(nspline).
310 : : !! The input value is incremented by the number of such points.
311 : : !!
312 : : !! SOURCE
313 : :
314 : 0 : 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 : 0 : my_err=0
335 : 0 : left=1
336 [ # # ]: 0 : do i=1,nfit
337 : 0 : dydxfit(i)=0._dp ! Initialize for the unlikely event that rmax exceed r(mesh)
338 [ # # ]: 0 : do k=left+1, nspline
339 [ # # ]: 0 : if(xspline(k) >= xfit(i)) then
340 [ # # ]: 0 : if(xspline(k-1) <= xfit(i)) then
341 : : right = k
342 : : left = k-1
343 : : else
344 [ # # ]: 0 : if (k-1.eq.1 .and. i.eq.1) then
345 : 0 : msg='xfit(1) < xspline(1)'
346 : : else
347 : 0 : msg='xfit not properly ordered'
348 : : end if
349 : 0 : LIBPAW_ERROR(msg)
350 : : end if
351 : 0 : delarg= xspline(right) - xspline(left)
352 : 0 : invdelarg= 1.0_dp/delarg
353 : 0 : aa= (xspline(right)-xfit(i))*invdelarg
354 : 0 : 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 : 0 : & -(3.0_dp*(bb*bb)-1.0_dp) *ysplin2(right) ) *delarg/6.0_dp
358 : 0 : exit
359 : : end if
360 : : end do ! k
361 [ # # ]: 0 : if (k==nspline+1) my_err=my_err+1 ! xfit not found
362 : : end do ! i
363 [ # # ]: 0 : if (present(ierr)) ierr=my_err
364 : :
365 : 0 : end subroutine paw_splint_der
366 : : !!***
367 : :
368 : : !----------------------------------------------------------------------
369 : :
370 : : !!****f* m_paw_numeric/paw_uniform_splfit
371 : : !! NAME
372 : : !! paw_uniform_splfit
373 : : !!
374 : : !! FUNCTION
375 : : !! Evaluate cubic spline fit to get function values on input set
376 : : !! of ORDERED, UNIFORMLY SPACED points.
377 : : !! Optionally gives derivatives (first and second) at those points too.
378 : : !! If point lies outside the range of arg, assign the extremal
379 : : !! point values to these points, and zero derivative.
380 : : !!
381 : : !! INPUTS
382 : : !! arg(numarg)=equally spaced arguments (spacing delarg) for data
383 : : !! to which spline was fit.
384 : : !! fun(numarg,2)=function values to which spline was fit and spline
385 : : !! fit to second derivatives (from Numerical Recipes spline).
386 : : !! ider= see above
387 : : !! newarg(numnew)=new values of arguments at which function is desired.
388 : : !! numarg=number of arguments at which spline was fit.
389 : : !! numnew=number of arguments at which function values are desired.
390 : : !!
391 : : !! OUTPUT
392 : : !! derfun(numnew)=(optional) values of first or second derivative of function.
393 : : !! This is only computed for ider=1 or 2; otherwise derfun not used.
394 : : !! newfun(numnew)=values of function at newarg(numnew).
395 : : !! This is only computed for ider=0 or 1.
396 : : !!
397 : : !! NOTES
398 : : !! if ider=0, compute only the function (contained in fun)
399 : : !! if ider=1, compute the function (contained in fun) and its first derivative (in derfun)
400 : : !! if ider=2, compute only the second derivative of the function (in derfun)
401 : : !!
402 : : !! SOURCE
403 : :
404 : 0 : 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 : 0 : argmin=arg(1)
425 : 0 : delarg=(arg(numarg)-argmin)/dble(numarg-1)
426 : :
427 [ # # ]: 0 : if(delarg<tol12)then
428 : 0 : write(msg,'(a,es16.8)') 'delarg should be strictly positive, while delarg= ',delarg
429 : 0 : LIBPAW_ERROR(msg)
430 : : endif
431 : :
432 : 0 : jspl=-1
433 : :
434 : : !Do one loop for no grads, other for grads:
435 [ # # ]: 0 : if (ider==0) then
436 : :
437 : : ! Spline index loop for no grads:
438 [ # # ]: 0 : do i=1,numnew
439 [ # # ]: 0 : if (newarg(i).ge.arg(numarg)) then
440 : 0 : newfun(i)=fun(numarg,1)
441 [ # # ]: 0 : else if (newarg(i).le.arg(1)) then
442 : 0 : newfun(i)=fun(1,1)
443 : : else
444 : 0 : jspl=1+int((newarg(i)-argmin)/delarg)
445 : 0 : d=newarg(i)-arg(jspl)
446 : 0 : bb = d/delarg
447 : 0 : aa = 1.0d0-bb
448 : 0 : cc = aa*(aa**2-1.0d0)*(delarg**2/6.0d0)
449 : 0 : dd = bb*(bb**2-1.0d0)*(delarg**2/6.0d0)
450 : 0 : 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 [ # # ]: 0 : else if(ider==1)then
455 : :
456 : : ! Spline index loop includes grads:
457 [ # # ]: 0 : do i=1,numnew
458 [ # # ]: 0 : if (newarg(i).ge.arg(numarg)) then
459 : 0 : newfun(i)=fun(numarg,1)
460 : 0 : derfun(i)=0.0d0
461 [ # # ]: 0 : else if (newarg(i).le.arg(1)) then
462 : 0 : newfun(i)=fun(1,1)
463 : 0 : derfun(i)=0.0d0
464 : : else
465 : : ! cubic spline interpolation:
466 : 0 : jspl=1+int((newarg(i)-arg(1))/delarg)
467 : 0 : d=newarg(i)-arg(jspl)
468 : 0 : bb = d/delarg
469 : 0 : aa = 1.0d0-bb
470 : 0 : cc = aa*(aa**2-1.0d0)*(delarg**2/6.0d0)
471 : 0 : dd = bb*(bb**2-1.0d0)*(delarg**2/6.0d0)
472 : 0 : 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 : 0 : & (3.d0*bb**2-1.d0)*fun(jspl+1,2)) * delarg/6.0d0
478 : : end if
479 : : end do
480 : :
481 [ # # ]: 0 : else if (ider==2) then
482 : :
483 [ # # ]: 0 : do i=1,numnew
484 [ # # ]: 0 : if (newarg(i).ge.arg(numarg)) then
485 : 0 : derfun(i)=0.0d0
486 [ # # ]: 0 : else if (newarg(i).le.arg(1)) then
487 : 0 : derfun(i)=0.0d0
488 : : else
489 : : ! cubic spline interpolation:
490 : 0 : jspl=1+int((newarg(i)-argmin)/delarg)
491 : 0 : d=newarg(i)-arg(jspl)
492 : 0 : bb = d/delarg
493 : 0 : aa = 1.0d0-bb
494 : : ! second derivative of spline (piecewise linear function)
495 : 0 : derfun(i) = aa*fun(jspl,2)+bb*fun(jspl+1,2)
496 : : end if
497 : : end do
498 : :
499 : : end if
500 : :
501 : 0 : end subroutine paw_uniform_splfit
502 : : !!***
503 : :
504 : : !----------------------------------------------------------------------
505 : :
506 : : !!****f* m_paw_numeric/paw_smooth
507 : : !! NAME
508 : : !! paw_smooth
509 : : !!
510 : : !! FUNCTION
511 : : !! Smooth an array of given ordinates (y's) that are in order of
512 : : !! increasing abscissas (x's), but without using the abscissas themselves
513 : : !! supposed to be equally spaced.
514 : : !!
515 : : !! INPUTS
516 : : !! it=number of abscissas to treat
517 : : !! mesh=size of the array (number of abscissas)
518 : : !!
519 : : !! OUTPUT
520 : : !!
521 : : !! SIDE EFFECTS
522 : : !! a(mesh)=array to be smoothed
523 : : !!
524 : : !! SOURCE
525 : :
526 : 0 : 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 : 0 : real(dp) :: asm(mesh)
539 : :
540 : : ! *************************************************************************
541 : :
542 [ # # ]: 0 : asm(1:4) = zero ! ?? Correct me ...
543 [ # # ]: 0 : do k=1,it
544 : 0 : 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 : 0 : & a(mesh-5)+a(mesh-6))
547 : : asm(mesh-3)=0.2_dp*(a(mesh-1)+a(mesh-2)+a(mesh-3)+&
548 : 0 : & a(mesh-4)+a(mesh-5))
549 : : asm(mesh-2)=0.2_dp*(a(mesh)+a(mesh-1)+a(mesh-2)+&
550 : 0 : & a(mesh-3)+a(mesh-4))
551 : 0 : asm(mesh-1)=0.25_dp*(a(mesh)+a(mesh-1)+a(mesh-2)+a(mesh-3))
552 : 0 : asm(mesh)=1.0_dp/3.0_dp*(a(mesh)+a(mesh-1)+a(mesh-2))
553 [ # # ]: 0 : 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 : & 0.05_dp*(a(i+5)+a(i-5))
559 : : end do
560 [ # # ]: 0 : do i=1,mesh
561 : 0 : a(i)=asm(i)
562 : : end do
563 : : end do
564 : :
565 : 0 : end subroutine paw_smooth
566 : : !!***
567 : :
568 : : !----------------------------------------------------------------------
569 : :
570 : : !!****f* m_paw_numeric/paw_sort_dp
571 : : !! NAME
572 : : !! paw_sort_dp
573 : : !!
574 : : !! FUNCTION
575 : : !! Sort real(dp) array list(n) into ascending numerical order using Heapsort
576 : : !! algorithm, while making corresponding rearrangement of the integer
577 : : !! array iperm. Consider that two real(dp) numbers
578 : : !! within tolerance tol are equal.
579 : : !!
580 : : !! INPUTS
581 : : !! n intent(in) dimension of the list
582 : : !! tol intent(in) numbers within tolerance are equal
583 : : !! list(n) intent(inout) list of real(dp) numbers to be sorted
584 : : !! iperm(n) intent(inout) iperm(i)=i (very important)
585 : : !!
586 : : !! OUTPUT
587 : : !! list(n) sorted list
588 : : !! iperm(n) index of permutation given the right ascending order
589 : : !!
590 : : !! SOURCE
591 : :
592 : 3 : 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 [ - + ]: 3 : if (n==1) return
613 : :
614 : : !Should not call with n<1
615 [ - + ]: 3 : if (n<1) then
616 : : write(msg,'(a,i12,2a)') &
617 : 0 : & 'paw_sort_dp has been called with array length n=',n, ch10, &
618 : 0 : & ' having a value less than 1. This is not allowed.'
619 : 0 : LIBPAW_ERROR(msg)
620 : : end if
621 : :
622 : : !Conduct the usual sort
623 : 3 : l=n/2+1 ; ir=n
624 : 6558 : do ! Infinite do-loop
625 [ + + ]: 6561 : if (l>1) then
626 : 2187 : l=l-1
627 : 2187 : ap=list(l)
628 : 2187 : iap=iperm(l)
629 : : else ! l<=1
630 : 4374 : ap=list(ir)
631 : 4374 : iap=iperm(ir)
632 : 4374 : list(ir)=list(1)
633 : 4374 : iperm(ir)=iperm(1)
634 : 4374 : ir=ir-1
635 [ + + ]: 4374 : if (ir==1) then
636 : 3 : list(1)=ap
637 : 3 : iperm(1)=iap
638 : : exit ! This is the end of this algorithm
639 : : end if
640 : : end if ! l>1
641 : 6558 : i=l
642 : 6558 : j=l+l
643 [ + + ]: 45699 : do while (j<=ir)
644 [ + + ]: 39141 : if (j<ir) then
645 [ + + + + : 39099 : if ( list(j)<list(j+1)-tol .or. &
+ + ]
646 : 39141 : & (list(j)<list(j+1)+tol.and.iperm(j)<iperm(j+1))) j=j+1
647 : : endif
648 [ + + + + : 45699 : if (ap<list(j)-tol.or.(ap<list(j)+tol.and.iap<iperm(j))) then
+ + ]
649 : 37401 : list(i)=list(j)
650 : 37401 : iperm(i)=iperm(j)
651 : 37401 : i=j
652 : 37401 : j=j+j
653 : : else
654 : 1740 : j=ir+1
655 : : end if
656 : : end do
657 : 6558 : list(i)=ap
658 : 6558 : iperm(i)=iap
659 : : end do ! End infinite do-loop
660 : :
661 : : end subroutine paw_sort_dp
662 : : !!***
663 : :
664 : : !----------------------------------------------------------------------
665 : :
666 : : !!****f* m_paw_numeric/paw_jbessel
667 : : !! NAME
668 : : !! paw_jbessel
669 : : !!
670 : : !! FUNCTION
671 : : !! Compute spherical Bessel function j_l(x) and derivative(s)
672 : : !!
673 : : !! INPUTS
674 : : !! ll=l-order of the Bessel function
675 : : !! order=1 if first derivative is requested
676 : : !! 2 if first and second derivatives are requested
677 : : !! xx=where to compute j_l
678 : : !!
679 : : !! OUTPUT
680 : : !! bes= Bessel function j_l at xx
681 : : !! besp= first derivative of j_l at xx (only if order>=1)
682 : : !! bespp= second derivative of j_l at xx (only if order=2)
683 : : !!
684 : : !! SOURCE
685 : :
686 : 427166668 : 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 [ - + ]: 427166668 : if (order>2) then
705 : 0 : msg='Wrong order in paw_jbessel!'
706 : 0 : LIBPAW_ERROR(msg)
707 : : end if
708 : :
709 [ + + ]: 427166668 : if (abs(xx)<prec) then
710 [ + + ]: 424672 : bes=zero;if (ll==0) bes=one
711 [ + + ]: 424672 : if (order>=1) then
712 [ + + ]: 21795 : besp=zero;if (ll==1) besp=third
713 : : end if
714 [ + + ]: 424672 : if (order==2) then
715 : 18 : bespp=zero
716 [ + + ]: 18 : if (ll==0) bespp=-third
717 [ + + ]: 18 : if (ll==2) bespp=2._dp/15._dp
718 : : end if
719 : 424672 : return
720 : : end if
721 : :
722 : 426741996 : xxinv=one/xx
723 : : if (order==0) then
724 : : factp=zero
725 : 426741996 : factpp=zero
726 : : jnp=zero
727 : 426741996 : jnpp=zero
728 : : end if
729 : :
730 [ + + ]: 426741996 : if (xx<one) then
731 : 155721139 : xx2=0.5_dp*xx*xx
732 : 155721139 : fact=one
733 [ + + ]: 751687890 : do il=1,ll
734 : 751687890 : fact=fact*xx/dble(2*il+1)
735 : : end do
736 : : jn=one;jr=one;ii=0
737 [ + + - + ]: 834497592 : do while(abs(jr)>=prec.and.ii<imax)
738 : 678776453 : ii=ii+1;jr=-jr*xx2/dble(ii*(2*(ll+ii)+1))
739 : 678776453 : jn=jn+jr
740 : : end do
741 : 155721139 : bes=jn*fact
742 [ - + ]: 155721139 : if (abs(jr)>prec) then
743 : 0 : msg='Bessel function did not converge!'
744 : 0 : LIBPAW_ERROR(msg)
745 : : end if
746 [ + + ]: 155721139 : if (order>=1) then
747 : 235421 : factp=fact*xx/dble(2*ll+3)
748 : 235421 : jnp=one;jr=one;ii=0
749 [ + + - + ]: 1869418 : do while(abs(jr)>=prec.AND.ii<imax)
750 : 1633997 : ii=ii+1;jr=-jr*xx2/dble(ii*(2*(ll+ii)+3))
751 : 1633997 : jnp=jnp+jr
752 : : end do
753 : 235421 : besp=-jnp*factp+jn*fact*xxinv*dble(ll)
754 [ - + ]: 235421 : if (abs(jr)>prec) then
755 : 0 : msg='1st der. of Bessel function did not converge!'
756 : 0 : LIBPAW_ERROR(msg)
757 : : end if
758 : : end if
759 [ + + ]: 155721139 : if (order==2) then
760 : 1321 : factpp=factp*xx/dble(2*ll+5)
761 : 1321 : jnpp=one;jr=one;ii=0
762 [ + + - + ]: 10984 : do while(abs(jr)>=prec.AND.ii<imax)
763 : 9663 : ii=ii+1;jr=-jr*xx2/dble(ii*(2*(ll+ii)+5))
764 : 9663 : jnpp=jnpp+jr
765 : : end do
766 : 1321 : besp1=-jnpp*factpp+jnp*factp*xxinv*dble(ll+1)
767 [ - + ]: 1321 : if (abs(jr)>prec) then
768 : 0 : msg='2nd der. of Bessel function did not converge !'
769 : 0 : LIBPAW_ERROR(msg)
770 : : end if
771 : : end if
772 : : else
773 : 271020857 : jn =sin(xx)*xxinv
774 : 271020857 : jnp=(-cos(xx)+jn)*xxinv
775 [ + + ]: 767060076 : do il=2,ll+1
776 : 496039219 : jr=-jn+dble(2*il-1)*jnp*xxinv
777 : 767060076 : jn=jnp;jnp=jr
778 : : end do
779 : 271020857 : bes=jn
780 [ + + ]: 271020857 : if (order>=1) besp =-jnp+jn *xxinv*dble(ll)
781 [ + + ]: 271020857 : if (order==2) besp1= jn -jnp*xxinv*dble(ll+2)
782 : : end if
783 : :
784 [ + + ]: 426741996 : if (order==2) bespp=-besp1+besp*ll*xxinv-bes*ll*xxinv*xxinv
785 : :
786 : : end subroutine paw_jbessel
787 : : !!***
788 : :
789 : : !----------------------------------------------------------------------
790 : :
791 : : !!****f* m_paw_numeric/paw_solvbes
792 : : !! NAME
793 : : !! paw_solvbes
794 : : !!
795 : : !! FUNCTION
796 : : !! Find nq first roots of instrinsic equation:
797 : : !! alpha.jl(Q) + beta.Q.djl/dr(Q) = 0
798 : : !!
799 : : !! INPUTS
800 : : !! alpha,beta= factors in intrinsic equation
801 : : !! ll= l quantum number
802 : : !! nq= number of roots to find
803 : : !!
804 : : !! OUTPUT
805 : : !! root(nq)= roots of instrinsic equation
806 : : !!
807 : : !! SOURCE
808 : :
809 : 430 : 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 : 430 : qq=dh;nroot=0
827 : :
828 [ + + ]: 1290 : do while (nroot<nq)
829 : 860 : call paw_jbessel(jbes,jbesp,dum,ll,1,qq)
830 : 860 : y1=alpha*jbes+beta*qq*jbesp
831 : 860 : qq=qq+dh
832 : 860 : call paw_jbessel(jbes,jbesp,dum,ll,1,qq)
833 : 860 : y2=alpha*jbes+beta*qq*jbesp
834 : :
835 [ + + ]: 35787 : do while (y1*y2>=zero)
836 : 34927 : qq=qq+dh
837 : 34927 : call paw_jbessel(jbes,jbesp,dum,ll,1,qq)
838 : 34927 : y2=alpha*jbes+beta*qq*jbesp
839 : : end do
840 : :
841 : 860 : hh=dh;qx=qq
842 [ + + ]: 38700 : do while (hh>tol)
843 : 37840 : hh=half*hh
844 [ + + ]: 37840 : if (y1*y2<zero) then
845 : 19983 : qx=qx-hh
846 : : else
847 : 17857 : qx=qx+hh
848 : : end if
849 : 37840 : call paw_jbessel(jbes,jbesp,dum,ll,1,qx)
850 : 37840 : y2=alpha*jbes+beta*qx*jbesp
851 : : end do
852 : 860 : nroot=nroot+1
853 : 860 : root(nroot)=qx
854 : :
855 : : end do
856 : :
857 : 430 : end subroutine paw_solvbes
858 : : !!***
859 : :
860 : : !----------------------------------------------------------------------
861 : :
862 : : !!****f* m_special_funcs/paw_jbessel_4spline
863 : : !! NAME
864 : : !! paw_jbessel_4spline
865 : : !!
866 : : !! FUNCTION
867 : : !! Compute spherical Bessel functions and derivatives.
868 : : !! A polynomial approximation is employed for q-->0.
869 : : !!
870 : : !! INPUTS
871 : : !! ll=l-order of the Bessel function
872 : : !! tol=tolerance below which a Polynomial approximation is employed
873 : : !! both for jl and its derivative (if required)
874 : : !! order=1 if only first derivative is requested
875 : : !! 2 if first and second derivatives are requested
876 : : !! xx=where to compute j_l
877 : : !!
878 : : !! OUTPUT
879 : : !! bes=Spherical Bessel function j_l at xx
880 : : !! besp= first derivative of j_l at xx (only if order>=1)
881 : : !!
882 : : !! TODO
883 : : !! Remove inline definitions, they are obsolete in F2003
884 : : !!
885 : : !! SOURCE
886 : :
887 : 9060029424 : 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 [ - + ]: 9060029424 : if (order>2) then
928 : 0 : msg='Wrong order in paw_jbessel_4spline'
929 : 0 : LIBPAW_ERROR(msg)
930 : : end if
931 : :
932 : 12736848456 : select case (ll)
933 : : case (0)
934 [ + + ]: 3676819032 : if (xx<TOL) then
935 : 265831267 : bes=1.0_dp-xx**2/6.0_dp*(1.0_dp-xx**2/20.0_dp)
936 [ + + ]: 265831267 : if (order>=1) besp=(-10.0_dp+xx*xx)*xx/30.0_dp
937 : : else
938 : 3410987765 : bes=sin(xx)/xx
939 [ + + ]: 3410987765 : if (order>=1) besp=-(sin(xx)-xx*cos(xx))/xx**2
940 : : end if
941 : :
942 : : case (1)
943 [ + + ]: 3630100392 : if (xx<TOL) then
944 : 263810582 : bes=(10.0_dp-xx*xx)*xx/30.0_dp
945 [ + + ]: 263810582 : if (order>=1) besp=(10.0_dp-3.0_dp*xx*xx)/30.0_dp
946 : : else
947 : 3366289810 : bes=(sin(xx)-xx*cos(xx))/xx**2
948 [ + + ]: 3366289810 : 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 [ + + ]: 1251390000 : if (xx<TOL) then
953 : 96800308 : bes=xx*xx/15.0_dp-xx**4/210.0_dp
954 [ + + ]: 96800308 : if (order>=1) besp=(1.0_dp-xx*xx/7.0_dp)*2.0_dp*xx/15.0_dp
955 : : else
956 : 1154589692 : bes=((3.0_dp-xx**2)*sin(xx)-3.0_dp*xx*cos(xx))/xx**3
957 [ + + ]: 1154589692 : 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 [ + + ]: 283083000 : if (xx<TOL) then
962 : 14251136 : bes=xx*xx*xx/105.0_dp-xx**5/1890.0_dp+xx**7/83160.0_dp
963 [ + + ]: 14251136 : if (order>=1) besp=(1.0_dp/35-xx*xx/378.0_dp+xx**4/11880.0_dp)*xx*xx
964 : : else
965 : 268831864 : 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 [ + + ]: 268831864 : 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 : 218637000 : call paw_jbessel(bes,besp,bespp,ll,order,xx)
971 : :
972 : : case default
973 : 0 : write(msg,'(a,i4)')' wrong value for ll = ',ll
974 [ + + + + : 9060029424 : LIBPAW_BUG(msg)
+ - ]
975 : : end select
976 : :
977 : 9060029424 : end subroutine paw_jbessel_4spline
978 : : !!***
979 : :
980 : : !----------------------------------------------------------------------
981 : :
982 : : !!****f* m_paw_numeric/paw_derfc
983 : : !! NAME
984 : : !! paw_derfc
985 : : !!
986 : : !! FUNCTION
987 : : !! Evaluates the complementary error function in real(dp).
988 : : !!
989 : : !! INPUTS
990 : : !! yy
991 : : !!
992 : : !! OUTPUT
993 : : !! derfc_yy=complementary error function of yy
994 : : !!
995 : : !! SOURCE
996 : :
997 : 264125101 : 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 : 264125101 : xx = yy
1042 : 264125101 : isw = 1
1043 : : !Here change the sign of xx, and keep track of it thanks to isw
1044 [ - + ]: 264125101 : if (xx<0.0e0_dp) then
1045 : 0 : isw = -1
1046 : 0 : xx = -xx
1047 : : end if
1048 : :
1049 : 264125101 : done=0
1050 : :
1051 : : !Residual value, if yy < -6.375e0_dp
1052 : 264125101 : res=2.0e0_dp
1053 : :
1054 : : !abs(yy) < .477, evaluate approximation for erfc
1055 [ + - ]: 264125101 : if (xx<0.477e0_dp) then
1056 : : ! xmin is a very small number
1057 [ - + ]: 264125101 : if (xx<xmin) then
1058 : 0 : res = xx*pp(3)/qq(3)
1059 : : else
1060 : 264125101 : xsq = xx*xx
1061 : 264125101 : xnum = pp(4)*xsq+pp(5)
1062 : 264125101 : xden = xsq+qq(4)
1063 [ + + ]: 1056500404 : do ii = 1,3
1064 : 792375303 : xnum = xnum*xsq+pp(ii)
1065 : 1056500404 : xden = xden*xsq+qq(ii)
1066 : : end do
1067 : 264125101 : res = xx*xnum/xden
1068 : : end if
1069 [ - + ]: 264125101 : if (isw==-1) res = -res
1070 : 264125101 : res = 1.0e0_dp-res
1071 : 264125101 : done=1
1072 : : end if
1073 : :
1074 : : !.477 < abs(yy) < 4.0 , evaluate approximation for erfc
1075 [ + - - + ]: 264125101 : if (xx<=4.0e0_dp .and. done==0 ) then
1076 : 0 : xsq = xx*xx
1077 : 0 : xnum = p1(8)*xx+p1(9)
1078 : 0 : xden = xx+q1(8)
1079 [ # # ]: 0 : do ii=1,7
1080 : 0 : xnum = xnum*xx+p1(ii)
1081 : 0 : xden = xden*xx+q1(ii)
1082 : : end do
1083 : 0 : res = xnum/xden
1084 : 0 : res = res* exp(-xsq)
1085 [ # # ]: 0 : if (isw.eq.-1) res = 2.0e0_dp-res
1086 : : done=1
1087 : : end if
1088 : :
1089 : : !y > 13.3e0_dp
1090 [ + - - + : 264125101 : if (isw > 0 .and. xx > xbig .and. done==0 ) then
- - ]
1091 : 0 : res = 0.0e0_dp
1092 : 0 : 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 [ - + - - : 264125101 : if ( ( isw > 0 .or. xx < xlarge ) .and. done==0 ) then
- + ]
1098 : 0 : xsq = xx*xx
1099 : 0 : xi = 1.0e0_dp/xsq
1100 : 0 : xnum= p2(5)*xi+p2(6)
1101 : 0 : xden = xi+q2(5)
1102 [ # # ]: 0 : do ii = 1,4
1103 : 0 : xnum = xnum*xi+p2(ii)
1104 : 0 : xden = xden*xi+q2(ii)
1105 : : end do
1106 : 0 : res = (sqrpi+xi*xnum/xden)/xx
1107 : 0 : res = res* exp(-xsq)
1108 [ # # ]: 0 : if (isw.eq.-1) res = 2.0e0_dp-res
1109 : : end if
1110 : :
1111 : : !All cases have been investigated
1112 : 264125101 : derfc_yy = res
1113 : :
1114 : 264125101 : end function paw_derfc
1115 : : !!***
1116 : :
1117 : : !----------------------------------------------------------------------
1118 : :
1119 : : end module m_paw_numeric
1120 : : !!***
|
