Line data Source code
1 : !!****m* ABINIT/m_special_funcs
2 : !! NAME
3 : !! m_special_funcs
4 : !!
5 : !! FUNCTION
6 : !! This module contains routines and functions used to
7 : !! evaluate special functions frequently needed in Abinit.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2008-2026 ABINIT group (MG, MT, FB, XG, MVer, FJ, NH, GZ, DRH)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_special_funcs
24 :
25 : use defs_basis
26 : use m_abicore
27 : use m_errors
28 : use m_splines
29 :
30 : use m_fstrings, only : sjoin, ftoa
31 : use m_numeric_tools, only : arth, simpson
32 :
33 : implicit none
34 :
35 : private
36 :
37 : public :: clp ! x-1, if x>1/2, x+1, if x<-1/2
38 : public :: factorial ! Calculates N! returning a real.
39 : public :: permutations ! Returns N!/(N-k) if N>=0 and N-k>0 else 0.
40 : public :: binomcoeff ! Binomial coefficient n!/(n-k)!
41 : public :: laguerre ! Laguerre Polynomial(x,n,a).
42 : public :: RadFnH ! Atomic radial function(r,n,l,Z).
43 : public :: iradfnh ! Norm of atomic radial function(a,b,n,l,Z).
44 : public :: gaussian ! Normalized Gaussian distribution.
45 : public :: lorentzian ! Approximate Dirac Delta with lorentzian
46 : public :: abi_derf ! Evaluates the error function in real(dp).
47 : public :: abi_derfc ! Evaluates the complementary error function in real(dp).
48 : public :: gamma_function ! Computes the gamma function
49 : public :: besjm ! Spherical bessel function of order nn. Handles nn=0,1,2,3,4, or 5 only.
50 : public :: sbf8 ! Computes set of spherical bessel functions using accurate algorithm
51 : public :: k_fermi ! Fermi wave vector corresponding to the local value of the real space density rhor.
52 : public :: k_thfermi ! Thomas-Fermi wave vector corresponding to the local value of the real space density rhor
53 : public :: levi_civita_3 ! Return Levi-Civita tensor of rank 3
54 : public :: fermi_dirac ! Fermi Dirac distribution
55 : public :: bose_einstein ! Bose Einstein distribution
56 : public :: dip12 ! Complete Fermi integral of order 1/2
57 : public :: dip32 ! Complete Fermi integral of order 3/2
58 : public :: djp12 ! Incomplete Fermi integral of order 1/2
59 : public :: djp32 ! Incomplete Fermi integral of order 3/2
60 : public :: tildeAx ! tilde Ax Pade fit and first and second derivatives
61 : public :: tildeBx ! tilde Bx Pade fit and first and second derivatives
62 : public :: tildeBc ! tilde Bc Pade fit and first and second derivatives
63 : !!***
64 :
65 : !!****t* m_special_funcs/jlspline_t
66 : !! NAME
67 : !! jlspline_t
68 : !!
69 : !! FUNCTION
70 : !! Object used to interpolate Bessel functions
71 : !!
72 : !! SOURCE
73 :
74 : type,public :: jlspline_t
75 :
76 : integer :: nx
77 : ! number of points on linear mesh used in spline.
78 :
79 : integer :: mlang
80 : ! mlang= max angular momentum + 1
81 :
82 : real(dp) :: delta
83 : ! Step of linear mesh.
84 :
85 : real(dp) :: maxarg
86 : ! max arg value.
87 :
88 : real(dp),allocatable :: xx(:)
89 : ! xx(nx)
90 : ! coordinates of points belonging to the grid
91 :
92 : real(dp),allocatable :: bess_spl(:,:)
93 : ! bess_spl(nx,mlang)
94 : ! bessel functions computed on the linear mesh
95 :
96 : real(dp),allocatable :: bess_spl_der(:,:)
97 : ! bess_spl_der(nx,mlang)
98 : ! the second derivatives of the cubic spline.
99 :
100 : contains
101 :
102 : procedure :: init => jlspline_init ! Create new object.
103 : procedure :: free => jlspline_free ! Free memory.
104 : procedure :: eval => jlspline_integral ! Compute integral.
105 :
106 : end type jlspline_t
107 : !!***
108 :
109 : !!****t* m_special_funcs/gspline_t
110 : !! NAME
111 : !! gspline_t
112 : !!
113 : !! FUNCTION
114 : !! Object used to interpolate the gaussian approximant and its primitive with cubic spline.
115 : !! Particularly useful if we are computing DOSes with many k-points/bands
116 : !! because one can significantly decrease the number of calls to exponential functions.
117 : !!
118 : !! SOURCE
119 :
120 : type,public :: gspline_t
121 :
122 : integer :: nspline
123 : ! Number of points used in spline table.
124 :
125 : real(dp) :: sigma
126 : ! Broadening parameter.
127 :
128 : real(dp) :: xmin, xmax
129 : ! Min and max x in spline mesh. Only positive xs are stored in memory
130 : ! The values at -x are reconstructed by symmetry.
131 : ! xmin is usually zero, xmax is the point where the gaussian == tol16.
132 : ! g(x) is set to zero if x > xmin.
133 :
134 : real(dp) :: step, stepm1, step2div6
135 : ! Step of the linear mesh used in spline and associated coefficients.
136 :
137 : real(dp),allocatable :: xvals(:)
138 : ! xvals(nspline)
139 : ! The xvalues used in the spline
140 :
141 : real(dp),allocatable :: svals(:,:)
142 : ! svals(nspline,4)
143 : ! Internal tables with spline data.
144 :
145 : contains
146 :
147 : procedure :: init => gspline_init ! Creation method.
148 : procedure :: eval => gspline_eval ! Evaluate interpolant
149 : procedure :: free => gspline_free ! Free memory.
150 : end type gspline_t
151 : !!***
152 :
153 : CONTAINS !===========================================================
154 : !!***
155 :
156 : !!****f* m_special_funcs/clp
157 : !! NAME
158 : !! clp
159 : !!
160 : !! FUNCTION
161 : !! clp(x)= x-1, if x>1/2
162 : !! x+1, if x<-1/2
163 : !!
164 : !! INPUTS
165 : !! x= input variable
166 : !!
167 : !! OUTPUT
168 : !! clp= resulting function
169 : !!
170 : !! SOURCE
171 :
172 0 : pure function clp(x)
173 :
174 : !Arguments ------------------------------------
175 : !scalars
176 : real(dp) :: clp
177 : real(dp),intent(in) :: x
178 : ! **********************************************************************
179 :
180 0 : if(x > half) then
181 0 : clp=x-one
182 0 : elseif(x < -half) then
183 0 : clp=x+one
184 : else
185 : clp=x
186 : end if
187 :
188 0 : end function clp
189 : !!***
190 :
191 : !!****f* m_special_funcs/factorial
192 : !! NAME
193 : !! factorial
194 : !!
195 : !! FUNCTION
196 : !! Calculates N!. Returns a (dp) real.
197 : !!
198 : !! INPUTS
199 : !! nn=number to use
200 : !!
201 : !! OUTPUT
202 : !! factorial= n! (real)
203 : !!
204 : !! SOURCE
205 :
206 47745 : elemental function factorial(nn)
207 :
208 : !Arguments ---------------------------------------------
209 : !scalars
210 : integer,intent(in) :: nn
211 : real(dp) :: factorial
212 :
213 : !Local variables ---------------------------------------
214 : !scalars
215 : integer :: ii
216 : real(dp) :: ff
217 : ! *********************************************************************
218 :
219 47745 : ff=one
220 115523 : do ii=2,nn
221 284483 : ff=ff*ii
222 : end do
223 :
224 47745 : factorial=ff
225 :
226 1665 : end function factorial
227 : !!***
228 :
229 : !!****f* m_special_funcs/permutations
230 : !! NAME
231 : !! permutations
232 : !!
233 : !! FUNCTION
234 : !! Returns N!/(N-k)! if N>=0 and N-k>0
235 : !! otherwise 0 is returned
236 : !! Output is real
237 : !!
238 : !! INPUTS
239 : !! kk=number k to use
240 : !! nn=number N to use
241 : !!
242 : !! OUTPUT
243 : !! permutations= n!/(n-k)! (real)
244 : !!
245 : !! SOURCE
246 :
247 1665 : pure function permutations(nn,kk)
248 :
249 : !Arguments ---------------------------------------------
250 : !scalars
251 : integer,intent(in) :: kk,nn
252 : real(dp) :: permutations
253 :
254 : !Local variables ---------------------------------------
255 : !scalars
256 : integer :: ii
257 : real(dp) :: pp
258 : ! *********************************************************************
259 :
260 1665 : if ((nn>=0).and.((nn-kk)>=0)) then
261 1665 : pp=one
262 9575 : do ii=nn-kk+1,nn
263 9575 : pp=pp*ii
264 : end do
265 : else
266 : pp=zero
267 : end if
268 :
269 1665 : permutations=pp
270 :
271 1665 : end function permutations
272 : !!***
273 :
274 : !----------------------------------------------------------------------
275 :
276 : !!****f* m_special_funcs/binomcoeff
277 : !! NAME
278 : !! factorial
279 : !!
280 : !! FUNCTION
281 : !! Calculates n!/( k!* (n-k)!). Returns a real (dp)
282 : !!
283 : !! INPUTS
284 : !! nn=number to use
285 : !!
286 : !! OUTPUT
287 : !! binomcoeff= n!/( k!* (n-k)!) (real dp)
288 : !!
289 : !! SOURCE
290 :
291 15360 : elemental function binomcoeff(n,k)
292 :
293 : !Arguments ---------------------------------------------
294 : !scalars
295 : integer,intent(in) :: n,k
296 : real(dp) :: binomcoeff
297 : ! *********************************************************************
298 :
299 107520 : binomcoeff=factorial(n)/(factorial(k)*factorial(n-k))
300 :
301 15360 : end function binomcoeff
302 : !!***
303 :
304 : !----------------------------------------------------------------------
305 :
306 : !!****f* m_special_funcs/laguerre
307 : !! NAME
308 : !! laguerre
309 : !!
310 : !! FUNCTION
311 : !! Laguerre(x,n,a). Returns a (dp) real.
312 : !!
313 : !! INPUTS
314 : !! x position
315 : !! n order of laguerre polynomial
316 : !! a
317 : !!
318 : !! OUTPUT
319 : !! Laguerre(x,n,a) (dp)
320 : !!
321 : !! SOURCE
322 :
323 15360 : function laguerre(x,n,a)
324 :
325 : !Arguments ---------------------------------------------
326 : !scalars
327 : integer,intent(in),optional :: n,a
328 : real(dp) :: laguerre
329 : real(dp),intent(in) :: x
330 :
331 : !Local variables ---------------------------------------
332 : !scalars
333 : integer :: ii, nn, aa
334 :
335 : !arrays
336 15360 : real(dp),allocatable :: ff(:)
337 : ! *********************************************************************
338 :
339 15360 : if (present(n)) then
340 15360 : nn=n
341 : else
342 : nn=1
343 : end if
344 :
345 15360 : if (present(a)) then
346 15360 : aa=a
347 : else
348 : aa=0
349 : end if
350 46080 : ABI_MALLOC(ff,(nn+1))
351 30720 : ff=0.0_dp
352 76800 : ff=(/ (binomcoeff(nn+aa,nn-ii)*((-1.0_dp)*x)**ii/factorial(ii) ,ii=0,nn) /)
353 30720 : laguerre=sum(ff)
354 :
355 15360 : ABI_FREE(ff)
356 :
357 15360 : end function laguerre
358 : !!***
359 :
360 : !----------------------------------------------------------------------
361 :
362 : !!****f* m_special_funcs/RadFnH
363 : !! NAME
364 : !! RadFnH
365 : !!
366 : !! FUNCTION
367 : !! RadFnH(r,n,l,Z) radial function of atomic wavefunction with nuclear charge Z.
368 : !! for quantum number n, and l.
369 : !! Default: Fe 3d function. Returns a (dp) real.
370 : !!
371 : !! INPUTS
372 : !! r radius
373 : !! n principal quantum number
374 : !! l quantum number
375 : !!
376 : !! OUTPUT
377 : !! RadFnH(r,n,l,Z) (dp)
378 : !!
379 : !! SOURCE
380 :
381 15360 : function RadFnH(r,n,l,Z)
382 :
383 : !Arguments ---------------------------------------------
384 : !scalars
385 : integer,intent(in),optional :: n,l
386 : real(dp) :: RadFnH
387 : real(dp),intent(in) :: r
388 : real(dp),intent(in),optional :: Z
389 :
390 : !Local variables ---------------------------------------
391 : !scalars
392 : integer :: nn,ll
393 : real(dp) :: ff,rr,ZZ
394 : ! *********************************************************************
395 :
396 15360 : if (present(n)) then
397 15360 : nn=n
398 : else
399 : nn=3
400 : end if
401 :
402 15360 : if (present(l)) then
403 15360 : ll=l
404 : else
405 : ll=2
406 : end if
407 :
408 15360 : if (present(Z)) then
409 15360 : ZZ=Z
410 : else
411 : ZZ=28.0_dp
412 : end if
413 :
414 15360 : rr=ZZ*r/nn
415 15360 : ff=exp(log(ZZ*1.0_dp)*(3.0_dp/2.0_dp))*2/nn**2
416 92160 : ff=ff*sqrt(factorial(nn-ll-1)/factorial(nn+ll))*(2*rr)**ll
417 15360 : RadFnH=ff*exp(-1*rr)*laguerre(2*rr,nn-ll-1,2*ll+1)
418 :
419 15360 : end function RadFnH
420 : !!***
421 :
422 : !----------------------------------------------------------------------
423 :
424 : !!****f* m_special_funcs/IRadFnH
425 : !! NAME
426 : !! IRadFnH
427 : !!
428 : !! FUNCTION
429 : !! IRadFnH(a,b,n,l,Z): Integral of radial function of atomic wavefunction between a and b.
430 : !! recursive programming using simpson's rule
431 : !! iteration depth of m=8 corresponds to relative error of 10^(-12).
432 : !!
433 : !! INPUTS
434 : !! a lower limit for integration
435 : !! b upper limit for integration
436 : !! n principal quantum number
437 : !! l quantum number
438 : ! Z nuclear charge
439 : !!
440 : !! OUTPUT
441 : !! IRadFnH(a,b,n,l,Z) (dp)
442 : !!
443 : !! SOURCE
444 :
445 10220 : recursive function IRadFnH(a,b,n,l,Z,m) result(x)
446 :
447 : !Arguments ---------------------------------------------
448 : !scalars
449 : integer,intent(in),optional :: n,l,m
450 : real(dp),intent(in):: a
451 : real(dp),intent(in),optional :: b,Z
452 :
453 : !Local variables ---------------------------------------
454 : !scalars
455 : integer :: nn,ll,mm
456 : real(dp) :: h,bb,ZZ,x
457 : ! *********************************************************************
458 :
459 10220 : if (present(n)) then
460 10220 : nn=n
461 : else
462 0 : nn=3
463 : end if
464 :
465 10220 : if (present(l)) then
466 10220 : ll=l
467 : else
468 0 : ll=2
469 : end if
470 :
471 10220 : if (present(Z)) then
472 10220 : ZZ=Z
473 : else
474 0 : ZZ=28
475 : end if
476 :
477 10220 : if (present(b)) then
478 10220 : bb=b
479 : else
480 0 : bb=100.0_dp
481 : end if
482 :
483 10220 : if (present(m)) then
484 10200 : mm=m
485 : else
486 : mm=0
487 : end if
488 :
489 10220 : h=(bb-a)/2.0_dp
490 10220 : if (mm<8) then
491 : !h=2*h/exp(1.0_dp)
492 5100 : x=IRadFnH(a,a+h,nn,ll,ZZ,mm+1)+IRadFnH(a+h,bb,nn,ll,ZZ,mm+1)
493 : else
494 5120 : x=RadFnH(a,nn,ll,ZZ)**2*a**2+4.0_dp*RadFnH(a+h,nn,ll,ZZ)**2*(a+h)**2
495 5120 : x=h/3.0_dp*(x+RadFnH(bb,nn,ll,ZZ)**2*bb**2)
496 : end if
497 :
498 10220 : end function IRadFnH
499 : !!***
500 :
501 : !----------------------------------------------------------------------
502 :
503 : !!****f* m_special_funcs/gaussian
504 : !! NAME
505 : !! gaussian
506 : !!
507 : !! FUNCTION
508 : !! Return the values of the normalized Gaussian distribution:
509 : !!
510 : !! Gauss(arg,sigma) = 1/(sigma SQRT(2*pi)) e^{-arg**2/(2*sigma**2)}
511 : !!
512 : !! INPUTS
513 : !! arg=Argument of the Gaussian.
514 : !! sigma=Standard deviation
515 : !!
516 : !! SOURCE
517 :
518 16907114 : elemental function gaussian(arg, sigma)
519 :
520 : !Arguments ---------------------------------------------
521 : !scalars
522 : real(dp),intent(in) :: arg,sigma
523 : real(dp) :: gaussian
524 :
525 : !Local variables ---------------------------------------
526 : real(dp) :: xx
527 : ! *********************************************************************
528 :
529 16907114 : xx = arg / (sqrt2 * sigma)
530 16907114 : gaussian = exp(-xx*xx) / (sigma * sqrt(two_pi))
531 :
532 16907114 : end function gaussian
533 : !!***
534 :
535 : !----------------------------------------------------------------------
536 :
537 : !!****f* m_special_funcs/lorentzian
538 : !! NAME
539 : !! lorentzian
540 : !!
541 : !! FUNCTION
542 : !! Lorentzian function.
543 : !!
544 : !! INPUTS
545 : !! arg=Argument of the lorentzian.
546 : !! sigma=Broadening factor
547 : !!
548 : !! SOURCE
549 :
550 0 : elemental function lorentzian(arg, sigma)
551 :
552 : !Arguments ---------------------------------------------
553 : real(dp),intent(in) :: arg, sigma
554 : real(dp) :: lorentzian
555 : ! *********************************************************************
556 :
557 0 : lorentzian = piinv * sigma / (arg ** 2 + sigma ** 2)
558 :
559 0 : end function lorentzian
560 : !!***
561 :
562 : !----------------------------------------------------------------------
563 :
564 : !!****f* m_special_funcs/abi_derf
565 : !! NAME
566 : !! abi_derf
567 : !!
568 : !! FUNCTION
569 : !! Evaluates the error function in real(dp).
570 : !! Same implementation as imsl.
571 : !! Simple mod of derfc.F90
572 : !!
573 : !! INPUTS
574 : !! yy
575 : !!
576 : !! OUTPUT
577 : !! derf_yy= error function of yy
578 : !!
579 : !! SOURCE
580 :
581 3278 : elemental function abi_derf(yy) result(derf_yy)
582 :
583 : !Arguments ------------------------------------
584 : !scalars
585 : real(dp),intent(in) :: yy
586 : real(dp) :: derf_yy
587 :
588 : !Local variables-------------------------------
589 : integer :: done,ii,isw
590 : ! coefficients for 0.0 <= yy < .477
591 : real(dp), parameter :: &
592 : & pp(5)=(/ 113.8641541510502e0_dp, 377.4852376853020e0_dp, &
593 : & 3209.377589138469e0_dp, .1857777061846032e0_dp, &
594 : & 3.161123743870566e0_dp /)
595 : real(dp), parameter :: &
596 : & qq(4)=(/ 244.0246379344442e0_dp, 1282.616526077372e0_dp, &
597 : & 2844.236833439171e0_dp, 23.60129095234412e0_dp/)
598 : ! coefficients for .477 <= yy <= 4.0
599 : real(dp), parameter :: &
600 : & p1(9)=(/ 8.883149794388376e0_dp, 66.11919063714163e0_dp, &
601 : & 298.6351381974001e0_dp, 881.9522212417691e0_dp, &
602 : & 1712.047612634071e0_dp, 2051.078377826071e0_dp, &
603 : & 1230.339354797997e0_dp, 2.153115354744038e-8_dp, &
604 : & .5641884969886701e0_dp /)
605 : real(dp), parameter :: &
606 : & q1(8)=(/ 117.6939508913125e0_dp, 537.1811018620099e0_dp, &
607 : & 1621.389574566690e0_dp, 3290.799235733460e0_dp, &
608 : & 4362.619090143247e0_dp, 3439.367674143722e0_dp, &
609 : & 1230.339354803749e0_dp, 15.74492611070983e0_dp/)
610 : ! coefficients for 4.0 < y,
611 : real(dp), parameter :: &
612 : & p2(6)=(/ -3.603448999498044e-01_dp, -1.257817261112292e-01_dp, &
613 : & -1.608378514874228e-02_dp, -6.587491615298378e-04_dp, &
614 : & -1.631538713730210e-02_dp, -3.053266349612323e-01_dp/)
615 : real(dp), parameter :: &
616 : & q2(5)=(/ 1.872952849923460e0_dp , 5.279051029514284e-01_dp, &
617 : & 6.051834131244132e-02_dp , 2.335204976268692e-03_dp, &
618 : & 2.568520192289822e0_dp /)
619 : real(dp), parameter :: &
620 : & sqrpi=.5641895835477563e0_dp, xbig=13.3e0_dp, xlarge=6.375e0_dp, xmin=1.0e-10_dp
621 : real(dp) :: res,xden,xi,xnum,xsq,xx
622 : ! ******************************************************************
623 :
624 3278 : xx = yy
625 3278 : isw = 1
626 : !Here change the sign of xx, and keep track of it thanks to isw
627 3278 : if (xx<0.0e0_dp) then
628 0 : isw = -1
629 0 : xx = -xx
630 : end if
631 :
632 3278 : done=0
633 :
634 : !Residual value, if yy < -6.375e0_dp
635 3278 : res=-1.0e0_dp
636 :
637 : !abs(yy) < .477, evaluate approximation for erfc
638 3278 : if (xx<0.477e0_dp) then
639 : ! xmin is a very small number
640 108 : if (xx<xmin) then
641 0 : res = xx*pp(3)/qq(3)
642 : else
643 108 : xsq = xx*xx
644 108 : xnum = pp(4)*xsq+pp(5)
645 108 : xden = xsq+qq(4)
646 432 : do ii = 1,3
647 324 : xnum = xnum*xsq+pp(ii)
648 432 : xden = xden*xsq+qq(ii)
649 : end do
650 108 : res = xx*xnum/xden
651 : end if
652 108 : if (isw==-1) res = -res
653 : done=1
654 : end if
655 :
656 : !.477 < abs(yy) < 4.0 , evaluate approximation for erfc
657 3278 : if (xx<=4.0e0_dp .and. done==0 ) then
658 1280 : xsq = xx*xx
659 1280 : xnum = p1(8)*xx+p1(9)
660 1280 : xden = xx+q1(8)
661 10240 : do ii=1,7
662 8960 : xnum = xnum*xx+p1(ii)
663 10240 : xden = xden*xx+q1(ii)
664 : end do
665 1280 : res = xnum/xden
666 1280 : res = res* exp(-xsq)
667 1280 : if (isw.eq.-1) then
668 0 : res = res-1.0e0_dp
669 : else
670 1280 : res=1.0e0_dp-res
671 : end if
672 1280 : done=1
673 : end if
674 :
675 : !y > 13.3e0_dp
676 3278 : if (isw > 0 .and. xx > xbig .and. done==0 ) then
677 0 : res = 1.0e0_dp
678 0 : done=1
679 : end if
680 :
681 : !4.0 < yy < 13.3e0_dp .or. -6.375e0_dp < yy < -4.0
682 : !evaluate minimax approximation for erfc
683 3278 : if ( ( isw > 0 .or. xx < xlarge ) .and. done==0 ) then
684 1890 : xsq = xx*xx
685 1890 : xi = 1.0e0_dp/xsq
686 1890 : xnum= p2(5)*xi+p2(6)
687 1890 : xden = xi+q2(5)
688 9450 : do ii = 1,4
689 7560 : xnum = xnum*xi+p2(ii)
690 9450 : xden = xden*xi+q2(ii)
691 : end do
692 1890 : res = (sqrpi+xi*xnum/xden)/xx
693 1890 : res = res* exp(-xsq)
694 1890 : if (isw.eq.-1) then
695 0 : res = res-1.0e0_dp
696 : else
697 1890 : res=1.0e0_dp-res
698 : end if
699 : end if
700 :
701 : !All cases have been investigated
702 3278 : derf_yy = res
703 :
704 3278 : end function abi_derf
705 : !!***
706 :
707 : !----------------------------------------------------------------------
708 :
709 : !!****f* m_special_funcs/abi_derfc
710 : !! NAME
711 : !! abi_derfc
712 : !!
713 : !! FUNCTION
714 : !! Evaluates the complementary error function in real(dp).
715 : !! Same implementation as imsl.
716 : !!
717 : !! INPUTS
718 : !! yy
719 : !!
720 : !! OUTPUT
721 : !! derfc_yy=complementary error function of yy
722 : !!
723 : !! SOURCE
724 :
725 1323947848 : elemental function abi_derfc(yy) result(derfc_yy)
726 :
727 : !Arguments ------------------------------------
728 : !scalars
729 : real(dp),intent(in) :: yy
730 : real(dp) :: derfc_yy
731 :
732 : !Local variables-------------------------------
733 : integer :: done,ii,isw
734 : ! coefficients for 0.0 <= yy < .477
735 : real(dp), parameter :: &
736 : & pp(5)=(/ 113.8641541510502e0_dp, 377.4852376853020e0_dp, &
737 : & 3209.377589138469e0_dp, .1857777061846032e0_dp, &
738 : & 3.161123743870566e0_dp /)
739 : real(dp), parameter :: &
740 : & qq(4)=(/ 244.0246379344442e0_dp, 1282.616526077372e0_dp, &
741 : & 2844.236833439171e0_dp, 23.60129095234412e0_dp/)
742 : ! coefficients for .477 <= yy <= 4.0
743 : real(dp), parameter :: &
744 : & p1(9)=(/ 8.883149794388376e0_dp, 66.11919063714163e0_dp, &
745 : & 298.6351381974001e0_dp, 881.9522212417691e0_dp, &
746 : & 1712.047612634071e0_dp, 2051.078377826071e0_dp, &
747 : & 1230.339354797997e0_dp, 2.153115354744038e-8_dp, &
748 : & .5641884969886701e0_dp /)
749 : real(dp), parameter :: &
750 : & q1(8)=(/ 117.6939508913125e0_dp, 537.1811018620099e0_dp, &
751 : & 1621.389574566690e0_dp, 3290.799235733460e0_dp, &
752 : & 4362.619090143247e0_dp, 3439.367674143722e0_dp, &
753 : & 1230.339354803749e0_dp, 15.74492611070983e0_dp/)
754 : ! coefficients for 4.0 < y,
755 : real(dp), parameter :: &
756 : & p2(6)=(/ -3.603448999498044e-01_dp, -1.257817261112292e-01_dp, &
757 : & -1.608378514874228e-02_dp, -6.587491615298378e-04_dp, &
758 : & -1.631538713730210e-02_dp, -3.053266349612323e-01_dp/)
759 : real(dp), parameter :: &
760 : & q2(5)=(/ 1.872952849923460e0_dp , 5.279051029514284e-01_dp, &
761 : & 6.051834131244132e-02_dp , 2.335204976268692e-03_dp, &
762 : & 2.568520192289822e0_dp /)
763 : real(dp), parameter :: &
764 : & sqrpi=.5641895835477563e0_dp, xbig=13.3e0_dp, xlarge=6.375e0_dp, xmin=1.0e-10_dp
765 : real(dp) :: res,xden,xi,xnum,xsq,xx
766 :
767 : !******************************************************************
768 :
769 1323947848 : xx = yy
770 1323947848 : isw = 1
771 : !Here change the sign of xx, and keep track of it thanks to isw
772 1323947848 : if (xx<0.0e0_dp) then
773 159608 : isw = -1
774 159608 : xx = -xx
775 : end if
776 :
777 1323947848 : done=0
778 :
779 : !Residual value, if yy < -6.375e0_dp
780 1323947848 : res=2.0e0_dp
781 :
782 : !abs(yy) < .477, evaluate approximation for erfc
783 1323947848 : if (xx<0.477e0_dp) then
784 : ! xmin is a very small number
785 123440 : if (xx<xmin) then
786 0 : res = xx*pp(3)/qq(3)
787 : else
788 123440 : xsq = xx*xx
789 123440 : xnum = pp(4)*xsq+pp(5)
790 123440 : xden = xsq+qq(4)
791 493760 : do ii = 1,3
792 370320 : xnum = xnum*xsq+pp(ii)
793 493760 : xden = xden*xsq+qq(ii)
794 : end do
795 123440 : res = xx*xnum/xden
796 : end if
797 123440 : if (isw==-1) res = -res
798 123440 : res = 1.0e0_dp-res
799 123440 : done=1
800 : end if
801 :
802 : !.477 < abs(yy) < 4.0 , evaluate approximation for erfc
803 1323947848 : if (xx<=4.0e0_dp .and. done==0 ) then
804 164750464 : xsq = xx*xx
805 164750464 : xnum = p1(8)*xx+p1(9)
806 164750464 : xden = xx+q1(8)
807 1318003712 : do ii=1,7
808 1153253248 : xnum = xnum*xx+p1(ii)
809 1318003712 : xden = xden*xx+q1(ii)
810 : end do
811 164750464 : res = xnum/xden
812 164750464 : res = res* exp(-xsq)
813 164750464 : if (isw.eq.-1) res = 2.0e0_dp-res
814 : done=1
815 : end if
816 :
817 : !y > 13.3e0_dp
818 1323923556 : if (isw > 0 .and. xx > xbig .and. done==0 ) then
819 5379318 : res = 0.0e0_dp
820 5379318 : done=1
821 : end if
822 :
823 : !4.0 < yy < 13.3e0_dp .or. -6.375e0_dp < yy < -4.0
824 : !evaluate minimax approximation for erfc
825 1323947848 : if ( ( isw > 0 .or. xx < xlarge ) .and. done==0 ) then
826 1153579474 : xsq = xx*xx
827 1153579474 : xi = 1.0e0_dp/xsq
828 1153579474 : xnum= p2(5)*xi+p2(6)
829 1153579474 : xden = xi+q2(5)
830 5767897370 : do ii = 1,4
831 4614317896 : xnum = xnum*xi+p2(ii)
832 5767897370 : xden = xden*xi+q2(ii)
833 : end do
834 1153579474 : res = (sqrpi+xi*xnum/xden)/xx
835 1153579474 : res = res* exp(-xsq)
836 1153579474 : if (isw.eq.-1) res = 2.0e0_dp-res
837 : end if
838 :
839 : !All cases have been investigated
840 1323947848 : derfc_yy = res
841 :
842 1323947848 : end function abi_derfc
843 : !!***
844 :
845 : !!****f* ABINIT/GAMMA_FUNCTION
846 : !! NAME
847 : !! GAMMA_FUNCTION
848 : !!
849 : !! FUNCTION
850 : !!
851 : !! INPUTS
852 : !!
853 : !! OUTPUT
854 : !!
855 : !! NOTES
856 : !!
857 : !! SOURCE
858 :
859 36 : subroutine GAMMA_FUNCTION(X,GA)
860 :
861 : ! ====================================================
862 : ! Purpose: This program computes the gamma function
863 : ! Gamma(x) using subroutine GAMMA
864 : ! Examples:
865 : ! x Gamma(x)
866 : ! ----------------------------
867 : ! 1/3 2.678938534708
868 : ! 0.5 1.772453850906
869 : ! -0.5 -3.544907701811
870 : ! -1.5 2.363271801207
871 : ! 5.0 24.000000000000
872 : ! ====================================================
873 : !
874 : ! This routine was downloaded from UIUC:
875 : ! http://jin.ece.uiuc.edu/routines/routines.html
876 : !
877 : ! The programs appear to accompany a book "Computation of Special
878 : ! Functions" (1996) John Wiley and Sons, but are distributed online
879 : ! by the authors. Exact copyright should be checked.
880 : !
881 : ! Authors / copyright:
882 : ! Shanjie Zhang and Jianming Jin
883 : ! Proposed contact is: j-jin1@uiuc.edu
884 : !
885 : ! 20 October 2008:
886 : ! Incorporated into ABINIT by M. Verstraete
887 : !
888 : !
889 : !
890 : ! ==================================================
891 : ! Purpose: Compute the gamma function Gamma(x)
892 : ! Input : x --- Argument of Gamma(x)
893 : ! ( x is not equal to 0,-1,-2, etc )
894 : ! Output: GA --- Gamma(x)
895 : ! ==================================================
896 : !
897 :
898 : ! arguments
899 :
900 : real(dp),intent(in) :: x
901 : real(dp),intent(out) :: ga
902 :
903 : ! local variables
904 : integer :: k,m
905 : real(dp) :: m1,z,r,gr
906 : real(dp) :: G(26)
907 :
908 : ! source code:
909 :
910 : ! initialization of reference data
911 : G=(/1.0D0,0.5772156649015329D0, &
912 : & -0.6558780715202538D0, -0.420026350340952D-1, &
913 : & 0.1665386113822915D0,-.421977345555443D-1, &
914 : & -.96219715278770D-2, .72189432466630D-2, &
915 : & -.11651675918591D-2, -.2152416741149D-3, &
916 : & .1280502823882D-3, -.201348547807D-4, &
917 : & -.12504934821D-5, .11330272320D-5, &
918 : & -.2056338417D-6, .61160950D-8, &
919 : & .50020075D-8, -.11812746D-8, &
920 : & .1043427D-9, .77823D-11, &
921 : & -.36968D-11, .51D-12, &
922 36 : & -.206D-13, -.54D-14, .14D-14, .1D-15/)
923 :
924 :
925 : ! for the integer case, do explicit factorial
926 36 : if (X==int(X)) then
927 0 : if (X > 0.0D0) then
928 0 : GA=1.0D0
929 0 : M1=X-1
930 0 : do K=2,int(M1)
931 0 : GA=GA*K
932 : end do
933 : else
934 0 : GA=1.0D+300
935 : end if
936 : ! for the integer case, do explicit factorial
937 : else
938 36 : if (abs(X) > 1.0D0) then
939 36 : Z=abs(X)
940 36 : M=int(Z)
941 36 : R=1.0D0
942 180 : do K=1,M
943 180 : R=R*(Z-K)
944 : end do
945 36 : Z=Z-M
946 : else
947 : Z=X
948 : end if
949 36 : GR=G(26)
950 936 : do K=25,1,-1
951 936 : GR=GR*Z+G(K)
952 : end do
953 36 : GA=1.0D0/(GR*Z)
954 36 : if (abs(X) > 1.0D0) then
955 36 : GA=GA*R
956 36 : if (X < 0.0D0) GA=-PI/(X*GA*SIN(PI*X))
957 : end if
958 : end if
959 36 : return
960 :
961 : end subroutine GAMMA_FUNCTION
962 : !!***
963 :
964 : !!****f* m_special_funcs/besjm
965 : !! NAME
966 : !! besjm
967 : !!
968 : !! FUNCTION
969 : !! Spherical bessel function of order nn. Handles nn=0,1,2,3,4, or 5 only.
970 : !!
971 : !! INPUTS
972 : !! arg= scaling to be applied to xx(nx)
973 : !! nn=order of spherical bessel function (only 0 through 5 allowed)
974 : !! cosx(1:nx)=cosines of arg*xx(1:nx)
975 : !! xx(1:nx)=set of dimensionless arguments of function
976 : !! nx=number of arguments
977 : !! sinx(1:nx)=sines of arg*xx(1:nx)
978 : !!
979 : !! OUTPUT
980 : !! besjx(1:nx)=returned values
981 : !!
982 : !! NOTES
983 : !! besj(nn,y)=$ j_{nn}(y) =(\frac{\pi}{2y})^{\frac{1}{2}}J(nn+\frac{1}{2},y)$
984 : !! where J=Bessel function of the first kind.
985 : !! besjm compute multiple values, and relies on precomputed values of sin and cos of y.
986 : !! The argument y is arg*xx(ix), for ix from 1 to nx
987 : !! The values of xx must be positive, and ordered by increasing order
988 : !! At small arg, the higher orders have so much cancellation that the
989 : !! analytic expression is very poor computationally. In that case we
990 : !! use a rational polynomial approximation.
991 : !!
992 : !! SOURCE
993 :
994 4807693 : subroutine besjm(arg,besjx,cosx,nn,nx,sinx,xx)
995 :
996 : !Arguments ------------------------------------
997 : !scalars
998 : integer,intent(in) :: nn,nx
999 : real(dp),intent(in) :: arg
1000 : !arrays
1001 : real(dp),intent(in) :: cosx(nx),sinx(nx),xx(nx)
1002 : real(dp),intent(out) :: besjx(nx)
1003 :
1004 : !Local variables-------------------------------
1005 : !scalars
1006 : integer :: ix,switchx
1007 : !Series or rational polynomial coefficients
1008 : real(dp),parameter :: b01=1.d0/6.d0,b02=1.d0/120.d0,b03=1.d0/5040.d0
1009 : real(dp),parameter :: b04=1.d0/362880.d0,b11=0.8331251468724171d-1
1010 : real(dp),parameter :: b12=0.2036961284395412d-2,b13=0.1932970379901801d-4
1011 : real(dp),parameter :: b14=0.6526053169009489d-7,b21=0.5867824627555163d-1
1012 : real(dp),parameter :: b22=0.1152501878595934d-2,b23=0.1011071389414764d-4
1013 : real(dp),parameter :: b24=0.4172322111421287d-7,b25=0.6790616688656543d-10
1014 : real(dp),parameter :: b31=0.439131885807176d-1,b32=0.6813139609887099d-3
1015 : real(dp),parameter :: b33=0.4899103784264755d-5,b34=0.17025590795625d-7
1016 : real(dp),parameter :: b35=0.2382642910613347d-10,b41=0.3587477991030971d-1
1017 : real(dp),parameter :: b42=0.4833719855268907d-3,b43=0.3238388977796242d-5
1018 : real(dp),parameter :: b44=0.1171802513125112d-7,b45=0.223261650431992d-10
1019 : real(dp),parameter :: b46=.1800045587335951d-13,b51=0.295232406376567d-1
1020 : real(dp),parameter :: b52=0.3359864457080573d-3,b53=0.19394750603618d-5
1021 : real(dp),parameter :: b54=0.6143166228216219d-8,b55=0.10378501636108d-10
1022 : real(dp),parameter :: b56=.749975122872713d-14
1023 : real(dp),parameter :: c11=0.1668748531275829d-1,c12=0.1342812442426702d-3
1024 : real(dp),parameter :: c13=0.6378249315355233d-6,c14=0.1573564527360138d-8
1025 : real(dp),parameter :: c21=0.127503251530198d-1,c22=0.7911240539893565d-4
1026 : real(dp),parameter :: c23=0.3044380758068054d-6,c24=0.7439837832363479d-9
1027 : real(dp),parameter :: c25=0.9515065658793124d-12,c31=0.1164236697483795d-1
1028 : real(dp),parameter :: c32=0.654858636312224d-4,c33=0.2265576367562734d-6
1029 : real(dp),parameter :: c34=0.4929905563217352d-9,c35=0.555120465710914d-12
1030 : real(dp),parameter :: c41=0.9579765544235745d-2,c42=0.4468999977536864d-4
1031 : real(dp),parameter :: c43=0.1315634305905896d-6,c44=0.2615492488301639d-9
1032 : real(dp),parameter :: c45=0.3387473312408129d-12,c46=.2280866204624012d-15
1033 : real(dp),parameter :: c51=0.8938297823881763d-2,c52=0.3874149021633025d-4
1034 : real(dp),parameter :: c53=0.1054692715135225d-6,c54=0.192879620987602d-9
1035 : real(dp),parameter :: c55=0.2284469423833734d-12,c56=0.139729234332572d-15
1036 : real(dp),parameter :: ffnth=1.d0/15.d0,o10395=1.d0/10395d0,oo105=1.d0/105.d0
1037 : real(dp),parameter :: oo945=1.d0/945.d0
1038 : real(dp) :: bot,rr,rsq,top
1039 : character(len=500) :: message
1040 : ! *************************************************************************
1041 :
1042 4807693 : if (nn==0) then
1043 :
1044 2370586 : switchx=nx+1
1045 644403935 : do ix=1,nx
1046 644394347 : rr=arg*xx(ix)
1047 644403935 : if (rr<=1.d-1) then
1048 642033349 : rsq=rr*rr
1049 642033349 : besjx(ix)=1.d0-rsq*(b01-rsq*(b02-rsq*(b03-rsq*b04)))
1050 : else
1051 : switchx=ix
1052 : exit
1053 : end if
1054 : end do
1055 :
1056 1321653242 : do ix=switchx,nx
1057 1319282656 : rr=arg*xx(ix)
1058 1321653242 : besjx(ix)=sinx(ix)/rr
1059 : end do
1060 :
1061 : else if (nn==1) then
1062 :
1063 1869755 : switchx=nx+1
1064 782041558 : do ix=1,nx
1065 781972309 : rr=arg*xx(ix)
1066 782041558 : if (rr<=1.d0) then
1067 780171803 : rsq=rr*rr
1068 780171803 : top=1.d0-rsq*(b11-rsq*(b12-rsq*(b13-rsq*b14)))
1069 780171803 : bot=1.d0+rsq*(c11+rsq*(c12+rsq*(c13+rsq*c14)))
1070 780171803 : besjx(ix)=third*rr*top/bot
1071 : else
1072 : switchx=ix
1073 : exit
1074 : end if
1075 : end do
1076 :
1077 783956221 : do ix=switchx,nx
1078 782086466 : rr=arg*xx(ix)
1079 782086466 : rsq=rr*rr
1080 783956221 : besjx(ix)=(sinx(ix)-rr*cosx(ix))/rsq
1081 : end do
1082 :
1083 : else if (nn==2) then
1084 :
1085 555190 : switchx=nx+1
1086 179572647 : do ix=1,nx
1087 179552563 : rr=arg*xx(ix)
1088 179572647 : if (rr<=2.d0) then
1089 179017457 : rsq=rr*rr
1090 179017457 : top=1.d0-rsq*(b21-rsq*(b22-rsq*(b23-rsq*(b24-rsq*b25))))
1091 179017457 : bot=1.d0+rsq*(c21+rsq*(c22+rsq*(c23+rsq*(c24+rsq*c25))))
1092 179017457 : besjx(ix)=ffnth*rsq*top/bot
1093 : else
1094 : switchx=ix
1095 : exit
1096 : end if
1097 : end do
1098 :
1099 263961161 : do ix=switchx,nx
1100 263405971 : rr=arg*xx(ix)
1101 263405971 : rsq=rr*rr
1102 263961161 : besjx(ix)=((3.d0-rsq)*sinx(ix)-3.d0*rr*cosx(ix))/(rr*rsq)
1103 : end do
1104 :
1105 : else if (nn==3) then
1106 :
1107 12081 : switchx=nx+1
1108 8377955 : do ix=1,nx
1109 8377180 : rr=arg*xx(ix)
1110 8377955 : if (rr<=2.d0) then
1111 8365874 : rsq=rr*rr
1112 8365874 : top=1.d0-rsq*(b31-rsq*(b32-rsq*(b33-rsq*(b34-rsq*b35))))
1113 8365874 : bot=1.d0+rsq*(c31+rsq*(c32+rsq*(c33+rsq*(c34+rsq*c35))))
1114 8365874 : besjx(ix)=rr*rsq*oo105*top/bot
1115 : else
1116 : switchx=ix
1117 : exit
1118 : end if
1119 : end do
1120 :
1121 3651507 : do ix=switchx,nx
1122 3639426 : rr=arg*xx(ix)
1123 3639426 : rsq=rr*rr
1124 : besjx(ix)=( (15.d0-6.d0*rsq)*sinx(ix)&
1125 3651507 : & + rr*(rsq-15.d0) *cosx(ix) ) /(rsq*rsq)
1126 : end do
1127 :
1128 : else if (nn==4) then
1129 :
1130 81 : switchx=nx+1
1131 3402 : do ix=1,nx
1132 3402 : rr=arg*xx(ix)
1133 3402 : if (rr<=4.d0) then
1134 3321 : rsq=rr*rr
1135 3321 : top=1.d0-rsq*(b41-rsq*(b42-rsq*(b43-rsq*(b44-rsq*(b45-rsq*b46)))))
1136 3321 : bot=1.d0+rsq*(c41+rsq*(c42+rsq*(c43+rsq*(c44+rsq*(c45+rsq*c46)))))
1137 3321 : besjx(ix)=rsq*rsq*oo945*top/bot
1138 : else
1139 : switchx=ix
1140 : exit
1141 : end if
1142 : end do
1143 :
1144 38060 : do ix=switchx,nx
1145 37979 : rr=arg*xx(ix)
1146 37979 : rsq=rr*rr
1147 : besjx(ix)=( (105.d0-rsq*(45.d0-rsq)) *sinx(ix)&
1148 38060 : & + rr * (10.d0*rsq-105.d0) *cosx(ix) ) /(rsq*rsq*rr)
1149 : end do
1150 :
1151 : else if (nn==5) then
1152 :
1153 0 : switchx=nx+1
1154 0 : do ix=1,nx
1155 0 : rr=arg*xx(ix)
1156 0 : if (rr<=4.d0) then
1157 0 : rsq=rr*rr
1158 0 : top=1.d0-rsq*(b51-rsq*(b52-rsq*(b53-rsq*(b54-rsq*(b55-rsq*b56)))))
1159 0 : bot=1.d0+rsq*(c51+rsq*(c52+rsq*(c53+rsq*(c54+rsq*(c55+rsq*c56)))))
1160 0 : besjx(ix)=rsq*rsq*rr*o10395*top/bot
1161 : else
1162 : switchx=ix
1163 : exit
1164 : end if
1165 : end do
1166 :
1167 0 : do ix=switchx,nx
1168 0 : rr=arg*xx(ix)
1169 0 : rsq=rr*rr
1170 : besjx(ix)=( (945.d0-rsq*(420.d0-rsq*15.d0)) *sinx(ix)&
1171 0 : & + rr * (945.d0-rsq*(105.d0-rsq)) *cosx(ix) ) /(rsq*rsq*rr)
1172 : end do
1173 :
1174 : else
1175 0 : write(message, '(a,i0,a)' )' besjm only defined for nn in [0,5]; input was nn=',nn,'.'
1176 0 : ABI_BUG(message)
1177 : end if
1178 :
1179 4807693 : end subroutine besjm
1180 : !!***
1181 :
1182 : !!****f* m_special_funcs/sbf8
1183 : !! NAME
1184 : !! sbf8
1185 : !!
1186 : !! FUNCTION
1187 : !! Computes set of spherical bessel functions using accurate algorithm
1188 : !! based on downward recursion in order and normalization sum.
1189 : !! Power series used at small arguments.
1190 : !!
1191 : !! INPUTS
1192 : !! nm=maximum angular momentum wanted + one
1193 : !! xx=argument of sbf
1194 : !!
1195 : !! OUTPUT
1196 : !! sb_out(nm)=values of spherical bessel functions for l=0,nm-1
1197 : !!
1198 : !! SOURCE
1199 :
1200 812315248 : subroutine sbf8(nm,xx,sb_out)
1201 :
1202 : !Arguments----------------------------------------------------------
1203 : !scalars
1204 : integer,intent(in) :: nm
1205 : real(dp),intent(in) :: xx
1206 : !arrays
1207 : real(dp),intent(out) :: sb_out(nm)
1208 :
1209 : !Local variables-------------------------------
1210 : !scalars
1211 : integer :: nlim,nn
1212 : real(dp) :: fn,sn,xi,xn,xs
1213 : !arrays
1214 812315248 : real(dp),allocatable :: sb(:)
1215 : ! *************************************************************************
1216 :
1217 812315248 : if(xx<= 1.0e-36_dp) then
1218 : ! zero argument section
1219 7827685 : sb_out(:)=zero
1220 2061565 : sb_out(1)=one
1221 810253683 : else if(xx<1.e-3_dp) then
1222 : ! small argument section
1223 15642225 : xn=one
1224 15642225 : xs=half*xx**2
1225 63748618 : do nn=1,nm
1226 48106393 : sb_out(nn)=xn*(one - xs*(one - xs/(4*nn+6))/(2*nn+1))
1227 63748618 : xn=xx*xn/(2*nn+1)
1228 : end do
1229 : else
1230 : ! recursion method
1231 794611458 : if(xx<one) then
1232 366068511 : nlim=nm+int(15.0e0_dp*xx)+1
1233 : else
1234 428542947 : nlim=nm+int(1.36e0_dp*xx)+15
1235 : end if
1236 2383834374 : ABI_MALLOC(sb,(nlim+1))
1237 794611458 : nn=nlim
1238 794611458 : xi=one/xx
1239 794611458 : sb(nn+1)=zero
1240 794611458 : sb(nn)=1.e-18_dp
1241 794611458 : sn=dble(2*nn-1)*1.e-36_dp
1242 14384210566 : do nn=nlim-1,1,-1
1243 14384210566 : sb(nn)=dble(2*nn+1)*xi*sb(nn+1) - sb(nn+2)
1244 : end do
1245 14384210566 : do nn=1,nlim-1
1246 14384210566 : sn=sn + dble(2*nn-1)*sb(nn)*sb(nn)
1247 : end do
1248 794611458 : fn=1.d0/sqrt(sn)
1249 3103206171 : sb_out(:)=fn*sb(1:nm)
1250 794611458 : ABI_FREE(sb)
1251 : end if
1252 :
1253 812315248 : end subroutine sbf8
1254 : !!***
1255 :
1256 : !----------------------------------------------------------------------
1257 :
1258 : !!****f* m_special_funcs/fermi_dirac
1259 : !! NAME
1260 : !! fermi_dirac
1261 : !!
1262 : !! FUNCTION
1263 : !! Returns the Fermi Dirac distribution for T and energy wrt Fermi level
1264 : !! presumes everything is in Hartree!!!! Not Kelvin for T
1265 : !!
1266 : !! INPUTS
1267 : !! energy = electron energy level
1268 : !! mu = chemical potential
1269 : !! temperature = T
1270 : !!
1271 : !! SOURCE
1272 :
1273 11046015 : function fermi_dirac(energy, mu, temperature)
1274 :
1275 : !Arguments ------------------------------------
1276 : !scalars
1277 : real(dp),intent(in) :: energy, mu, temperature
1278 : real(dp) :: fermi_dirac
1279 :
1280 : !Local variables-------------------------------
1281 : !scalars
1282 : real(dp) :: arg
1283 : ! *************************************************************************
1284 :
1285 11046015 : fermi_dirac = zero
1286 11046015 : if (temperature > tol12) then
1287 11039622 : arg = (energy-mu)/temperature
1288 11039622 : if(arg < -600._dp)then ! far below Ef
1289 : fermi_dirac = one
1290 11039622 : else if (arg < 600._dp)then ! around Ef
1291 11039622 : fermi_dirac = one / (exp(arg) + one)
1292 : end if
1293 : else ! T is too small - just step function
1294 6393 : if (mu-energy > tol12) fermi_dirac = one
1295 : end if
1296 :
1297 11046015 : end function fermi_dirac
1298 : !!***
1299 :
1300 : !----------------------------------------------------------------------
1301 :
1302 : !!****f* m_special_funcs/bose_einstein
1303 : !! NAME
1304 : !! bose_einstein
1305 : !!
1306 : !! FUNCTION
1307 : !! Returns the Bose Einstein distribution for T and energy
1308 : !! presumes everything is in Hartree!!!! Not Kelvin for T
1309 : !!
1310 : !! INPUTS
1311 : !! energy = electron energy level
1312 : !! temperature = T
1313 : !!
1314 : !! SOURCE
1315 :
1316 1247510 : real(dp) function bose_einstein(energy, temperature)
1317 :
1318 : !Arguments ------------------------------------
1319 : !scalars
1320 : real(dp),intent(in) :: energy, temperature
1321 :
1322 : !Local variables-------------------------------
1323 : real(dp) :: arg
1324 : ! *************************************************************************
1325 :
1326 1247510 : bose_einstein = zero
1327 1247510 : if (temperature > tol12) then
1328 1241117 : arg = energy/temperature
1329 1241117 : if(arg > tol12 .and. arg < 600._dp)then
1330 1102913 : bose_einstein = one / (exp(arg) - one)
1331 138204 : else if (arg < tol12) then
1332 138204 : ABI_WARNING('No Bose Einstein for negative energies')
1333 : end if
1334 : else
1335 6393 : ABI_WARNING('No Bose Einstein for negative or 0 T')
1336 : end if
1337 :
1338 :
1339 1247510 : end function bose_einstein
1340 : !!***
1341 :
1342 : !----------------------------------------------------------------------
1343 :
1344 : !!****f* m_special_funcs/dip12
1345 : !! NAME
1346 : !! dip12
1347 : !!
1348 : !! FUNCTION
1349 : !! Returns the complete Fermi integral of order 1/2.
1350 : !! Based on an analytical approximation.
1351 : !!
1352 : !! INPUTS
1353 : !! gamma=complete Fermi integral argument
1354 : !!
1355 : !! OUTPUT
1356 : !! dip12=resulting function
1357 : !!
1358 : !! SOURCE
1359 :
1360 158 : real(dp) function dip12(gamma)
1361 :
1362 : ! Arguments -------------------------------
1363 : ! Scalars
1364 : real(dp),intent(in) :: gamma
1365 :
1366 : ! Local variables -------------------------
1367 : ! Scalars
1368 : real(dp) :: d,dy
1369 : ! *********************************************************************
1370 :
1371 158 : if (gamma.lt.3.) then
1372 79 : dy=exp(gamma)
1373 79 : if (gamma+1.9375.LE.0) then
1374 : dip12=dy*&
1375 : & (1.-dy*(0.35355283-dy*(0.19242767-dy*(0.12456909-dy*&
1376 79 : & (8.5114507E-02-dy*4.551794E-02)))))
1377 : else
1378 0 : d=gamma-0.5
1379 : dip12=dy*(0.677695804-d*(0.187773135+d*(2.16197521E-02-d*&
1380 : & (9.23703807E-03+d*&
1381 : & (1.71735167E-03-d*(6.07913775E-04+d*&
1382 : & (1.1448629E-04-d*&
1383 : & (4.544432E-05+d*(6.4719368E-06-d*(3.794983E-06+d*&
1384 : & (1.7338029E-07-d*&
1385 : & (3.5546516E-07-d*(3.7329191E-08+d*&
1386 : & (3.3097822E-08-d*&
1387 : & (8.3190193E-09+d*(2.2752769E-09-d*(7.836005E-10+d*&
1388 0 : & (7.519551E-11-d*2.960006E-11))))))))))))))))))
1389 : end if
1390 79 : else if (gamma.lt.20.) then
1391 0 : if (gamma.lt.10.) then
1392 0 : d=gamma-6.5
1393 : dip12=12.839811+d*&
1394 : & (2.844774+d*(0.114920926-d*(3.43733039E-03-d*&
1395 : & (2.3980356E-04-d*&
1396 : & (2.0201888E-05-d*(1.5219883E-06-d*&
1397 : & (6.2770524E-08+d*&
1398 : & (4.8830336E-09-d*(2.1031164E-09-d*(5.785753E-10-d*&
1399 0 : & (7.233066E-11-d*1.230727E-12)))))))))))
1400 : else
1401 0 : d=gamma-14.5
1402 : dip12=41.7799227+d*&
1403 : & (4.2881461+d*(7.45407825E-02-d*(8.79243296E-04-d*&
1404 : & (2.38288861E-05-d*&
1405 : & (8.82474867E-07-d*(3.82865217E-08-d*&
1406 0 : & (1.9274292E-09-d*(1.42248669E-10-d*8.17019813E-12))))))))
1407 : end if
1408 : else
1409 79 : d=1./gamma
1410 79 : dy=gamma*dsqrt(gamma)/1.329340388
1411 : dip12=dy*(1.-d*(9.354E-07-d*(1.2338391-d*(6.77931E-03-d*&
1412 79 : & 1.17871643))))
1413 : end if
1414 158 : dip12=dip12*0.88622692
1415 158 : end function dip12
1416 : !!***
1417 :
1418 : !----------------------------------------------------------------------
1419 :
1420 : !!****f* m_special_funcs/dip32
1421 : !! NAME
1422 : !! dip32
1423 : !!
1424 : !! FUNCTION
1425 : !! Returns the complete Fermi integral of order 3/2.
1426 : !! Based on an analytical approximation.
1427 : !!
1428 : !! INPUTS
1429 : !! gamma=complete Fermi integral argument
1430 : !!
1431 : !! OUTPUT
1432 : !! dip32=resulting function
1433 : !!
1434 : !! SOURCE
1435 :
1436 79 : real(dp) function dip32(gamma)
1437 :
1438 : ! Arguments -------------------------------
1439 : real(dp),intent(in) :: gamma
1440 :
1441 : ! Local variables -------------------------
1442 : ! Scalars
1443 : real(dp) :: d,dval
1444 : ! *********************************************************************
1445 :
1446 79 : if (gamma.GT.1.75) then
1447 0 : dval=gamma*gamma*dsqrt(gamma)
1448 0 : if (gamma.LT.4.5) then
1449 0 : d=gamma-3.125
1450 : dip32=(1.27623+0.596065*gamma+0.3*dval)*&
1451 : & (1.0055558385-d*(5.23889494E-03+d*&
1452 : & (3.13523144E-03-d*(3.06124286E-03-d*&
1453 : & (1.3644667E-03-d*&
1454 : & (4.1528384E-04-d*(8.901188E-05-d*(1.079979E-05+d*&
1455 0 : & (2.29058E-06-d*(2.58985E-06-d*7.30909E-07))))))))))
1456 0 : else if (gamma.LT.12.) then
1457 0 : if (gamma.LT.8.) then
1458 0 : d=gamma-6.25
1459 : dip32=(2.01508+0.425775*gamma+0.3*dval)*&
1460 : & (1.000387131-d*(3.93626295E-04+d*&
1461 : & (2.55710115E-04-d*&
1462 : & (1.57383494E-04-d*(5.0286036E-05-d*&
1463 : & (1.2073559865E-05-d*&
1464 : & (2.4909523213E-06-d*(5.244328548E-07-d*&
1465 0 : & 8.0884033896E-08))))))))
1466 : else
1467 0 : d=gamma-10.
1468 : dip32=0.3*dval*&
1469 : & (1.064687247-d*(1.22972303E-02-d*(1.8362121E-03-d*&
1470 : & (2.433558E-04-d*(3.018186E-05-d*(3.5694E-06-d*&
1471 0 : & (4.11212E-07-d*(5.2151E-08-d*5.8424E-09))))))))
1472 : end if
1473 : else
1474 0 : d=1./gamma
1475 : dip32=0.30090111127*dval*&
1476 : & (1.-d*(2.863E-06-d*(6.168876549-d*&
1477 0 : & (1.740553E-02+d*(1.425257+d*2.95887)))))
1478 : end if
1479 79 : else if (gamma+0.75.LE.0) then
1480 79 : d=EXP(gamma)
1481 : dip32=d*&
1482 : & (1.-d*(1.76775246E-01-d*(6.4124584E-02-d*&
1483 : & (3.1027055E-02-d*(1.6797637E-02-d*&
1484 79 : & (8.212636E-03-d*(2.384106E-03)))))))
1485 : else
1486 0 : d=gamma-0.5
1487 : dip32=EXP(gamma)*(0.846691-0.128948*gamma)*&
1488 : & (1.034064158+d*(2.778947E-03-d*&
1489 : & (3.572502805E-02+d*(3.0411645E-03-d*&
1490 : & (1.7380548E-03+d*(2.7756776E-04-d*&
1491 : & (8.08302E-05+d*(1.59606E-05-d*&
1492 0 : & (3.8144E-06+d*7.4446E-07)))))))))
1493 : end if
1494 79 : dip32=dip32*1.32934038
1495 :
1496 79 : end function dip32
1497 : !!***
1498 :
1499 : !----------------------------------------------------------------------
1500 :
1501 : !!****f* m_special_funcs/djp12
1502 : !! NAME
1503 : !! djp12
1504 : !!
1505 : !! FUNCTION
1506 : !! Returns the incomplete Fermi integral of order 1/2.
1507 : !! Based on an analytical approximation.
1508 : !!
1509 : !! INPUTS
1510 : !! xcut=lower bound of the incomplete Fermi integral
1511 : !! gamma=incomplete Fermi integral argument
1512 : !!
1513 : !! OUTPUT
1514 : !! djp12=resulting function
1515 : !!
1516 : !! SOURCE
1517 :
1518 4149 : real(dp) function djp12(xcut, gamma)
1519 :
1520 : ! Arguments -------------------------------
1521 : ! Scalars
1522 : real(dp),intent(in) :: xcut,gamma
1523 :
1524 : ! Local variables -------------------------
1525 : ! Scalars
1526 : real(dp) :: d2h,db,dc,dd,de,df1,df2,df3,dh
1527 : real(dp) :: ds,dt,dv,dw,dxm,dxp
1528 : integer :: i,ind,iq,k,nm,np,nq
1529 : ! Arrays
1530 : real(dp) :: dq(5),df(101),dy(101)
1531 : ! *********************************************************************
1532 :
1533 4149 : dh=0.2D+0
1534 4149 : d2h=0.4D+0
1535 4149 : nm=101
1536 4149 : ind=0
1537 4149 : dq=(/1.D+0,2.828427124D+0,5.196152423D+0,8.D+0,1.118033989D+1/)
1538 :
1539 4149 : djp12=0.D+0
1540 4149 : dxm=gamma-1.5D+1
1541 4149 : if (xcut.gt.dxm) then
1542 : if (ind.eq.0) then
1543 415140 : do i=1,nm
1544 411070 : dy(i)=-1.5D+1+(i-1)*dh
1545 415140 : df(i)=1.D+0+dexp(dy(i))
1546 : end do
1547 4070 : ind=1
1548 : end if
1549 4070 : dxp=gamma+5.D+0
1550 : if (xcut.lt.dxp) then
1551 : dc=dxp
1552 : else
1553 : dc=xcut
1554 : end if
1555 4070 : db=dexp(gamma-dc)
1556 4070 : dt=db
1557 15892 : do iq=1,5
1558 15892 : dd=iq*dc
1559 15892 : ds=dsqrt(dd)
1560 15892 : dw=1.+.3275911*ds
1561 15892 : dw=1.D+0/dw
1562 : dv=dw*(.2258368458D+0+&
1563 : & dw*(-.2521286676D+0+dw*(1.2596951294D+0+&
1564 15892 : & dw*(-1.2878224530D+0+dw*(.9406460699D+0)))))
1565 15892 : dv=dv+ds
1566 15892 : de=dt*dv/dq(iq)
1567 15892 : djp12=djp12+de
1568 15892 : if (dabs(de).lt.(1.D-07*djp12)) exit
1569 15892 : dt=-dt*db
1570 : end do
1571 4070 : if (xcut.ge.dxp) return
1572 3507 : np=(dxp-xcut)/dh
1573 3507 : np=2*(np/2)
1574 3507 : np=nm-np
1575 3507 : nq=(15.-gamma)/dh
1576 3507 : nq=1+2*(nq/2)
1577 3507 : if (np.lt.nq) np=nq
1578 3507 : if (np.le.nm) then
1579 3507 : df3=0.D+0
1580 3507 : dt=dy(np)+gamma
1581 3507 : dv=(dt-xcut)/2.D+0
1582 3507 : df3=0.D0
1583 3507 : if (dt.ge.1.D-13) df3=dsqrt(dt)
1584 3507 : if (dabs(dv).ge.1.D-13) then
1585 3507 : df1=dsqrt(xcut)
1586 3507 : dt=df1+df3
1587 3507 : dw=(dv+dv)/(dt*dt)
1588 3507 : df2=dw*dw
1589 3507 : df2=df2+df2
1590 3507 : db=df2*(df2+7.D+0)+7.D+1
1591 3507 : dc=7.D+0*(1.D+1-df2)
1592 3507 : dc=dc*dw
1593 3507 : dd=-df2*(df2-2.8D+1)+1.4D+2
1594 3507 : dd=dd+dd
1595 : ds=dt*((db-dc)/(1.D+0+dexp(xcut-gamma))+&
1596 3507 : & dd/(1.D+0+dexp(xcut+dv-gamma))+(db+dc)/df(np))
1597 3507 : ds=ds*dv/4.2D+2
1598 3507 : djp12=djp12+ds
1599 : end if
1600 3507 : if (np.ne.nm) then
1601 3507 : ds=0.D+0
1602 3507 : np=np+2
1603 3507 : do k=np,nm,2
1604 7851 : df1=df3
1605 7851 : df3=dsqrt(dy(k)+gamma)
1606 7851 : dt=df1+df3
1607 7851 : dw=d2h/(dt*dt)
1608 7851 : df2=dw*dw
1609 7851 : df2=df2+df2
1610 7851 : db=df2*(df2+7.D+0)+7.D+1
1611 7851 : dc=7.D+0*(1.D+1-df2)
1612 7851 : dc=dc*dw
1613 7851 : dd=-df2*(df2-2.8D+1)+1.4D+2
1614 7851 : dd=dd+dd
1615 7851 : ds=ds+dt*((db-dc)/df(k-2)+dd/df(k-1)+(db+dc)/df(k))
1616 : end do
1617 3507 : ds=ds*dh/4.2D+2
1618 3507 : djp12=djp12+ds
1619 : end if
1620 3507 : if (xcut.ge.dxm) return
1621 : end if
1622 : end if
1623 79 : djp12=dip12(gamma)-xcut*dsqrt(xcut)/1.5D+0
1624 :
1625 79 : end function djp12
1626 : !!***
1627 :
1628 : !----------------------------------------------------------------------
1629 :
1630 : !!****f* m_special_funcs/djp32
1631 : !! NAME
1632 : !! djp32
1633 : !!
1634 : !! FUNCTION
1635 : !! Returns the incomplete Fermi integral of order 3/2.
1636 : !! Based on an analytical approximation.
1637 : !!
1638 : !! INPUTS
1639 : !! xcut=lower bound of the incomplete Fermi integral
1640 : !! gamma=incomplete Fermi integral argument
1641 : !!
1642 : !! OUTPUT
1643 : !! djp32=resulting function
1644 : !!
1645 : !! SOURCE
1646 :
1647 79 : real(dp) function djp32(xcut,gamma)
1648 :
1649 : ! Arguments -------------------------------
1650 : ! Scalars
1651 : real(dp),intent(in) :: xcut,gamma
1652 :
1653 : ! Local variables -------------------------
1654 : ! Scalars
1655 : real(dp) :: d2h,db,dc,dd,de,df1,df2,df3,dh
1656 : real(dp) :: ds,dt,dv,dw,dx1,dx2
1657 : real(dp) :: dx3,dxm,dxp
1658 : integer :: i,ind,iq,k,nm,np,nq
1659 : ! Arrays
1660 : real(dp) :: dq(5),df(101),dy(101)
1661 : ! *********************************************************************
1662 :
1663 79 : dh=0.2D+0
1664 79 : d2h=0.4D+0
1665 79 : nm=101
1666 79 : ind=0
1667 79 : dq=(/1.D+0,5.656854228D+0,1.558845727D+1,3.2D+1,5.590169945D+0/)
1668 :
1669 79 : djp32=0.D+0
1670 79 : dxm=gamma-1.5D+1
1671 79 : if (xcut.GT.dxm) then
1672 : if (ind.EQ.0) then
1673 8058 : do i=1,nm
1674 7979 : dy(i)=-1.5D+1+(i-1)*dh
1675 8058 : df(i)=1.D+0+DEXP(dy(i))
1676 : end do
1677 79 : ind=1
1678 : end if
1679 79 : dxp=gamma+5.D+0
1680 : if (xcut.LT.dxp) then
1681 : dc=dxp
1682 : else
1683 : dc=xcut
1684 : end if
1685 79 : db=DEXP(gamma-dc)
1686 79 : dt=db
1687 309 : do iq=1,5
1688 309 : dd=iq*dc
1689 309 : ds=DSQRT(dd)
1690 309 : dw=1.+.3275911*ds
1691 309 : dw=1.D+0/dw
1692 : dv=dw*(.2258368458D+0+&
1693 : & dw*(-.2521286676D+0+dw*(1.2596951294D+0+&
1694 309 : & dw*(-1.2878224530D+0+dw*(.9406460699D+0)))))
1695 309 : dv=dv+ds
1696 309 : dv=1.5D+0*dv+ds*dd
1697 309 : de=dt*dv/dq(iq)
1698 309 : djp32=djp32+de
1699 309 : if (DABS(de).LT.(1.D-07*djp32)) exit
1700 309 : dt=-dt*db
1701 : end do
1702 79 : if (xcut.GE.dxp) return
1703 72 : np=(dxp-xcut)/dh
1704 72 : np=2*(np/2)
1705 72 : np=nm-np
1706 72 : nq=(15.-gamma)/dh
1707 72 : nq=1+2*(nq/2)
1708 72 : if (np.LT.nq) np=nq
1709 72 : if (np.LE.nm) then
1710 72 : df3=0.D+0
1711 72 : dt=dy(np)+gamma
1712 72 : dv=(dt-xcut)/2.D+0
1713 72 : df3=DSQRT(dt)
1714 72 : dx3=dt
1715 72 : if (DABS(dv).GE.1.D-13) then
1716 72 : df1=DSQRT(xcut)
1717 72 : dt=df1+df3
1718 72 : dw=(dv+dv)/(dt*dt)
1719 72 : df2=dw*dw
1720 72 : df2=df2+df2
1721 72 : db=df2*(df2+7.D+0)+7.D+1
1722 72 : dc=7.D+0*(1.D+1-df2)
1723 72 : dc=dc*dw
1724 72 : dd=-df2*(df2-2.8D+1)+1.4D+2
1725 72 : dd=dd+dd
1726 : ds=dt*((db-dc)*xcut/(1.D+0+DEXP(xcut-gamma))+dd*(xcut+dv)&
1727 72 : & /(1.D+0+DEXP(xcut+dv-gamma))+(db+dc)*(dy(np)+gamma)/df(np))
1728 72 : ds=ds*dv/4.2D+2
1729 72 : djp32=djp32+ds
1730 : end if
1731 72 : if (np.NE.nm) then
1732 72 : ds=0.D+0
1733 72 : np=np+2
1734 72 : do k=np,nm,2
1735 144 : dx1=dx3
1736 144 : df1=df3
1737 144 : dx2=dy(k-1)+gamma
1738 144 : dx3=dy(k)+gamma
1739 144 : df3=DSQRT(dx3)
1740 144 : dt=df1+df3
1741 144 : dw=d2h/(dt*dt)
1742 144 : df2=dw*dw
1743 144 : df2=df2+df2
1744 144 : db=df2*(df2+7.D+0)+7.D+1
1745 144 : dc=7.D+0*(1.D+1-df2)
1746 144 : dc=dc*dw
1747 144 : dd=-df2*(df2-2.8D+1)+1.4D+2
1748 144 : dd=dd+dd
1749 : ds=ds+dt*((db-dc)*dx1/df(k-2)+dd*dx2/df(k-1)&
1750 144 : & +(db+dc)*dx3/df(k))
1751 : end do
1752 72 : ds=ds*dh/4.2D+2
1753 72 : djp32=djp32+ds
1754 : end if
1755 72 : if (xcut.GE.dxm) return
1756 : end if
1757 : end if
1758 0 : djp32=dip32(gamma)-xcut*xcut*DSQRT(xcut)/2.5D+0
1759 :
1760 0 : end function djp32
1761 : !!***
1762 :
1763 : !!****f* m_special_funcs/tildeAx
1764 : !! NAME
1765 : !! tildeAx
1766 : !!
1767 : !! FUNCTION
1768 : !! Returns tilde Ax Pade fit and first and second derivatives
1769 : !! w.r.t. reduced temperature t.
1770 : !!
1771 : !! NOTES
1772 : !! Karasiev-Sjostrom-Dufty-Trickey (KSDT) TLDA xc-functional
1773 : !! V.V. Karasiev, T. Sjostrom, J. Dufty, and S.B. Trickey, PRL 112, 076403 (2014) [[cite:Karasiev2014]]
1774 : !!
1775 : !! INPUTS
1776 : !! t=reduced temperature, t=T/T_F
1777 : !!
1778 : !! OUTPUT
1779 : !! Ax=tilde Ax(t)
1780 : !! dAx=dAx(t)/dt
1781 : !! d2Ax=d^2Ax(t)/dt^2
1782 : !!
1783 : !! SOURCE
1784 :
1785 129516 : subroutine tildeAx(t,Ax,dAx,d2Ax)
1786 :
1787 : !Arguments ------------------------------------
1788 : !scalars
1789 : real(dp),intent(in) :: t
1790 : real(dp),intent(out) :: Ax,dAx,d2Ax
1791 : !Local variables ------------------------------
1792 : !scalars
1793 : real(dp),parameter :: aln= -0.0475410604245741_dp
1794 : real(dp),parameter :: a52= -0.1065378473507800_dp
1795 : real(dp),parameter :: a1 = 0.5823869764908659_dp
1796 : real(dp),parameter :: a2 = -0.0068339509356661_dp
1797 : real(dp),parameter :: a3 = 11.5469239288490009_dp
1798 : real(dp),parameter :: a4 = -0.8465428870889800_dp
1799 : real(dp),parameter :: a5 = -0.1212525366470300_dp
1800 : real(dp),parameter :: a6 = 1.9902818786101000_dp
1801 : real(dp),parameter :: a7 = 0.0000000000000000_dp
1802 : real(dp),parameter :: a8 = 0.0744389046707120_dp
1803 : real(dp),parameter :: b1 = 19.9256144707979992_dp
1804 : real(dp),parameter :: b2 = 5.1663994545590004_dp
1805 : real(dp),parameter :: b3 = 2.0463164858237000_dp
1806 : real(dp),parameter :: b4 = 0.0744389046707120_dp
1807 : real(dp),parameter :: onethird = 1._dp/3._dp
1808 : real(dp),parameter :: twothird = 2._dp/3._dp
1809 : real(dp),parameter :: fourthird = 4._dp/3._dp
1810 : real(dp),parameter :: threehalf = 3._dp/2._dp
1811 : real(dp),parameter :: fivehalf = 5._dp/2._dp
1812 : real(dp),parameter :: sevenhalf = 7._dp/2._dp
1813 : real(dp) :: y,u,du,d2u
1814 : real(dp) :: v,dv,d2v
1815 : real(dp) :: dydt,d2ydt2
1816 : real(dp) :: num,den,fit,dnum,d2num,dden,d2den,dfit,d2fit
1817 : ! *************************************************************************
1818 :
1819 129516 : y = twothird/t**threehalf
1820 129516 : u = y**twothird
1821 129516 : du = twothird/y**onethird
1822 129516 : d2u = -onethird*du/y
1823 :
1824 129516 : v = y**fourthird
1825 129516 : dv = fourthird *y**onethird
1826 129516 : d2v = onethird*dv/y
1827 :
1828 129516 : dydt = -1._dp/t**fivehalf
1829 129516 : d2ydt2 = fivehalf/t**sevenhalf
1830 :
1831 : num = a52*u**fivehalf &
1832 : +a1*u+a2*u**2+a3*u**3+a4*u**4 &
1833 : +a5*u**5+a6*u**6+a7*u**7+a8*u**8 &
1834 129516 : +aln*log(y)*y**4
1835 129516 : den = 1._dp+b1*v+b2*v**2+b3*v**3+b4*v**4
1836 129516 : fit = num/den
1837 : dnum = du*(fivehalf*a52*u**threehalf &
1838 : +a1+2._dp*a2*u+3._dp*a3*u**2+4._dp*a4*u**3 &
1839 : +5._dp*a5*u**4+6._dp*a6*u**5+7._dp*a7*u**6+8._dp*a8*u**7) &
1840 129516 : +aln*y**3+4._dp*aln*log(y)*y**3
1841 : d2num = d2u*(fivehalf*a52*u**threehalf &
1842 : +a1+2._dp*a2*u+3._dp*a3*u**2+4._dp*a4*u**3 &
1843 : +5._dp*a5*u**4+6._dp*a6*u**5+7._dp*a7*u**6+8._dp*a8*u**7) &
1844 : +du*du*(fivehalf*threehalf*a52*u**half &
1845 : +2._dp*a2+2._dp*3._dp*a3*u+3._dp*4._dp*a4*u**2 &
1846 : +4._dp*5._dp*a5*u**3+5._dp*6._dp*a6*u**4 &
1847 : +6._dp*7._dp*a7*u**5+7._dp*8._dp*a8*u**6) &
1848 129516 : +7._dp*aln*y**2+12._dp*aln*log(y)*y**2
1849 129516 : dden = dv*(b1+2._dp*b2*v+3._dp*b3*v**2+4._dp*b4*v**3)
1850 : d2den = d2v*(b1+2._dp*b2*v+3._dp*b3*v**2+4._dp*b4*v**3) &
1851 129516 : +dv*dv*(2._dp*b2+2._dp*3._dp*b3*v+3._dp*4._dp*b4*v**2)
1852 : ! derivatives w.r.t. y
1853 129516 : dfit = dnum/den - (num/den)*(dden/den)
1854 : d2fit = d2num/den - (dnum/den)*(dden/den) &
1855 : - (dnum/den)*(dden/den) + 2._dp*(num/den)*(dden/den)*(dden/den) &
1856 129516 : - (num/den)*(d2den/den)
1857 : ! Ax, and derivatives w.r.t. t
1858 129516 : Ax = fit
1859 129516 : dAx = dfit * dydt
1860 129516 : d2Ax = d2fit*dydt**2 + dfit*d2ydt2
1861 129516 : if(Ax/=Ax) Ax=zero
1862 129516 : if(dAx/=dAx) dAx=zero
1863 129516 : if(d2Ax/=d2Ax) d2Ax=zero
1864 :
1865 129516 : end subroutine tildeAx
1866 : !!***
1867 :
1868 : !!****f* m_special_funcs/tildeBx
1869 : !! NAME
1870 : !! tildeBx
1871 : !!
1872 : !! FUNCTION
1873 : !! Returns tilde Bx Pade fit and first and second derivatives
1874 : !! w.r.t. reduced temperature t.
1875 : !!
1876 : !! NOTES
1877 : !! Karasiev-Dufty-Trickey (KDT16) TGGA xc-functional
1878 : !! V.V. Karasiev, J.W. Dufty, and S.B. Trickey, PRL 120(7), 076401 (2018) [[cite:Karasiev2018]]
1879 : !!
1880 : !! INPUTS
1881 : !! t=reduced temperature, t=T/T_F
1882 : !!
1883 : !! OUTPUT
1884 : !! Bx=tilde Bx(t)
1885 : !! dBx=dBx(t)/dt
1886 : !! d2Bx=d^2Bx(t)/dt^2
1887 : !!
1888 : !! SOURCE
1889 :
1890 64758 : subroutine tildeBx(t,Bx,dBx,d2Bx)
1891 :
1892 : !Arguments ------------------------------------
1893 : !scalars
1894 : real(dp),intent(in) :: t
1895 : real(dp),intent(out) :: Bx,dBx,d2Bx
1896 :
1897 : !Local variables ------------------------------
1898 : !scalars
1899 : real(dp),parameter :: a2 = -3.4341427276599950_dp
1900 : real(dp),parameter :: a3 = -0.9066069544311700_dp
1901 : real(dp),parameter :: a4 = 2.2386316137237001_dp
1902 : real(dp),parameter :: a5 = 2.4232553178542000_dp
1903 : real(dp),parameter :: a6 = -0.1339278564306200_dp
1904 : real(dp),parameter :: a7 = 0.4392739633708200_dp
1905 : real(dp),parameter :: a8 = -0.0497109675177910_dp
1906 : real(dp),parameter :: a9 = 0.0000000000000000_dp
1907 : real(dp),parameter :: a10= 0.0028609701106953_dp
1908 : real(dp),parameter :: b1 = 0.7098198258073800_dp
1909 : real(dp),parameter :: b2 = 4.6311326377185997_dp
1910 : real(dp),parameter :: b3 = -2.9243190977647000_dp
1911 : real(dp),parameter :: b4 = 6.1688157841895004_dp
1912 : real(dp),parameter :: b5 = -1.3435764191535999_dp
1913 : real(dp),parameter :: b6 = 0.1576046383295400_dp
1914 : real(dp),parameter :: b7 = 0.4365792821186800_dp
1915 : real(dp),parameter :: b8 = -0.0620444574606262_dp
1916 : real(dp),parameter :: b9 = 0.0000000000000000_dp
1917 : real(dp),parameter :: b10= 0.0028609701106953_dp
1918 : real(dp),parameter :: onethird = 1._dp/3._dp
1919 : real(dp),parameter :: twothird = 2._dp/3._dp
1920 : real(dp),parameter :: threehalf = 3._dp/2._dp
1921 : real(dp),parameter :: fivehalf = 5._dp/2._dp
1922 : real(dp),parameter :: sevenhalf = 7._dp/2._dp
1923 : real(dp) :: y,u,du,d2u
1924 : real(dp) :: v,dv,d2v
1925 : real(dp) :: dydt,d2ydt2
1926 : real(dp) :: num,den,fit,dnum,d2num,dden,d2den,dfit,d2fit
1927 : ! *************************************************************************
1928 :
1929 64758 : y = twothird/t**threehalf
1930 64758 : u = y**twothird
1931 64758 : du = twothird/y**onethird
1932 64758 : d2u = -onethird*du/y
1933 :
1934 64758 : v = u
1935 64758 : dv = du
1936 64758 : d2v = d2u
1937 :
1938 64758 : dydt = -1._dp/t**fivehalf
1939 64758 : d2ydt2 = fivehalf/t**sevenhalf
1940 :
1941 : num = a2*u**2+a3*u**3+a4*u**4 &
1942 : +a5*u**5+a6*u**6+a7*u**7+a8*u**8 &
1943 64758 : +a9*u**9+a10*u**10
1944 : den = 1._dp+b1*v+b2*v**2+b3*v**3+b4*v**4 &
1945 64758 : +b5*v**5+b6*v**6+b7*v**7+b8*v**8+b9*v**9+b10*v**10
1946 64758 : fit = num/den
1947 :
1948 : dnum = du*(2._dp*a2*u+3._dp*a3*u**2+4._dp*a4*u**3 &
1949 : +5._dp*a5*u**4+6._dp*a6*u**5+7._dp*a7*u**6+8._dp*a8*u**7 &
1950 64758 : +9._dp*a9*u**8+10._dp*a10*u**9)
1951 :
1952 : d2num = d2u*(2._dp*a2*u+3._dp*a3*u**2+4._dp*a4*u**3 &
1953 : +5._dp*a5*u**4+6._dp*a6*u**5+7._dp*a7*u**6+8._dp*a8*u**7 &
1954 : +9._dp*a9*u**8+10._dp*a10*u**9) &
1955 : +du*du*(2._dp*a2+2._dp*3._dp*a3*u+3._dp*4._dp*a4*u**2 &
1956 : +4._dp*5._dp*a5*u**3+5._dp*6._dp*a6*u**4 &
1957 : +6._dp*7._dp*a7*u**5+7._dp*8._dp*a8*u**6 &
1958 64758 : +8._dp*9._dp*a9*u**7+9._dp*10._dp*a10*u**8)
1959 :
1960 : dden = dv*(b1+2._dp*b2*v+3._dp*b3*v**2+4._dp*b4*v**3 &
1961 : +5._dp*b5*v**4+6._dp*b6*v**5+7._dp*b7*v**6+8._dp*b8*v**7 &
1962 64758 : +9._dp*b9*v**8+10._dp*b10*v**9)
1963 :
1964 : d2den = d2v*(b1+2._dp*b2*v+3._dp*b3*v**2+4._dp*b4*v**3 &
1965 : +5._dp*b5*v**4+6._dp*b6*v**5+7._dp*b7*v**6+8._dp*b8*v**7 &
1966 : +9._dp*b9*v**8+10._dp*b10*v**9) &
1967 : + dv*dv*(2._dp*b2+2._dp*3._dp*b3*v+3._dp*4._dp*b4*v**2 &
1968 : +4._dp*5._dp*b5*v**3+5._dp*6._dp*b6*v**4+6._dp*7._dp*b7*v**5+7._dp*8._dp*b8*v**6 &
1969 64758 : +8._dp*9._dp*b9*v**7+9._dp*10._dp*b10*v**8)
1970 :
1971 : ! derivatives w.r.t. y
1972 64758 : dfit = dnum/den - (num/den)*(dden/den)
1973 : d2fit = d2num/den - (dnum/den)*(dden/den) &
1974 : - (dnum/den)*(dden/den) + 2._dp*(num/den)*(dden/den)*(dden/den) &
1975 64758 : - (num/den)*(d2den/den)
1976 :
1977 : ! Bx, and derivatives w.r.t. t
1978 64758 : Bx = fit
1979 64758 : dBx = dfit * dydt
1980 64758 : d2Bx = d2fit*dydt**2 + dfit*d2ydt2
1981 64758 : if(Bx/=Bx) Bx=zero
1982 64758 : if(dBx/=dBx) dBx=zero
1983 64758 : if(d2Bx/=d2Bx) d2Bx=zero
1984 :
1985 64758 : end subroutine tildeBx
1986 : !!***
1987 :
1988 : !!****f* m_special_funcs/tildeBc
1989 : !! NAME
1990 : !! tildeBc
1991 : !!
1992 : !! FUNCTION
1993 : !! Returns tilde Bc Pade fit and first and second derivatives
1994 : !! w.r.t. reduced temperature t.
1995 : !!
1996 : !! NOTES
1997 : !! Karasiev-Dufty-Trickey (KDT16) TGGA xc-functional
1998 : !! V.V. Karasiev, J.W. Dufty, and S.B. Trickey, PRL 120(7), 076401 (2018) [[cite:Karasiev2018]]
1999 : !!
2000 : !! INPUTS
2001 : !! iflag=flag selector integer
2002 : !! rs=Wigner-Seitz radius (bohr)
2003 : !! t=reduced temperature, t=T/T_F
2004 : !!
2005 : !! OUTPUT
2006 : !! Bc=tilde Bc(t)
2007 : !! dBcdrs=dBxc(rs,t)/drs
2008 : !! dBcdt=dBc(rs,t)/dt
2009 : !!
2010 : !! SOURCE
2011 :
2012 64758 : subroutine tildeBc(iflag,rs,t,Bc,dBcdrs,dBcdt)
2013 :
2014 : !Arguments ------------------------------------
2015 : !scalars
2016 : integer,intent(in) :: iflag
2017 : real(dp),intent(in) :: rs,t
2018 : real(dp),intent(out) :: Bc,dBcdrs,dBcdt
2019 :
2020 : !Local variables ------------------------------
2021 : !scalars
2022 : real(dp),parameter :: alpha_n = 0.50000000000000D+00
2023 : real(dp),parameter :: alpha_d = 0.15000000000000D+01
2024 : real(dp),parameter :: alpha_t = 0.32500000000000D+01
2025 : real(dp),parameter :: a1 = 0.30047772904141D+03
2026 : real(dp),parameter :: b1 = -0.11166043894641D+03
2027 : real(dp),parameter :: a2 = -0.38706401119284D+03
2028 : real(dp),parameter :: b2 = -0.45327974938936D+02
2029 : real(dp),parameter :: a3 = 0.25112236519758D+04
2030 : real(dp),parameter :: b3 = -0.14507109325068D+04
2031 : real(dp),parameter :: a4 = 0.52243427453456D+03
2032 : real(dp),parameter :: b4 = -0.30665095324907D+02
2033 : real(dp),parameter :: c1 = 0.11077393333429D+03
2034 : real(dp),parameter :: d1 = 0.12854960224127D+01
2035 : real(dp),parameter :: c2 = 0.32355494275181D+03
2036 : real(dp),parameter :: d2 = 0.13482659120012D+02
2037 : real(dp),parameter :: c3 = 0.45509212104516D+03
2038 : real(dp),parameter :: d3 = 0.23416017878226D+02
2039 : real(dp),parameter :: c4 = 0.10884351801356D+04
2040 : real(dp),parameter :: d4 = 0.24480831491950D+02
2041 : real(dp),parameter :: c5 = 0.36112604933128D+00
2042 : real(dp),parameter :: d5 = 0.32161372287131D-08
2043 : real(dp),parameter :: e1 = 0.32175261286726D+02
2044 : real(dp),parameter :: e2 = 0.61853047558212D+02
2045 : real(dp),parameter :: e3 = 0.33585054134674D+03
2046 : real(dp),parameter :: e4 = 0.12874240529185D+03
2047 : real(dp),parameter :: f1 = 0.41006056761680D-02
2048 : real(dp),parameter :: f2 = 0.18933118065366D-01
2049 : real(dp),parameter :: f3 = 0.24295412676204D-04
2050 : real(dp),parameter :: f4 = 0.18369775992299D-07
2051 : real(dp),parameter :: f5 = 0.69274680951701D-10
2052 : real(dp),parameter :: onethird = 1.d0/3.d0
2053 : real(dp),parameter :: twothird = 2.d0/3.d0
2054 : real(dp),parameter :: threehalf = 3.d0/2.d0
2055 : real(dp),parameter :: fivehalf = 5.d0/2.d0
2056 : real(dp),parameter :: sevenhalf = 7.d0/2.d0
2057 : real(dp) :: rsn,rsd,u,du
2058 : real(dp) :: num,den,dnumdrs,dnumdt,ddendrs,ddendt
2059 : ! *************************************************************************
2060 :
2061 64758 : if(iflag==5.or.iflag==6.or.iflag==7.or.iflag==8) then
2062 : ! Bc(rs,t) = 1
2063 0 : Bc = 1._dp
2064 0 : dBcdrs = 0._dp
2065 0 : dBcdt = 0._dp
2066 64758 : elseif(iflag==1.or.iflag==2.or.iflag==3.or.iflag==4) then
2067 : !
2068 : ! Bc(rs,t) = Pade Fit
2069 : !
2070 64758 : rsn = rs**alpha_n
2071 64758 : rsd = rs**alpha_d
2072 64758 : u = t**alpha_t
2073 64758 : du = alpha_t*t**(alpha_t-1.d0)
2074 : !
2075 64758 : num = 1.d0+(a1+b1*rsn+e1*rsn**2)*u+(a2+b2*rsn+e2*rsn**2)*u**2+(a3+b3*rsn+e3*rsn**2)*u**3+(a4+b4*rsn+e4*rsn**2)*u**4
2076 64758 : dnumdrs = (b1+2.d0*e1*rsn)*u+(b2+2.d0*e2*rsn)*u**2+(b3+2.d0*e3*rsn)*u**3+(b4+2.d0*e4*rsn)*u**4
2077 64758 : dnumdrs = dnumdrs * alpha_n*rs**(alpha_n-1.d0)
2078 64758 : dnumdt = (a1+b1*rsn+e1*rsn**2)+2.d0*(a2+b2*rsn+e2*rsn**2)*u+3.d0*(a3+b3*rsn+e3*rsn**2)*u**2+4.d0*(a4+b4*rsn+e4*rsn**2)*u**3
2079 64758 : dnumdt = dnumdt * du
2080 : !
2081 64758 : den = 1.d0+(c1+d1*rsd+f1*rsd**2)*u+(c2+d2*rsd+f2*rsd**2)*u**2+(c3+d3*rsd+f3*rsd**2)*u**3+(c4+d4*rsd+f4*rsd**2)*u**4+(c5+d5*rsd+f5*rsd**2)*u**5
2082 64758 : ddendrs = (d1+2.d0*f1*rsd)*u+(d2+2.d0*f2*rsd)*u**2+(d3+2.d0*f3*rsd)*u**3+(d4+2.d0*f4*rsd)*u**4+(d5+2.d0*f5*rsd)*u**5
2083 64758 : ddendrs = ddendrs * alpha_d*rs**(alpha_d-1.d0)
2084 64758 : ddendt = (c1+d1*rsd+f1*rsd**2)+2.d0*(c2+d2*rsd+f2*rsd**2)*u+3.d0*(c3+d3*rsd+f3*rsd**2)*u**2+4.d0*(c4+d4*rsd+f4*rsd**2)*u**3+5.d0*(c5+d5*rsd+f5*rsd**2)*u**4
2085 64758 : ddendt = ddendt * du
2086 : !
2087 64758 : Bc = num/den
2088 64758 : dBcdrs = dnumdrs/den - (num/den)*ddendrs/den
2089 64758 : dBcdt = dnumdt/den - (num/den)*ddendt/den
2090 64758 : if(Bc/=Bc) Bc=zero
2091 64758 : if(dBcdrs/=dBcdrs) dBcdrs=zero
2092 64758 : if(dBcdt/=dBcdt) dBcdt=zero
2093 : !
2094 : endif
2095 :
2096 64758 : end subroutine tildeBc
2097 : !!***
2098 :
2099 : !----------------------------------------------------------------------
2100 :
2101 : !!****f* m_special_funcs/k_fermi
2102 : !! NAME
2103 : !! k_fermi
2104 : !!
2105 : !! FUNCTION
2106 : !! Returns the Fermi wave vector corresponding to the local value of the real space density rhor.
2107 : !!
2108 : !! INPUTS
2109 : !! rhor=Local density in real space.
2110 : !!
2111 : !! SOURCE
2112 :
2113 280434564 : elemental function k_fermi(rhor)
2114 :
2115 : !Arguments ------------------------------------
2116 : !scalars
2117 : real(dp),intent(in) :: rhor
2118 : real(dp) :: k_fermi
2119 :
2120 : !Local variables-------------------------------
2121 : !scalars
2122 : real(dp),parameter :: pisq=pi**2
2123 : ! *************************************************************************
2124 :
2125 280434564 : k_fermi = (three*pisq*rhor)**third
2126 :
2127 140216832 : end function k_fermi
2128 : !!***
2129 :
2130 : !----------------------------------------------------------------------
2131 :
2132 : !!****f* m_special_funcs/k_thfermi
2133 : !! NAME
2134 : !! k_thfermi
2135 : !!
2136 : !! FUNCTION
2137 : !! Returns the Thomas-Fermi wave vector corresponding to the local value of the real space density rhor.
2138 : !!
2139 : !! INPUTS
2140 : !! rhor=Local density in real space.
2141 : !!
2142 : !! SOURCE
2143 :
2144 140217732 : elemental function k_thfermi(rhor)
2145 :
2146 : !Arguments ------------------------------------
2147 : !scalars
2148 : real(dp),intent(in) :: rhor
2149 : real(dp) :: k_thfermi
2150 :
2151 : !Local variables-------------------------------
2152 : !scalars
2153 : real(dp),parameter :: pisq=pi**2
2154 : ! *************************************************************************
2155 :
2156 140217732 : k_thfermi = SQRT(four*k_fermi(rhor)*piinv)
2157 :
2158 140217732 : end function k_thfermi
2159 : !!***
2160 :
2161 : !----------------------------------------------------------------------
2162 :
2163 : !!****f* m_special_funcs/levi_civita_3
2164 : !! NAME
2165 : !! levi_civita_3
2166 : !!
2167 : !! FUNCTION
2168 : !! Return Levi-Civita tensor of rank 3
2169 : !!
2170 : !! SOURCE
2171 :
2172 0 : pure function levi_civita_3() result(ee)
2173 :
2174 : !Arguments ------------------------------------
2175 : integer :: ee(3,3,3)
2176 : ! *************************************************************************
2177 :
2178 0 : ee = 0
2179 0 : ee(1,2,3) = 1
2180 0 : ee(2,3,1) = 1
2181 0 : ee(3,1,2) = 1
2182 : !
2183 0 : ee(3,2,1) = -1
2184 0 : ee(1,3,2) = -1
2185 0 : ee(2,1,3) = -1
2186 :
2187 0 : end function levi_civita_3
2188 : !!***
2189 :
2190 : !!****f* m_special_funcs/jlspline_init
2191 : !! NAME
2192 : !! jlspline_init
2193 : !!
2194 : !! FUNCTION
2195 : !! Pre-calculate the j_v(y) for recip_ylm on regular grid
2196 : !! NOTE: spherical Bessel function small j!
2197 : !!
2198 : !! INPUTS
2199 : !! nx = max number of points on grid for integral
2200 : !! delta = space between integral arguments
2201 : !! mlang= max angular momentum
2202 : !!
2203 : !! OUTPUT
2204 : !! bess_spl=array of integrals
2205 : !! bess_spl_der=array of derivatives of integrals
2206 : !! xx=coordinates of points belonging to the grid
2207 : !!
2208 : !! SOURCE
2209 :
2210 81 : subroutine jlspline_init(new, nx, delta, mlang)
2211 :
2212 : !Arguments ------------------------------------
2213 : !scalars
2214 : class(jlspline_t),intent(inout) :: new
2215 : integer,intent(in) :: nx,mlang
2216 : real(dp),intent(in) :: delta
2217 :
2218 : !Local variables -------------------------
2219 : !scalars
2220 : integer :: ix,ll
2221 : real(dp) :: yp1,ypn
2222 : !arrays
2223 81 : real(dp),allocatable :: cosbessx(:),sinbessx(:)
2224 : ! *********************************************************************
2225 :
2226 81 : if (nx < 2) then
2227 0 : ABI_ERROR('need more than one point for the interpolation routines')
2228 : end if
2229 :
2230 81 : new%nx = nx; new%mlang = mlang; new%delta = delta; new%maxarg = (nx-1) * delta
2231 243 : ABI_MALLOC(new%xx, (nx))
2232 324 : ABI_MALLOC(new%bess_spl, (nx, mlang))
2233 243 : ABI_MALLOC(new%bess_spl_der, (nx, mlang))
2234 :
2235 : !-----------------------------------------------------------------
2236 : !Bessel function into array
2237 : !-----------------------------------------------------------------
2238 : ! integration grid is nfiner times finer than the interpolation grid
2239 162 : ABI_MALLOC(sinbessx, (nx))
2240 162 : ABI_MALLOC(cosbessx, (nx))
2241 :
2242 : ! could be done by chain rule for cos sin (is it worth it?) but
2243 : ! precision problems as numerical errors are propagated.
2244 41381 : do ix=1,nx
2245 41300 : new%xx(ix) = (ix-1) * delta
2246 41300 : sinbessx(ix) = sin(new%xx(ix))
2247 41381 : cosbessx(ix) = cos(new%xx(ix))
2248 : end do
2249 :
2250 : ! fill bess_spl array
2251 486 : do ll=0,mlang-1
2252 405 : call besjm(one,new%bess_spl(:,ll+1),cosbessx,ll,nx,sinbessx,new%xx)
2253 :
2254 : ! call spline to get 2nd derivative (reuse in splint later)
2255 405 : yp1 = zero; ypn = zero
2256 486 : call spline(new%xx, new%bess_spl(:,ll+1), nx, yp1, ypn, new%bess_spl_der(:,ll+1))
2257 : end do
2258 :
2259 : !write(std_out,*) ' bess funct 0 1 2 3 4'
2260 : !do ix=1,nx
2261 : !write(std_out,*) xx(ix), (new%bess_spl(ix,ll),ll=1,mlang)
2262 : !end do
2263 :
2264 81 : ABI_FREE(sinbessx)
2265 81 : ABI_FREE(cosbessx)
2266 :
2267 81 : end subroutine jlspline_init
2268 : !!***
2269 :
2270 : !----------------------------------------------------------------------
2271 :
2272 : !!****f* m_special_funcs/jlspline_free
2273 : !! NAME
2274 : !! jlspline_free
2275 : !!
2276 : !! FUNCTION
2277 : !! deallocate memory
2278 : !!
2279 : !! SOURCE
2280 :
2281 81 : subroutine jlspline_free(jlspl)
2282 :
2283 : !Arguments ------------------------------------
2284 : class(jlspline_t),intent(inout) :: jlspl
2285 : ! *********************************************************************
2286 :
2287 81 : ABI_SFREE(jlspl%xx)
2288 81 : ABI_SFREE(jlspl%bess_spl)
2289 81 : ABI_SFREE(jlspl%bess_spl_der)
2290 :
2291 81 : end subroutine jlspline_free
2292 : !!***
2293 :
2294 : !----------------------------------------------------------------------
2295 :
2296 : !!****f* m_special_funcs/jlspline_integral
2297 : !! NAME
2298 : !! jlspline_integral
2299 : !!
2300 : !! INPUTS
2301 : !!
2302 : !! OUTPUT
2303 : !!
2304 : !! FUNCTION
2305 : !!
2306 : !! SOURCE
2307 :
2308 0 : real(dp) function jlspline_integral(jlspl, il, qq, powr, nr, rcut) result(res)
2309 :
2310 : !Arguments ------------------------------------
2311 : class(jlspline_t),intent(in) :: jlspl
2312 : integer,intent(in) :: il,nr,powr
2313 : real(dp),intent(in) :: qq, rcut
2314 :
2315 : !Local variables ---------------------------------------
2316 : integer :: ierr
2317 : real(dp) :: step
2318 : !arrays
2319 0 : real(dp):: xfit(nr),yfit(nr),rr(nr)
2320 : ! *********************************************************************
2321 :
2322 0 : step = rcut / (nr - 1)
2323 0 : rr = arth(zero, step, nr)
2324 0 : xfit = qq * rr
2325 0 : call splint(jlspl%nx, jlspl%xx, jlspl%bess_spl(:,il), jlspl%bess_spl_der(:,il), nr, xfit, yfit, ierr=ierr)
2326 :
2327 0 : if (ierr /= 0) then
2328 0 : write(std_out,*)"qq, rcut, qq*rcut, maxarg", qq, rcut, qq*rcut, jlspl%maxarg
2329 0 : write(std_out,*)"x[0], x[-1]",jlspl%xx(1),jlspl%xx(jlspl%nx)
2330 0 : write(std_out,*)"minval xfit: ",minval(xfit)
2331 0 : write(std_out,*)"maxval xfit: ",maxval(xfit)
2332 0 : ABI_ERROR("splint returned ierr != 0")
2333 : end if
2334 :
2335 0 : if (powr /= 1) yfit = yfit * (rr ** powr)
2336 0 : res = simpson(step, yfit)
2337 :
2338 0 : end function jlspline_integral
2339 : !!***
2340 :
2341 : !!****f* m_special_funcs/gspline_init
2342 : !! NAME
2343 : !! gspline_init
2344 : !!
2345 : !! FUNCTION
2346 : !! Build object to spline the gaussian approximant and its primitive.
2347 : !!
2348 : !! INPUTS
2349 : !! sigma=Broadening parameter.
2350 : !!
2351 : !! SOURCE
2352 :
2353 0 : subroutine gspline_init(new, sigma)
2354 :
2355 : !Arguments ------------------------------------
2356 : !scalars
2357 : class(gspline_t),intent(out) :: new
2358 : real(dp),intent(in) :: sigma
2359 :
2360 : !Local variables ------------------------------
2361 : integer :: ii
2362 : real(dp) :: ybcbeg, ybcend
2363 : ! *************************************************************************
2364 :
2365 0 : new%nspline = 5 * 1024; new%sigma = sigma
2366 0 : ABI_CHECK(sigma > zero, sjoin("invalid sigma:", ftoa(sigma)))
2367 0 : new%xmin = zero
2368 0 : new%xmax = sigma * sqrt(-log(sigma * sqrt(pi) * tol12)) ! gauss(xmax) = tol12
2369 0 : new%step = (new%xmax - new%xmin) / (new%nspline - 1)
2370 0 : new%stepm1 = one / new%step; new%step2div6 = new%step**2 / six
2371 :
2372 0 : ABI_MALLOC(new%xvals, (new%nspline))
2373 0 : do ii=1,new%nspline
2374 0 : new%xvals(ii) = new%xmin + (ii-1) * new%step
2375 : end do
2376 0 : new%xmax = new%xvals(new%nspline)
2377 :
2378 : ! Spline the gaussian approximant.
2379 0 : ABI_MALLOC(new%svals, (new%nspline, 4))
2380 0 : new%svals(:, 1) = gaussian(new%xvals, sigma)
2381 0 : ybcbeg = - (two * new%xmin / sigma**2) * new%svals(1,1)
2382 0 : ybcend = - (two * new%xmax / sigma**2) * new%svals(new%nspline,1)
2383 0 : call spline(new%xvals, new%svals(:,1), new%nspline, ybcbeg, ybcend, new%svals(:,2))
2384 :
2385 : ! Spline the primitive: 1/2 [1 + erf(x/sigma)]
2386 0 : new%svals(:, 3) = half * (one + abi_derf(new%xvals / new%sigma))
2387 0 : call spline(new%xvals, new%svals(:,3), new%nspline, new%svals(1,1), new%svals(new%nspline, 1), new%svals(:,4))
2388 : !do ii=1,new%nspline; write(98,*)new%xvals(ii),new%svals(ii,3),new%svals(ii,4); end do
2389 :
2390 0 : end subroutine gspline_init
2391 : !!***
2392 :
2393 : !!****f* m_special_funcs/gspline_eval
2394 : !! NAME
2395 : !! gspline_eval
2396 : !!
2397 : !! FUNCTION
2398 : !! Evaluate the gaussian approximant and its primitive at (xmesh - x0)
2399 : !!
2400 : !! INPUTS
2401 : !! self<gspline_t>=Object used to spline the gaussian approximant
2402 : !! x0=Shift to be given to xmesh
2403 : !! nx=Number of points in input mesh.
2404 : !! xmesh(nx)=Frequency points (not necessary linear).
2405 : !!
2406 : !! OUTPUT
2407 : !! weights(nx,2)=First slice contains the gaussian approximant on xmesh.
2408 : !! The second slice stores the primitive.
2409 : !!
2410 : !! SOURCE
2411 :
2412 0 : pure subroutine gspline_eval(self, x0, nx, xmesh, weights)
2413 :
2414 : !Arguments ------------------------------------
2415 : !scalars
2416 : class(gspline_t),intent(in) :: self
2417 : integer,intent(in) :: nx
2418 : real(dp),intent(in) :: x0
2419 : !arrays
2420 : real(dp),intent(in) :: xmesh(nx)
2421 : real(dp),intent(out) :: weights(nx,2)
2422 :
2423 : !Local variables ------------------------------
2424 : !scalars
2425 : integer :: ix,jspl
2426 : real(dp) :: xx,absx,aa,bb,cc,dd
2427 : logical :: isneg
2428 : !real(dp) :: int_values(nx)
2429 : ! *************************************************************************
2430 :
2431 0 : do ix=1,nx
2432 0 : xx = xmesh(ix) - x0; absx = abs(xx); isneg = xx < zero
2433 0 : if (absx >= self%xmax) then
2434 : ! Region in which gauss(x) is negligible.
2435 0 : weights(ix,1) = zero
2436 0 : if (isneg) then
2437 0 : weights(ix,2) = zero
2438 : else
2439 0 : weights(ix,2) = one
2440 : end if
2441 : else
2442 : ! Spline functions at |x| and recover the value at x:
2443 : ! g(x) = g(-x); G(-x) = 1 - G(x)
2444 0 : jspl = 1 + int((absx - self%xmin) * self%stepm1); dd = absx - self%xvals(jspl)
2445 0 : bb = dd * self%stepm1
2446 0 : aa = one - bb
2447 0 : cc = aa*(aa**2-one) * self%step2div6
2448 0 : dd = bb*(bb**2-one) * self%step2div6
2449 :
2450 0 : weights(ix,1) = aa*self%svals(jspl,1) + bb*self%svals(jspl+1,1) + cc*self%svals(jspl,2) + dd*self%svals(jspl+1,2)
2451 0 : weights(ix,2) = aa*self%svals(jspl,3) + bb*self%svals(jspl+1,3) + cc*self%svals(jspl,4) + dd*self%svals(jspl+1,4)
2452 0 : if (isneg) weights(ix,2) = one - weights(ix,2)
2453 : end if
2454 : end do
2455 :
2456 : !call simpson_int(nx,xmesh(2) - xmesh(1),weights(:,1),int_values)
2457 : !do ix=1,nx
2458 : ! write(99,*)xmesh(ix), weights(ix,1), gaussian(xx, self%sigma), weights(ix,2), int_values(ix)
2459 : !end do
2460 :
2461 0 : end subroutine gspline_eval
2462 : !!***
2463 :
2464 : !!****f* m_special_funcs/gspline_free
2465 : !! NAME
2466 : !! gspline_free
2467 : !!
2468 : !! FUNCTION
2469 : !! Free dynamic memory
2470 : !!
2471 : !! SOURCE
2472 :
2473 0 : subroutine gspline_free(self)
2474 :
2475 : !Arguments ------------------------------------
2476 : class(gspline_t),intent(inout) :: self
2477 : ! *************************************************************************
2478 :
2479 0 : ABI_SFREE(self%xvals)
2480 0 : ABI_SFREE(self%svals)
2481 :
2482 0 : end subroutine gspline_free
2483 : !!***
2484 :
2485 0 : end module m_special_funcs
2486 : !!***
|
