Branch data 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 : 1288313644 : 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 : 1288313644 : xx = yy
770 : 1288313644 : isw = 1
771 : : !Here change the sign of xx, and keep track of it thanks to isw
772 [ + + ]: 1288313644 : if (xx<0.0e0_dp) then
773 : 159608 : isw = -1
774 : 159608 : xx = -xx
775 : : end if
776 : :
777 : 1288313644 : done=0
778 : :
779 : : !Residual value, if yy < -6.375e0_dp
780 : 1288313644 : res=2.0e0_dp
781 : :
782 : : !abs(yy) < .477, evaluate approximation for erfc
783 [ + + ]: 1288313644 : 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 [ + + + + ]: 1288313644 : if (xx<=4.0e0_dp .and. done==0 ) then
804 : 160134854 : xsq = xx*xx
805 : 160134854 : xnum = p1(8)*xx+p1(9)
806 : 160134854 : xden = xx+q1(8)
807 [ + + ]: 1281078832 : do ii=1,7
808 : 1120943978 : xnum = xnum*xx+p1(ii)
809 : 1281078832 : xden = xden*xx+q1(ii)
810 : : end do
811 : 160134854 : res = xnum/xden
812 : 160134854 : res = res* exp(-xsq)
813 [ + + ]: 160134854 : if (isw.eq.-1) res = 2.0e0_dp-res
814 : : done=1
815 : : end if
816 : :
817 : : !y > 13.3e0_dp
818 [ + + + + : 1288289352 : 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 [ + + + + : 1288313644 : if ( ( isw > 0 .or. xx < xlarge ) .and. done==0 ) then
+ + ]
826 : 1122560880 : xsq = xx*xx
827 : 1122560880 : xi = 1.0e0_dp/xsq
828 : 1122560880 : xnum= p2(5)*xi+p2(6)
829 : 1122560880 : xden = xi+q2(5)
830 [ + + ]: 5612804400 : do ii = 1,4
831 : 4490243520 : xnum = xnum*xi+p2(ii)
832 : 5612804400 : xden = xden*xi+q2(ii)
833 : : end do
834 : 1122560880 : res = (sqrpi+xi*xnum/xden)/xx
835 : 1122560880 : res = res* exp(-xsq)
836 [ + + ]: 1122560880 : if (isw.eq.-1) res = 2.0e0_dp-res
837 : : end if
838 : :
839 : : !All cases have been investigated
840 : 1288313644 : derfc_yy = res
841 : :
842 : 1288313644 : 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 : 4816696 : 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 [ + + + + : 4816696 : if (nn==0) then
+ - - ]
1043 : :
1044 : 2376587 : switchx=nx+1
1045 [ + + ]: 646744437 : do ix=1,nx
1046 : 646734810 : rr=arg*xx(ix)
1047 [ + + ]: 646744437 : if (rr<=1.d-1) then
1048 : 644367850 : rsq=rr*rr
1049 : 644367850 : 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 [ + + ]: 1324584649 : do ix=switchx,nx
1057 : 1322208062 : rr=arg*xx(ix)
1058 : 1324584649 : besjx(ix)=sinx(ix)/rr
1059 : : end do
1060 : :
1061 : : else if (nn==1) then
1062 : :
1063 : 1872757 : switchx=nx+1
1064 [ + + ]: 783899660 : do ix=1,nx
1065 : 783830243 : rr=arg*xx(ix)
1066 [ + + ]: 783899660 : if (rr<=1.d0) then
1067 : 782026903 : rsq=rr*rr
1068 : 782026903 : top=1.d0-rsq*(b11-rsq*(b12-rsq*(b13-rsq*b14)))
1069 : 782026903 : bot=1.d0+rsq*(c11+rsq*(c12+rsq*(c13+rsq*c14)))
1070 : 782026903 : besjx(ix)=third*rr*top/bot
1071 : : else
1072 : : switchx=ix
1073 : : exit
1074 : : end if
1075 : : end do
1076 : :
1077 [ + + ]: 784826876 : do ix=switchx,nx
1078 : 782954119 : rr=arg*xx(ix)
1079 : 782954119 : rsq=rr*rr
1080 : 784826876 : 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 : 4816696 : 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 : 815397275 : 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 : 815397275 : real(dp),allocatable :: sb(:)
1215 : : ! *************************************************************************
1216 : :
1217 [ + + ]: 815397275 : if(xx<= 1.0e-36_dp) then
1218 : : ! zero argument section
1219 [ + + ]: 7903793 : sb_out(:)=zero
1220 : 2080592 : sb_out(1)=one
1221 [ + + ]: 813316683 : else if(xx<1.e-3_dp) then
1222 : : ! small argument section
1223 : 15642588 : xn=one
1224 : 15642588 : xs=half*xx**2
1225 [ + + ]: 63750070 : do nn=1,nm
1226 : 48107482 : sb_out(nn)=xn*(one - xs*(one - xs/(4*nn+6))/(2*nn+1))
1227 : 63750070 : xn=xx*xn/(2*nn+1)
1228 : : end do
1229 : : else
1230 : : ! recursion method
1231 [ + + ]: 797674095 : if(xx<one) then
1232 : 366539106 : nlim=nm+int(15.0e0_dp*xx)+1
1233 : : else
1234 : 431134989 : nlim=nm+int(1.36e0_dp*xx)+15
1235 : : end if
1236 [ + - - + : 2393022285 : ABI_MALLOC(sb,(nlim+1))
+ - + - +
- - + -
+ ]
1237 : 797674095 : nn=nlim
1238 : 797674095 : xi=one/xx
1239 : 797674095 : sb(nn+1)=zero
1240 : 797674095 : sb(nn)=1.e-18_dp
1241 : 797674095 : sn=dble(2*nn-1)*1.e-36_dp
1242 [ + + ]: 14463745855 : do nn=nlim-1,1,-1
1243 : 14463745855 : sb(nn)=dble(2*nn+1)*xi*sb(nn+1) - sb(nn+2)
1244 : : end do
1245 [ + + ]: 14463745855 : do nn=1,nlim-1
1246 : 14463745855 : sn=sn + dble(2*nn-1)*sb(nn)*sb(nn)
1247 : : end do
1248 : 797674095 : fn=1.d0/sqrt(sn)
1249 [ + + ]: 3115456719 : sb_out(:)=fn*sb(1:nm)
1250 : 797674095 : ABI_FREE(sb)
1251 : : end if
1252 : :
1253 : 815397275 : 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 : : !!***
|
