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 6980 : 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 6980 : real(dp),allocatable :: tmp(:)
145 :
146 : ! *************************************************************************
147 :
148 : !Check
149 6980 : 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 20940 : LIBPAW_ALLOCATE(tmp,(n))
158 :
159 21076074 : do i=1,n-1
160 21076074 : 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 6980 : ibcbeg=1;if(ybcbeg>1.0d+30)ibcbeg=0
171 6980 : ibcend=1;if(ybcend>1.0d+30)ibcend=0
172 :
173 : !Set the first and last equations
174 6980 : if (ibcbeg==0) then
175 0 : ypp(1) = 0._dp
176 0 : tmp(1) = 0._dp
177 : else if ( ibcbeg == 1 ) then
178 6980 : ypp(1) = -0.5_dp
179 6980 : tmp(1) = (3._dp/(t(2)-t(1)))*((y(2)-y(1))/(t(2)-t(1))-ybcbeg)
180 : end if
181 6980 : if (ibcend==0) then
182 0 : ypp(n) = 0._dp
183 0 : tmp(n) = 0._dp
184 : else if ( ibcend == 1 ) then
185 6980 : ypp(n) = 0.5_dp
186 6980 : 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 21069094 : do i=2,n-1
191 21062114 : ratio=(t(i)-t(i-1))/(t(i+1)-t(i-1))
192 21062114 : pinv = 1.0_dp/(ratio*ypp(i-1) + 2.0_dp)
193 21062114 : 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 21062114 : & /(t(i)-t(i-1)))/(t(i+1)-t(i-1))-ratio*tmp(i-1))*pinv
196 21069094 : if (abs(tmp(i))<1.d5*tiny(0._dp)) tmp(i)=0._dp
197 : end do
198 :
199 : !Solve the equations
200 6980 : ypp(n) = (tmp(n)-ypp(n)*tmp(n-1))/(ypp(n)*ypp(n-1)+1.0_dp)
201 21076074 : do k=n-1,1,-1
202 21076074 : ypp(k)=ypp(k)*ypp(k+1)+tmp(k)
203 : end do
204 :
205 6980 : LIBPAW_DEALLOCATE(tmp)
206 :
207 6980 : 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 72 : 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 72 : my_err=0
256 72 : left=1
257 120302 : do i=1,nfit
258 120230 : yfit(i)=0._dp ! Initialize for the unlikely event that rmax exceed r(mesh)
259 253878 : do k=left+1, nspline
260 253878 : if(xspline(k) >= xfit(i)) then
261 120230 : 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 120230 : delarg= xspline(right) - xspline(left)
273 120230 : invdelarg= 1.0_dp/delarg
274 120230 : aa= (xspline(right)-xfit(i))*invdelarg
275 120230 : bb= (xfit(i)-xspline(left))*invdelarg
276 : yfit(i) = aa*yspline(left)+bb*yspline(right) &
277 : & +( (aa*aa*aa-aa)*ysplin2(left) + &
278 120230 : & (bb*bb*bb-bb)*ysplin2(right) ) *delarg*delarg/6.0_dp
279 120230 : exit
280 : end if
281 : end do ! k
282 120302 : if (k==nspline+1) my_err=my_err+1 ! xfit not found
283 : end do ! i
284 72 : if (present(ierr)) ierr=my_err
285 :
286 72 : 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 413423668 : 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 413423668 : if (order>2) then
705 0 : msg='Wrong order in paw_jbessel!'
706 0 : LIBPAW_ERROR(msg)
707 : end if
708 :
709 413423668 : if (abs(xx)<prec) then
710 415672 : bes=zero;if (ll==0) bes=one
711 415672 : if (order>=1) then
712 21792 : besp=zero;if (ll==1) besp=third
713 : end if
714 415672 : 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 415672 : return
720 : end if
721 :
722 413007996 : xxinv=one/xx
723 : if (order==0) then
724 : factp=zero
725 413007996 : factpp=zero
726 : jnp=zero
727 413007996 : jnpp=zero
728 : end if
729 :
730 413007996 : if (xx<one) then
731 146555749 : xx2=0.5_dp*xx*xx
732 146555749 : fact=one
733 705860940 : do il=1,ll
734 705860940 : fact=fact*xx/dble(2*il+1)
735 : end do
736 : jn=one;jr=one;ii=0
737 787025829 : do while(abs(jr)>=prec.and.ii<imax)
738 640470080 : ii=ii+1;jr=-jr*xx2/dble(ii*(2*(ll+ii)+1))
739 640470080 : jn=jn+jr
740 : end do
741 146555749 : bes=jn*fact
742 146555749 : if (abs(jr)>prec) then
743 0 : msg='Bessel function did not converge!'
744 0 : LIBPAW_ERROR(msg)
745 : end if
746 146555749 : if (order>=1) then
747 232844 : factp=fact*xx/dble(2*ll+3)
748 232844 : jnp=one;jr=one;ii=0
749 1855522 : do while(abs(jr)>=prec.AND.ii<imax)
750 1622678 : ii=ii+1;jr=-jr*xx2/dble(ii*(2*(ll+ii)+3))
751 1622678 : jnp=jnp+jr
752 : end do
753 232844 : besp=-jnp*factp+jn*fact*xxinv*dble(ll)
754 232844 : 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 146555749 : 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 266452247 : jn =sin(xx)*xxinv
774 266452247 : jnp=(-cos(xx)+jn)*xxinv
775 744217026 : do il=2,ll+1
776 477764779 : jr=-jn+dble(2*il-1)*jnp*xxinv
777 744217026 : jn=jnp;jnp=jr
778 : end do
779 266452247 : bes=jn
780 266452247 : if (order>=1) besp =-jnp+jn *xxinv*dble(ll)
781 266452247 : if (order==2) besp1= jn -jnp*xxinv*dble(ll+2)
782 : end if
783 :
784 413007996 : 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 8730043512 : 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 8730043512 : if (order>2) then
928 0 : msg='Wrong order in paw_jbessel_4spline'
929 0 : LIBPAW_ERROR(msg)
930 : end if
931 :
932 12279645528 : select case (ll)
933 : case (0)
934 3549602016 : if (xx<TOL) then
935 258097396 : bes=1.0_dp-xx**2/6.0_dp*(1.0_dp-xx**2/20.0_dp)
936 258097396 : if (order>=1) besp=(-10.0_dp+xx*xx)*xx/30.0_dp
937 : else
938 3291504620 : bes=sin(xx)/xx
939 3291504620 : if (order>=1) besp=-(sin(xx)-xx*cos(xx))/xx**2
940 : end if
941 :
942 : case (1)
943 3502514496 : if (xx<TOL) then
944 256076711 : bes=(10.0_dp-xx*xx)*xx/30.0_dp
945 256076711 : if (order>=1) besp=(10.0_dp-3.0_dp*xx*xx)/30.0_dp
946 : else
947 3246437785 : bes=(sin(xx)-xx*cos(xx))/xx**2
948 3246437785 : 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 1206729000 : if (xx<TOL) then
953 94628755 : bes=xx*xx/15.0_dp-xx**4/210.0_dp
954 94628755 : if (order>=1) besp=(1.0_dp-xx*xx/7.0_dp)*2.0_dp*xx/15.0_dp
955 : else
956 1112100245 : bes=((3.0_dp-xx**2)*sin(xx)-3.0_dp*xx*cos(xx))/xx**3
957 1112100245 : 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 266304000 : if (xx<TOL) then
962 13265407 : bes=xx*xx*xx/105.0_dp-xx**5/1890.0_dp+xx**7/83160.0_dp
963 13265407 : if (order>=1) besp=(1.0_dp/35-xx*xx/378.0_dp+xx**4/11880.0_dp)*xx*xx
964 : else
965 253038593 : 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 253038593 : 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 204894000 : 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 8730043512 : LIBPAW_BUG(msg)
975 : end select
976 :
977 8730043512 : 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 : !!***
|
