Line data Source code
1 : !!****m* ABINIT/m_vtorhotf
2 : !! NAME
3 : !! m_vtorhotf
4 : !!
5 : !! FUNCTION
6 : !! Computes the new density from a fixed potential (vtrial) using the Thomas-Fermi functional
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, MF, AR, MM)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_vtorhotf
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_xmpi
28 : use m_dtset
29 :
30 : use defs_abitypes, only : MPI_type
31 : use m_time, only : timab
32 : use m_spacepar, only : symrhg
33 :
34 : implicit none
35 :
36 : private
37 : !!***
38 :
39 : public :: vtorhotf
40 : !!***
41 :
42 : contains
43 : !!***
44 :
45 : !!****f* ABINIT/vtorhotf
46 : !! NAME
47 : !! vtorhotf
48 : !!
49 : !! FUNCTION
50 : !! This routine computes the new density from a fixed potential (vtrial)
51 : !! using the Thomas-Fermi functional
52 : !!
53 : !! INPUTS
54 : !! dtset <type(dataset_type)>=all input variables for this dataset
55 : !! gprimd(3,3)=dimensional reciprocal space primitive translations
56 : !! irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
57 : !! mpi_enreg=information about MPI parallelization
58 : !! natom=number of atoms in cell.
59 : !! nfft=number of fft grid points
60 : !! nspden=number of spin-density components
61 : !! nsppol=1 for unpolarized, 2 for spin-polarized
62 : !! nsym=number of symmetry elements in space group
63 : !! phnons(2,nfft**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases
64 : !! gprimd(3,3)=dimensional real space primitive translations
65 : !! ucvol=unit cell volume in bohr**3.
66 : !! vtrial(nfft,nspden)=INPUT Vtrial(r).
67 : !!
68 : !! OUTPUT
69 : !! ek=kinetic energy part of total energy.
70 : !! enlx=nonlocal psp + potential Fock ACE part of total energy.
71 : !! entropy=entropy due to the occupation number smearing (if metal)
72 : !! fermie=fermi energy (Hartree)
73 : !! grnl(3*natom)=stores grads of nonlocal energy wrt length scales
74 : !! (3x3 tensor) and grads wrt atomic coordinates (3*natom)
75 : !!
76 : !! SIDE EFFECTS
77 : !! rhog(2,nfft)=array for Fourier transform of electron density
78 : !! rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
79 : !!
80 : !! SOURCE
81 :
82 277 : subroutine vtorhotf(dtset,ek,enlx,entropy,fermie,gprimd,grnl,&
83 277 : & irrzon,mpi_enreg,natom,nfft,nspden,nsppol,nsym,phnons,rhog,rhor,rprimd,ucvol,vtrial)
84 :
85 : !Arguments ------------------------------------
86 : !scalars
87 : integer,intent(in) :: natom,nfft,nspden,nsppol,nsym
88 : real(dp),intent(in) :: ucvol
89 : real(dp),intent(out) :: ek,enlx,entropy,fermie
90 : type(MPI_type),intent(in) :: mpi_enreg
91 : type(dataset_type),intent(in) :: dtset
92 : !arrays
93 : integer,intent(in) :: irrzon((dtset%ngfft(1)*dtset%ngfft(1)*dtset%ngfft(1))**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))
94 : real(dp),intent(in) :: gprimd(3,3)
95 : real(dp),intent(in) :: phnons(2,(dtset%ngfft(1)*dtset%ngfft(1)*dtset%ngfft(1))**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))
96 : real(dp),intent(in) :: rprimd(3,3),vtrial(nfft,nspden)
97 : real(dp),intent(inout) :: rhog(2,nfft),rhor(nfft,nspden)
98 : real(dp),intent(out) :: grnl(3*natom)
99 :
100 : !Local variables-------------------------------
101 : !scalars
102 : integer,parameter :: jdichomax=20,level=111
103 : integer :: i1,i2,i3,ierr,ifft,ii,ir,iscf,jdicho
104 : integer :: me_fft,n1,n2,n3,nfftot,nproc_fft,prtvol
105 : real(dp),save :: cktf,fermie_tol,nelect_mid
106 : real(dp) :: dnelect_mid_dx,dxrtnewt,eektemp,eektf,feektemp,feektf
107 : real(dp) :: rtnewt,sum_rhor_mid,sum_rhor_middx
108 : logical,save :: lfirst_time_tf=.true.
109 : logical :: lnewtonraphson
110 : character(len=500) :: message
111 : !arrays
112 : real(dp) :: tsec(2)
113 277 : real(dp),allocatable :: betamumoinsV(:),rhor_mid(:),rhor_middx(:)
114 :
115 : ! *************************************************************************
116 :
117 : !Keep track of total time spent in vtorho
118 277 : call timab(21,1,tsec)
119 :
120 : !Structured debugging if prtvol==-level
121 277 : prtvol=dtset%prtvol
122 277 : if(prtvol==-level)then
123 0 : write(message,'(80a,a,a)') ('=',ii=1,80),ch10,' vtorho : enter '
124 0 : call wrtout(std_out,message,'COLL')
125 : end if
126 :
127 277 : n1=dtset%ngfft(1) ; n2=dtset%ngfft(2) ; n3=dtset%ngfft(3)
128 277 : me_fft=dtset%ngfft(11) ; nproc_fft=dtset%ngfft(10)
129 277 : iscf=dtset%iscf
130 : !Debugging : print vtrial and rhor
131 277 : if(prtvol==-level)then
132 0 : write(message,'(a)') ' ir vtrial(ir) rhor(ir) '
133 0 : call wrtout(std_out,message,'COLL')
134 0 : do ir=1,nfft
135 : ! if(ir<=11 .or. mod(ir,301)==0 )then
136 0 : i3=(ir-1)/n1/(n2/nproc_fft)
137 0 : i2=(ir-1-i3*n1*n2/nproc_fft)/n1
138 0 : i1=ir-1-i3*n1*n2/nproc_fft-(i2-me_fft)*n1
139 0 : write(message,'(i5,3i3,a,2d13.6)')ir,i1,i2,i3,' ',vtrial(ir,1),rhor(ir,1)
140 0 : call wrtout(std_out,message,'COLL')
141 0 : if(nspden>=2)then
142 0 : write(message,'(a,2d13.6)')' ',vtrial(ir,2),rhor(ir,2)
143 0 : call wrtout(std_out,message,'COLL')
144 : end if
145 : ! end if
146 : end do
147 : end if
148 :
149 277 : ek=zero
150 277 : enlx=zero
151 26869 : grnl(:)=zero
152 :
153 : !Initialize rhor if needed
154 12255962 : if(iscf>0) rhor(:,:)=zero
155 :
156 : !call Thomas-Fermi for the density
157 277 : call tf
158 : !Compute energy terms
159 277 : call tfek
160 :
161 277 : call timab(21,2,tsec)
162 : !End thomas fermi
163 :
164 : contains
165 : !!***
166 :
167 : !!****f* vtorhotf/tf
168 : !! NAME
169 : !! tf
170 : !!
171 : !! FUNCTION
172 : !!
173 : !! INPUTS
174 : !!
175 : !! OUTPUT
176 : !!
177 : !! SOURCE
178 277 : subroutine tf()
179 :
180 : ! *************************************************************************
181 :
182 831 : ABI_MALLOC(rhor_mid,(nfft))
183 831 : ABI_MALLOC(rhor_middx,(nfft))
184 277 : fermie_tol=1.e-10_dp
185 277 : cktf=one/two/pi**2*(two*dtset%tphysel)**1.5_dp
186 :
187 : ! Should be made an input variable, if TF really needed for production
188 : ! rtnewt=dtset%userra
189 277 : rtnewt=zero
190 :
191 : ! Newton Raphson
192 277 : if (lfirst_time_tf) then
193 4 : lfirst_time_tf=.false.
194 : end if
195 277 : jdicho=0
196 277 : lnewtonraphson=.false.
197 4641 : do while (.not.lnewtonraphson)
198 2182 : jdicho=jdicho+1
199 : ! do ifft=1,nfft
200 : ! rhor_mid(ifft)=cktf*zfermi12((rtnewt-vtrial(ifft,1))/dtset%tphysel)
201 : ! rhor_middx(ifft)=cktf*zfermim12((rtnewt-vtrial(ifft,1))/dtset%tphysel)
202 : ! end do
203 : call fm12a1t(cktf,rtnewt,dtset%tphysel,vtrial(:,1),rhor_middx,rhor_mid,&
204 2182 : & nfft)
205 97127446 : sum_rhor_mid=sum(rhor_mid(:))
206 97127446 : sum_rhor_middx=sum(rhor_middx(:))
207 : call xmpi_sum(sum_rhor_mid,mpi_enreg%comm_fft ,ierr)
208 : call xmpi_sum(sum_rhor_middx,mpi_enreg%comm_fft ,ierr)
209 2182 : nelect_mid=sum_rhor_mid*ucvol/(nfft*nproc_fft)-dtset%nelect
210 2182 : dnelect_mid_dx=sum_rhor_middx*ucvol/(nfft*nproc_fft)/dtset%tphysel/two
211 2182 : dxrtnewt=nelect_mid/dnelect_mid_dx
212 2182 : rtnewt=rtnewt-dxrtnewt
213 2182 : if (abs(nelect_mid) < fermie_tol/2._dp) then
214 277 : lnewtonraphson=.true.
215 : end if
216 2459 : if (jdicho > jdichomax) then
217 0 : ABI_ERROR('NEWTON RAPHSON NOT CONVERGED')
218 : end if
219 : end do
220 277 : fermie=rtnewt
221 12255685 : rhor(:,1)=rhor_mid(:)
222 277 : ABI_FREE(rhor_mid)
223 277 : ABI_FREE(rhor_middx)
224 :
225 : ! DEBUG
226 : ! write(std_out,*)'fmid,nmid,jdicho',fermie,nelect_mid,jdicho
227 : ! ENDDEBUG
228 :
229 : ! Compute rhog
230 277 : call timab(70,1,tsec)
231 :
232 277 : nfftot=dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3)
233 : call symrhg(1,gprimd,irrzon,mpi_enreg,nfft,nfftot,dtset%ngfft,nspden,nsppol,nsym,phnons,&
234 277 : & rhog,rhor,rprimd,dtset%symafm,dtset%symrel,dtset%tnons)
235 :
236 : ! We now have both rho(r) and rho(G), symmetrized, and if nsppol=2
237 : ! we also have the spin-up density, symmetrized, in rhor(:,2).
238 :
239 277 : call timab(70,2,tsec)
240 :
241 277 : end subroutine tf
242 : !!***
243 :
244 : !!****f* vtorhotf/tfek
245 : !! NAME
246 : !! tfek
247 : !!
248 : !! FUNCTION
249 : !! This is the calculation of the kinetic energy for Thomas Fermi
250 : !! Energy and free energy must be distinguished
251 : !!
252 : !! INPUTS
253 : !!
254 : !! OUTPUT
255 : !!
256 : !! SOURCE
257 :
258 277 : subroutine tfek()
259 :
260 : ! *************************************************************************
261 :
262 831 : ABI_MALLOC(betamumoinsV,(nfft))
263 277 : cktf=one/two/pi**2*(two*dtset%tphysel)**1.5_dp
264 277 : eektf=zero
265 277 : feektf=zero
266 12255685 : do ifft=1,nfft
267 :
268 : ! betamumoinsv(ifft)=ifermi12(rhor(ifft,1)/cktf)
269 12255408 : betamumoinsv(ifft)=(rtnewt-vtrial(ifft,1))/dtset%tphysel
270 : ! eektemp=zfermi32(betamumoinsV(ifft))/zfermi12(betamumoinsV(ifft))
271 12255408 : eektemp=fp32a1(betamumoinsV(ifft))/rhor(ifft,1)*cktf
272 12255408 : feektemp=betamumoinsV(ifft)-two/three*eektemp
273 12255408 : feektf=feektf+feektemp*rhor(ifft,1)
274 12255685 : eektf=eektf+eektemp*rhor(ifft,1)
275 : end do
276 : ! Init mpi_comm
277 277 : call timab(48,1,tsec)
278 : call xmpi_sum(eektf,mpi_enreg%comm_fft ,ierr)
279 : call xmpi_sum(feektf,mpi_enreg%comm_fft ,ierr)
280 277 : call timab(48,2,tsec)
281 277 : eektf=eektf*dtset%tphysel
282 277 : eektf=eektf*ucvol/dble(nfft*nproc_fft)
283 277 : feektf=feektf*dtset%tphysel
284 277 : feektf=feektf*ucvol/dble(nfft*nproc_fft)
285 : ! DEBUG
286 : ! write(std_out,*)'eektf',eektf
287 : ! stop ('vtorhotf')
288 : ! ENDDEBUG
289 277 : ek=eektf
290 277 : entropy=(eektf-feektf)/dtset%tphysel
291 277 : ABI_FREE(betamumoinsV)
292 277 : end subroutine tfek
293 : !!***
294 :
295 : !!****f* ABINIT/zfermim12
296 : !! NAME
297 : !! zfermim12
298 : !!
299 : !! FUNCTION
300 : !!..file contains fermi-dirac integral routines:
301 : !!..
302 : !!..function zfermim12 does a rational function fit for the order -1/2 integral
303 : !!..function zfermi12 does a rational function fit for the order 1/2 integral
304 : !!..function zfermi1 does a rational function fit for the order 1 integral
305 : !!..function zfermi32 does a rational function fit for the order 3/2 integral
306 : !!..function zfermi2 does a rational function fit for the order 2 integral
307 : !!..function zfermi52 does a rational function fit for the order 5/2 integral
308 : !!..function zfermi3 does a rational function fit for the order 3 integral
309 : !!..
310 : !!..function ifermim12 is a rational function fit for the inverse of order -1/2
311 : !!..function ifermi12 is a rational function fit for the inverse of order 1/2
312 : !!..function ifermi32 is a rational function fit for the inverse of order 3/2
313 : !!..function ifermi52 is a rational function fit for the inverse of order 5/2
314 : !!
315 : !!..this routine applies a rational function expansion to get the fermi-dirac
316 : !!..integral of order -1/2 evaluated at x. maximum error is 1.23d-12.
317 : !!..reference: antia apjs 84,101 1993
318 : !!
319 : !! INPUTS
320 : !!
321 : !! OUTPUT
322 : !!
323 : !! SOURCE
324 :
325 : function zfermim12(xx)
326 :
327 : !Arguments -------------------------------
328 : real(dp), intent(in) :: xx
329 : real(dp) :: zfermim12
330 :
331 : !Local variables-------------------------------
332 : integer :: ii,m1,k1,m2,k2
333 : real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
334 : !..load the coefficients of the expansion
335 : data an,m1,k1,m2,k2 /-0.5e0_dp, 7, 7, 11, 11/
336 : data (a1(ii),ii=1,8)/ 1.71446374704454e7_dp, 3.88148302324068e7_dp,&
337 : & 3.16743385304962e7_dp, 1.14587609192151e7_dp,&
338 : & 1.83696370756153E6_dp, 1.14980998186874e5_dp,&
339 : & 1.98276889924768e3_dp, 1.0e0_dp/
340 : data (b1(ii),ii=1,8)/ 9.67282587452899e6_dp, 2.87386436731785e7_dp,&
341 : & 3.26070130734158e7_dp, 1.77657027846367e7_dp,&
342 : & 4.81648022267831e6_dp, 6.13709569333207e5_dp,&
343 : & 3.13595854332114e4_dp, 4.35061725080755e2_dp/
344 : data (a2(ii),ii=1,12)/-4.46620341924942e-15_dp, -1.58654991146236e-12_dp,&
345 : & -4.44467627042232e-10_dp, -6.84738791621745e-8_dp,&
346 : & -6.64932238528105e-6_dp, -3.69976170193942e-4_dp,&
347 : & -1.12295393687006e-2_dp, -1.60926102124442e-1_dp,&
348 : & -8.52408612877447e-1_dp, -7.45519953763928e-1_dp,&
349 : & 2.98435207466372e0_dp, 1.0e0_dp/
350 : data (b2(ii),ii=1,12)/-2.23310170962369e-15_dp, -7.94193282071464e-13_dp,&
351 : & -2.22564376956228e-10_dp, -3.43299431079845e-8_dp,&
352 : & -3.33919612678907e-6_dp, -1.86432212187088e-4_dp,&
353 : & -5.69764436880529e-3_dp, -8.34904593067194e-2_dp,&
354 : & -4.78770844009440e-1_dp, -4.99759250374148e-1_dp,&
355 : & 1.86795964993052e0_dp, 4.16485970495288e-1_dp/
356 :
357 : ! *************************************************************************
358 :
359 : if (xx .lt. 2.0e0_dp) then
360 : xx1 = exp(xx)
361 : rn = xx1 + a1(m1)
362 : do ii=m1-1,1,-1
363 : rn = rn*xx1 + a1(ii)
364 : end do
365 : den = b1(k1+1)
366 : do ii=k1,1,-1
367 : den = den*xx1 + b1(ii)
368 : end do
369 : zfermim12 = xx1 * rn/den
370 : ! ..
371 : else
372 : xx1 = one/(xx*xx)
373 : rn = xx1 + a2(m2)
374 : do ii=m2-1,1,-1
375 : rn = rn*xx1 + a2(ii)
376 : end do
377 : den = b2(k2+1)
378 : do ii=k2,1,-1
379 : den = den*xx1 + b2(ii)
380 : end do
381 : zfermim12 = sqrt(xx)*rn/den
382 : end if
383 :
384 : end function zfermim12
385 : !!***
386 :
387 : !!****f* ABINIT/zfermi12
388 : !! NAME
389 : !! zfermi12
390 : !!
391 : !! FUNCTION
392 : !!
393 : !! INPUTS
394 : !!
395 : !! OUTPUT
396 : !!
397 : !! SOURCE
398 :
399 : function zfermi12(xx)
400 : !..
401 : !..this routine applies a rational function expansion to get the fermi-dirac
402 : !..integral of order 1/2 evaluated at x. maximum error is 5.47d-13.
403 : !..reference: antia apjs 84,101 1993
404 : !..
405 : !..declare
406 :
407 : !Arguments -------------------------------
408 : real(dp), intent(in) :: xx
409 : real(dp):: zfermi12
410 :
411 : !Local variables-------------------------------
412 : integer :: ii,m1,k1,m2,k2
413 : real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
414 :
415 : !..load the coefficients of the expansion
416 : data an,m1,k1,m2,k2 /0.5e0_dp, 7, 7, 10, 11/
417 : data (a1(ii),ii=1,8)/5.75834152995465e6_dp, 1.30964880355883e7_dp,&
418 : & 1.07608632249013e7_dp, 3.93536421893014e6_dp,&
419 : & 6.42493233715640e5_dp, 4.16031909245777e4_dp,&
420 : & 7.77238678539648e2_dp, 1.0e0_dp/
421 : data (b1(ii),ii=1,8)/6.49759261942269e6_dp, 1.70750501625775e7_dp,&
422 : & 1.69288134856160e7_dp, 7.95192647756086e6_dp,&
423 : & 1.83167424554505e6_dp, 1.95155948326832e5_dp,&
424 : & 8.17922106644547e3_dp, 9.02129136642157e1_dp/
425 : data (a2(ii),ii=1,11)/4.85378381173415e-14_dp, 1.64429113030738e-11_dp,&
426 : & 3.76794942277806e-9_dp, 4.69233883900644e-7_dp,&
427 : & 3.40679845803144e-5_dp, 1.32212995937796e-3_dp,&
428 : & 2.60768398973913e-2_dp, 2.48653216266227e-1_dp,&
429 : & 1.08037861921488e0_dp, 1.91247528779676e0_dp,&
430 : & 1.0e0_dp/
431 : data (b2(ii),ii=1,12)/7.28067571760518e-14_dp, 2.45745452167585e-11_dp,&
432 : & 5.62152894375277e-9_dp, 6.96888634549649e-7_dp,&
433 : & 5.02360015186394e-5_dp, 1.92040136756592e-3_dp,&
434 : & 3.66887808002874e-2_dp, 3.24095226486468e-1_dp,&
435 : & 1.16434871200131e0_dp, 1.34981244060549e0_dp,&
436 : & 2.01311836975930e-1_dp, -2.14562434782759e-2_dp/
437 :
438 : ! *************************************************************************
439 :
440 : if (xx .lt. two) then
441 : xx1 = exp(xx)
442 : rn = xx1 + a1(m1)
443 : do ii=m1-1,1,-1
444 : rn = rn*xx1 + a1(ii)
445 : end do
446 : den = b1(k1+1)
447 : do ii=k1,1,-1
448 : den = den*xx1 + b1(ii)
449 : end do
450 : zfermi12 = xx1 * rn/den
451 :
452 : else
453 : xx1 = one/(xx*xx)
454 : rn = xx1 + a2(m2)
455 : do ii=m2-1,1,-1
456 : rn = rn*xx1 + a2(ii)
457 : end do
458 : den = b2(k2+1)
459 : do ii=k2,1,-1
460 : den = den*xx1 + b2(ii)
461 : end do
462 : zfermi12 = xx*sqrt(xx)*rn/den
463 : end if
464 :
465 : end function zfermi12
466 : !!***
467 :
468 : !!****f* ABINIT/zfermi1
469 : !! NAME
470 : !! zfermi1
471 : !!
472 : !! FUNCTION
473 : !!
474 : !! INPUTS
475 : !!
476 : !! OUTPUT
477 : !!
478 : !! SOURCE
479 :
480 : function zfermi1(xx)
481 : !..
482 : !..this routine applies a rational function expansion to get the fermi-dirac
483 : !..integral of order 1 evaluated at x. maximum error is 1.0e-8.
484 : !..reference: antia priv comm. 11sep94
485 : !..
486 : !..declare
487 :
488 : !Arguments -------------------------------
489 : real(dp), intent(in) :: xx
490 : real(dp):: zfermi1
491 : !Local variables-------------------------------
492 : integer :: ii,m1,k1,m2,k2
493 : real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
494 :
495 : !..load the coefficients of the expansion
496 : data an,m1,k1,m2,k2 /1.0_dp, 7, 4, 9, 5/
497 : data (a1(ii),ii=1,8)/-7.606458638543e7_dp, -1.143519707857e8_dp,&
498 : & -5.167289383236e7_dp, -7.304766495775e6_dp,&
499 : & -1.630563622280e5_dp, 3.145920924780e3_dp,&
500 : & -7.156354090495e1_dp, 1.0_dp/
501 : data (b1(ii),ii=1,5)/-7.606458639561e7_dp, -1.333681162517e8_dp,&
502 : & -7.656332234147e7_dp, -1.638081306504e7_dp,&
503 : & -1.044683266663e6_dp/
504 : data (a2(ii),ii=1,10)/-3.493105157219e-7_dp, -5.628286279892e-5_dp,&
505 : & -5.188757767899e-3_dp, -2.097205947730e-1_dp,&
506 : & -3.353243201574_dp, -1.682094530855e1_dp,&
507 : & -2.042542575231e1_dp, 3.551366939795_dp,&
508 : & -2.400826804233_dp, 1.0_dp/
509 : data (b2(ii),ii=1,6)/-6.986210315105e-7_dp, -1.102673536040e-4_dp,&
510 : & -1.001475250797e-2_dp, -3.864923270059e-1_dp,&
511 : & -5.435619477378_dp, -1.563274262745e1_dp/
512 :
513 : ! *************************************************************************
514 :
515 : if (xx .lt. 2.0_dp) then
516 : xx1 = exp(xx)
517 : rn = xx1 + a1(m1)
518 : do ii=m1-1,1,-1
519 : rn = rn*xx1 + a1(ii)
520 : end do
521 : den = b1(k1+1)
522 : do ii=k1,1,-1
523 : den = den*xx1 + b1(ii)
524 : end do
525 : zfermi1 = xx1 * rn/den
526 :
527 : else
528 : xx1 = 1.0_dp/(xx*xx)
529 : rn = xx1 + a2(m2)
530 : do ii=m2-1,1,-1
531 : rn = rn*xx1 + a2(ii)
532 : end do
533 : den = b2(k2+1)
534 : do ii=k2,1,-1
535 : den = den*xx1 + b2(ii)
536 : end do
537 : zfermi1 = xx*xx*rn/den
538 : end if
539 :
540 : end function zfermi1
541 : !!***
542 :
543 : !!****f* ABINIT/zfermi32
544 : !! NAME
545 : !! zfermi32
546 : !!
547 : !! FUNCTION
548 : !! this routine applies a rational function expansion to get the fermi-dirac
549 : !! integral of order 3/2 evaluated at x. maximum error is 5.07d-13.
550 : !! reference: antia apjs 84,101 1993
551 : !!
552 : !! INPUTS
553 : !!
554 : !! OUTPUT
555 : !!
556 : !! SOURCE
557 :
558 : function zfermi32(xx)
559 :
560 : !Arguments -------------------------------
561 : real(dp), intent(in) :: xx
562 : real(dp) :: zfermi32
563 :
564 : !Local variables-------------------------------
565 : integer :: ii,m1,k1,m2,k2
566 : real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
567 :
568 : !..load the coefficients of the expansion
569 : data an,m1,k1,m2,k2 /1.5e0_dp, 6, 7, 9, 10/
570 : data (a1(ii),ii=1,7)/4.32326386604283e4_dp, 8.55472308218786e4_dp,&
571 : & 5.95275291210962e4_dp, 1.77294861572005e4_dp,&
572 : & 2.21876607796460e3_dp, 9.90562948053193e1_dp,&
573 : & 1.0e0_dp/
574 : data (b1(ii),ii=1,8)/3.25218725353467e4_dp, 7.01022511904373e4_dp,&
575 : & 5.50859144223638e4_dp, 1.95942074576400e4_dp,&
576 : & 3.20803912586318e3_dp, 2.20853967067789e2_dp,&
577 : & 5.05580641737527e0_dp, 1.99507945223266e-2_dp/
578 : data (a2(ii),ii=1,10)/2.80452693148553e-13_dp, 8.60096863656367e-11_dp,&
579 : & 1.62974620742993e-8_dp, 1.63598843752050e-6_dp,&
580 : & 9.12915407846722e-5_dp, 2.62988766922117e-3_dp,&
581 : & 3.85682997219346e-2_dp, 2.78383256609605e-1_dp,&
582 : & 9.02250179334496e-1_dp, 1.0e0_dp/
583 : data (b2(ii),ii=1,11)/7.01131732871184e-13_dp, 2.10699282897576e-10_dp,&
584 : & 3.94452010378723e-8_dp, 3.84703231868724e-6_dp,&
585 : & 2.04569943213216e-4_dp, 5.31999109566385e-3_dp,&
586 : & 6.39899717779153e-2_dp, 3.14236143831882e-1_dp,&
587 : & 4.70252591891375e-1_dp, -2.15540156936373e-2_dp,&
588 : & 2.34829436438087e-3_dp/
589 :
590 : ! *************************************************************************
591 :
592 : if (xx .lt. 2.0e0_dp) then
593 : xx1 = exp(xx)
594 : rn = xx1 + a1(m1)
595 : do ii=m1-1,1,-1
596 : rn = rn*xx1 + a1(ii)
597 : end do
598 : den = b1(k1+1)
599 : do ii=k1,1,-1
600 : den = den*xx1 + b1(ii)
601 : end do
602 : zfermi32 = xx1 * rn/den
603 :
604 : else
605 : xx1 = one/(xx*xx)
606 : rn = xx1 + a2(m2)
607 : do ii=m2-1,1,-1
608 : rn = rn*xx1 + a2(ii)
609 : end do
610 : den = b2(k2+1)
611 : do ii=k2,1,-1
612 : den = den*xx1 + b2(ii)
613 : end do
614 : zfermi32 = xx*xx*sqrt(xx)*rn/den
615 : end if
616 :
617 : end function zfermi32
618 : !!***
619 :
620 : !!****f* ABINIT/zfermi2
621 : !! NAME
622 : !! zfermi2
623 : !!
624 : !! FUNCTION
625 : !!
626 : !! this routine applies a rational function expansion to get the fermi-dirac
627 : !! integral of order 2 evaluated at x. maximum error is 1.0e-8.
628 : !! reference: antia priv comm. 11sep94
629 : !!
630 : !! INPUTS
631 : !!
632 : !! OUTPUT
633 : !!
634 : !! SOURCE
635 :
636 : function zfermi2(xx)
637 :
638 : !Arguments -------------------------------
639 : real(dp), intent(in) :: xx
640 : real(dp) :: zfermi2
641 : !Local variables-------------------------------
642 : integer :: ii,m1,k1,m2,k2
643 : real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
644 :
645 : !..load the coefficients of the expansion
646 : data an,m1,k1,m2,k2 /2.0_dp, 7, 4, 5, 9/
647 : data (a1(ii),ii=1,8)/-1.434885992395e8_dp, -2.001711155617e8_dp,&
648 : & -8.507067153428e7_dp, -1.175118281976e7_dp,&
649 : & -3.145120854293e5_dp, 4.275771034579e3_dp,&
650 : & -8.069902926891e1_dp, 1.0e0_dp/
651 : data (b1(ii),ii=1,5)/-7.174429962316e7_dp, -1.090535948744e8_dp,&
652 : & -5.350984486022e7_dp, -9.646265123816e6_dp,&
653 : & -5.113415562845e5_dp/
654 : data (a2(ii),ii=1,6)/ 6.919705180051e-8_dp, 1.134026972699e-5_dp,&
655 : & 7.967092675369e-4_dp, 2.432500578301e-2_dp,&
656 : & 2.784751844942e-1_dp, 1.0e0_dp/
657 : data (b2(ii),ii=1,10)/ 2.075911553728e-7_dp, 3.197196691324e-5_dp,&
658 : & 2.074576609543e-3_dp, 5.250009686722e-2_dp,&
659 : & 3.171705130118e-1_dp, -1.147237720706e-1_dp,&
660 : & 6.638430718056e-2_dp, -1.356814647640e-2_dp,&
661 : & -3.648576227388e-2_dp, 3.621098757460e-2_dp/
662 :
663 : ! *************************************************************************
664 :
665 : if (xx .lt. 2.0e0_dp) then
666 : xx1 = exp(xx)
667 : rn = xx1 + a1(m1)
668 : do ii=m1-1,1,-1
669 : rn = rn*xx1 + a1(ii)
670 : end do
671 : den = b1(k1+1)
672 : do ii=k1,1,-1
673 : den = den*xx1 + b1(ii)
674 : end do
675 : zfermi2 = xx1 * rn/den
676 :
677 : else
678 : xx1 = one/(xx*xx)
679 : rn = xx1 + a2(m2)
680 : do ii=m2-1,1,-1
681 : rn = rn*xx1 + a2(ii)
682 : end do
683 : den = b2(k2+1)
684 : do ii=k2,1,-1
685 : den = den*xx1 + b2(ii)
686 : end do
687 : zfermi2 = xx*xx*xx*rn/den
688 : end if
689 :
690 : end function zfermi2
691 : !!***
692 :
693 : !!****f* ABINIT/zfermi52
694 : !! NAME
695 : !! zfermi52
696 : !!
697 : !! FUNCTION
698 : !! this routine applies a rational function expansion to get the fermi-dirac
699 : !! integral of order 5/2 evaluated at x. maximum error is 2.47d-13.
700 : !! reference: antia apjs 84,101 1993
701 : !!
702 : !! INPUTS
703 : !!
704 : !! OUTPUT
705 : !!
706 : !! SOURCE
707 :
708 : function zfermi52(xx)
709 :
710 : !Arguments -------------------------------
711 : real(dp), intent(in) :: xx
712 : real(dp) :: zfermi52
713 :
714 : !Local variables-------------------------------
715 : integer :: ii,m1,k1,m2,k2
716 : real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
717 :
718 : !..load the coefficients of the expansion
719 : data an,m1,k1,m2,k2 /2.5e0_dp, 6, 7, 10, 9/
720 : data (a1(ii),ii=1,7)/6.61606300631656e4_dp, 1.20132462801652e5_dp,&
721 : & 7.67255995316812e4_dp, 2.10427138842443e4_dp,&
722 : & 2.44325236813275e3_dp, 1.02589947781696e2_dp,&
723 : & 1.0e0_dp/
724 : data (b1(ii),ii=1,8)/1.99078071053871e4_dp, 3.79076097261066e4_dp,&
725 : & 2.60117136841197e4_dp, 7.97584657659364e3_dp,&
726 : & 1.10886130159658e3_dp, 6.35483623268093e1_dp,&
727 : & 1.16951072617142e0_dp, 3.31482978240026e-3_dp/
728 : data (a2(ii),ii=1,11)/8.42667076131315e-12_dp, 2.31618876821567e-9_dp,&
729 : & 3.54323824923987e-7_dp, 2.77981736000034e-5_dp,&
730 : & 1.14008027400645e-3_dp, 2.32779790773633e-2_dp,&
731 : & 2.39564845938301e-1_dp, 1.24415366126179e0_dp,&
732 : & 3.18831203950106e0_dp, 3.42040216997894e0_dp,&
733 : & 1.0e0_dp/
734 : data (b2(ii),ii=1,10)/2.94933476646033e-11_dp, 7.68215783076936e-9_dp,&
735 : & 1.12919616415947e-6_dp, 8.09451165406274e-5_dp,&
736 : & 2.81111224925648e-3_dp, 3.99937801931919e-2_dp,&
737 : & 2.27132567866839e-1_dp, 5.31886045222680e-1_dp,&
738 : & 3.70866321410385e-1_dp, 2.27326643192516e-2_dp/
739 :
740 : ! *************************************************************************
741 :
742 : if (xx .lt. two) then
743 : xx1 = exp(xx)
744 : rn = xx1 + a1(m1)
745 : do ii=m1-1,1,-1
746 : rn = rn*xx1 + a1(ii)
747 : end do
748 : den = b1(k1+1)
749 : do ii=k1,1,-1
750 : den = den*xx1 + b1(ii)
751 : end do
752 : zfermi52 = xx1 * rn/den
753 :
754 : else
755 : xx1 = one/(xx*xx)
756 : rn = xx1 + a2(m2)
757 : do ii=m2-1,1,-1
758 : rn = rn*xx1 + a2(ii)
759 : end do
760 : den = b2(k2+1)
761 : do ii=k2,1,-1
762 : den = den*xx1 + b2(ii)
763 : end do
764 : zfermi52 = xx*xx*xx*sqrt(xx)*rn/den
765 : end if
766 :
767 : end function zfermi52
768 : !!***
769 :
770 : !!****f* ABINIT/zfermi3
771 : !! NAME
772 : !! zfermi3
773 : !!
774 : !! FUNCTION
775 : !! this routine applies a rational function expansion to get the fermi-dirac
776 : !! integral of order 3 evaluated at x. maximum error is 1.0e-8.
777 : !! reference: antia priv comm. 11sep94
778 : !!
779 : !! INPUTS
780 : !!
781 : !! OUTPUT
782 : !!
783 : !! SOURCE
784 :
785 : function zfermi3(xx)
786 :
787 : !Arguments -------------------------------
788 : real(dp), intent(in) :: xx
789 : real(dp):: zfermi3
790 :
791 : !Local variables-------------------------------
792 : integer :: ii,m1,k1,m2,k2
793 : real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx1
794 :
795 : !..load the coefficients of the expansion
796 : data an,m1,k1,m2,k2 /3.0, 4, 6, 7, 7/
797 : data (a1(ii),ii=1,5)/ 6.317036716422e2_dp, 7.514163924637e2_dp,&
798 : & 2.711961035750e2_dp, 3.274540902317e1_dp,&
799 : & 1.0_dp/
800 : data (b1(ii),ii=1,7)/ 1.052839452797e2_dp, 1.318163114785e2_dp,&
801 : & 5.213807524405e1_dp, 7.500064111991_dp,&
802 : & 3.383020205492e-1_dp, 2.342176749453e-3_dp,&
803 : & -8.445226098359e-6_dp/
804 : data (a2(ii),ii=1,8)/ 1.360999428425e-8_dp, 1.651419468084e-6_dp,&
805 : & 1.021455604288e-4_dp, 3.041270709839e-3_dp,&
806 : & 4.584298418374e-2_dp, 3.440523212512e-1_dp,&
807 : & 1.077505444383_dp, 1.0_dp/
808 : data (b2(ii),ii=1,8)/ 5.443997714076e-8_dp, 5.531075760054e-6_dp,&
809 : & 2.969285281294e-4_dp, 6.052488134435e-3_dp,&
810 : & 5.041144894964e-2_dp, 1.048282487684e-1_dp,&
811 : & 1.280969214096e-2_dp, -2.851555446444e-3_dp/
812 :
813 : ! *************************************************************************
814 :
815 : if (xx .lt. two) then
816 : xx1 = exp(xx)
817 : rn = xx1 + a1(m1)
818 : do ii=m1-1,1,-1
819 : rn = rn*xx1 + a1(ii)
820 : end do
821 : den = b1(k1+1)
822 : do ii=k1,1,-1
823 : den = den*xx1 + b1(ii)
824 : end do
825 : zfermi3 = xx1 * rn/den
826 :
827 : else
828 : xx1 = one/(xx*xx)
829 : rn = xx1 + a2(m2)
830 : do ii=m2-1,1,-1
831 : rn = rn*xx1 + a2(ii)
832 : end do
833 : den = b2(k2+1)
834 : do ii=k2,1,-1
835 : den = den*xx1 + b2(ii)
836 : end do
837 : zfermi3 = xx*xx*xx*xx*rn/den
838 : end if
839 :
840 : end function zfermi3
841 : !!***
842 :
843 : !!****f* ABINIT/ifermim12
844 : !! NAME
845 : !! ifermim12
846 : !!
847 : !! FUNCTION
848 : !! this routine applies a rational function expansion to get the inverse
849 : !! fermi-dirac integral of order -1/2 when it is equal to f.
850 : !! maximum error is 3.03d-9. reference: antia apjs 84,101 1993
851 : !!
852 : !! INPUTS
853 : !!
854 : !! OUTPUT
855 : !!
856 : !! SOURCE
857 :
858 : function ifermim12(ff)
859 :
860 : !Arguments -------------------------------
861 : real(dp), intent(in) :: ff
862 : real(dp) :: ifermim12
863 : !Local variables-------------------------------
864 : integer :: ii,m1,k1,m2,k2
865 : real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,ff1
866 :
867 : !..load the coefficients of the expansion
868 : data an,m1,k1,m2,k2 /-0.5e0_dp, 5, 6, 6, 6/
869 : data (a1(ii),ii=1,6)/-1.570044577033e4_dp, 1.001958278442e4_dp,&
870 : & -2.805343454951e3_dp, 4.121170498099e2_dp,&
871 : & -3.174780572961e1_dp, 1.0e0_dp/
872 : data (b1(ii),ii=1,7)/-2.782831558471e4_dp, 2.886114034012e4_dp,&
873 : & -1.274243093149e4_dp, 3.063252215963e3_dp,&
874 : & -4.225615045074e2_dp, 3.168918168284e1_dp,&
875 : & -1.008561571363e0_dp/
876 : data (a2(ii),ii=1,7)/ 2.206779160034e-8_dp, -1.437701234283e-6_dp,&
877 : & 6.103116850636e-5_dp, -1.169411057416e-3_dp,&
878 : & 1.814141021608e-2_dp, -9.588603457639e-2_dp,&
879 : & 1.0e0_dp/
880 : data (b2(ii),ii=1,7)/ 8.827116613576e-8_dp, -5.750804196059e-6_dp,&
881 : & 2.429627688357e-4_dp, -4.601959491394e-3_dp,&
882 : & 6.932122275919e-2_dp, -3.217372489776e-1_dp,&
883 : & 3.124344749296e0_dp/
884 :
885 : ! *************************************************************************
886 :
887 : if (ff .lt. 4.0e0_dp) then
888 : rn = ff + a1(m1)
889 : do ii=m1-1,1,-1
890 : rn = rn*ff + a1(ii)
891 : end do
892 : den = b1(k1+1)
893 : do ii=k1,1,-1
894 : den = den*ff + b1(ii)
895 : end do
896 : ifermim12 = log(ff * rn/den)
897 :
898 : else
899 : ff1 = one/ff**(one/(one + an))
900 : rn = ff1 + a2(m2)
901 : do ii=m2-1,1,-1
902 : rn = rn*ff1 + a2(ii)
903 : end do
904 : den = b2(k2+1)
905 : do ii=k2,1,-1
906 : den = den*ff1 + b2(ii)
907 : end do
908 : ifermim12 = rn/(den*ff1)
909 : end if
910 :
911 : end function ifermim12
912 : !!***
913 :
914 : !!****f* ABINIT/ifermi12
915 : !! NAME
916 : !! ifermi12
917 : !!
918 : !! FUNCTION
919 : !! this routine applies a rational function expansion to get the inverse
920 : !! fermi-dirac integral of order 1/2 when it is equal to f.
921 : !! maximum error is 4.19d-9. reference: antia apjs 84,101 1993
922 : !!
923 : !! INPUTS
924 : !!
925 : !! OUTPUT
926 : !!
927 : !! SOURCE
928 :
929 : function ifermi12(ff)
930 :
931 : !Arguments -------------------------------
932 : real(dp), intent(in) :: ff
933 : real(dp) :: ifermi12
934 : !Local variables-------------------------------
935 : integer :: ii,m1,k1,m2,k2
936 : real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,ff1
937 :
938 : !..load the coefficients of the expansion
939 : data an,m1,k1,m2,k2 /0.5e0_dp, 4, 3, 6, 5/
940 : data (a1(ii),ii=1,5)/ 1.999266880833e4_dp, 5.702479099336e3_dp,&
941 : & 6.610132843877e2_dp, 3.818838129486e1_dp,&
942 : & 1.0e0_dp/
943 : data (b1(ii),ii=1,4)/ 1.771804140488e4_dp, -2.014785161019e3_dp,&
944 : & 9.130355392717e1_dp, -1.670718177489e0_dp/
945 : data (a2(ii),ii=1,7)/-1.277060388085e-2_dp, 7.187946804945e-2_dp,&
946 : & -4.262314235106e-1_dp, 4.997559426872e-1_dp,&
947 : & -1.285579118012e0_dp, -3.930805454272e-1_dp,&
948 : & 1.0e0_dp/
949 : data (b2(ii),ii=1,6)/-9.745794806288e-3_dp, 5.485432756838e-2_dp,&
950 : & -3.299466243260e-1_dp, 4.077841975923e-1_dp,&
951 : & -1.145531476975e0_dp, -6.067091689181e-2_dp/
952 :
953 : ! *************************************************************************
954 :
955 : if (ff .lt. 4.0e0_dp) then
956 : rn = ff + a1(m1)
957 : do ii=m1-1,1,-1
958 : rn = rn*ff + a1(ii)
959 : end do
960 : den = b1(k1+1)
961 : do ii=k1,1,-1
962 : den = den*ff + b1(ii)
963 : end do
964 : ifermi12 = log(ff * rn/den)
965 :
966 : else
967 : ff1 = one/ff**(one/(one + an))
968 : rn = ff1 + a2(m2)
969 : do ii=m2-1,1,-1
970 : rn = rn*ff1 + a2(ii)
971 : end do
972 : den = b2(k2+1)
973 : do ii=k2,1,-1
974 : den = den*ff1 + b2(ii)
975 : end do
976 : ifermi12 = rn/(den*ff1)
977 : end if
978 :
979 : end function ifermi12
980 : !!***
981 :
982 : !!****f* ABINIT/ifermi32
983 : !! NAME
984 : !! ifermi32
985 : !!
986 : !! FUNCTION
987 : !! this routine applies a rational function expansion to get the inverse
988 : !! fermi-dirac integral of order 3/2 when it is equal to f.
989 : !! maximum error is 2.26d-9. reference: antia apjs 84,101 1993
990 : !!
991 : !! INPUTS
992 : !!
993 : !! OUTPUT
994 : !!
995 : !! SOURCE
996 :
997 : function ifermi32(ff)
998 :
999 : !Arguments -------------------------------
1000 : real(dp), intent(in) :: ff
1001 : real(dp) :: ifermi32
1002 : !Local variables-------------------------------
1003 : integer :: ii,m1,k1,m2,k2
1004 : real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,ff1
1005 :
1006 : !..load the coefficients of the expansion
1007 : data an,m1,k1,m2,k2 /1.5e0_dp, 3, 4, 6, 5/
1008 : data (a1(ii),ii=1,4)/ 1.715627994191e2_dp, 1.125926232897e2_dp,&
1009 : & 2.056296753055e1_dp, 1.0e0_dp/
1010 : data (b1(ii),ii=1,5)/ 2.280653583157e2_dp, 1.193456203021e2_dp,&
1011 : & 1.167743113540e1_dp, -3.226808804038e-1_dp,&
1012 : & 3.519268762788e-3_dp/
1013 : data (a2(ii),ii=1,7)/-6.321828169799e-3_dp, -2.183147266896e-2_dp,&
1014 : & -1.057562799320e-1_dp, -4.657944387545e-1_dp,&
1015 : & -5.951932864088e-1_dp, 3.684471177100e-1_dp,&
1016 : & 1.0e0_dp/
1017 : data (b2(ii),ii=1,6)/-4.381942605018e-3_dp, -1.513236504100e-2_dp,&
1018 : & -7.850001283886e-2_dp, -3.407561772612e-1_dp,&
1019 : & -5.074812565486e-1_dp, -1.387107009074e-1_dp/
1020 :
1021 : ! *************************************************************************
1022 :
1023 : if (ff .lt. 4.0e0_dp) then
1024 : rn = ff + a1(m1)
1025 : do ii=m1-1,1,-1
1026 : rn = rn*ff + a1(ii)
1027 : end do
1028 : den = b1(k1+1)
1029 : do ii=k1,1,-1
1030 : den = den*ff + b1(ii)
1031 : end do
1032 : ifermi32 = log(ff * rn/den)
1033 :
1034 : else
1035 : ff1 = one/ff**(one/(one + an))
1036 : rn = ff1 + a2(m2)
1037 : do ii=m2-1,1,-1
1038 : rn = rn*ff1 + a2(ii)
1039 : end do
1040 : den = b2(k2+1)
1041 : do ii=k2,1,-1
1042 : den = den*ff1 + b2(ii)
1043 : end do
1044 : ifermi32 = rn/(den*ff1)
1045 : end if
1046 :
1047 : end function ifermi32
1048 : !!***
1049 :
1050 : !!****f* ABINIT/ifermi52
1051 : !! NAME
1052 : !! ifermi52
1053 : !!
1054 : !! FUNCTION
1055 : !! this routine applies a rational function expansion to get the inverse
1056 : !! fermi-dirac integral of order 5/2 when it is equal to f.
1057 : !! maximum error is 6.17d-9. reference: antia apjs 84,101 1993
1058 : !!
1059 : !! INPUTS
1060 : !!
1061 : !! OUTPUT
1062 : !!
1063 : !! SOURCE
1064 :
1065 : function ifermi52(ff)
1066 :
1067 : !Arguments -------------------------------
1068 : real(dp), intent(in) :: ff
1069 : real(dp) :: ifermi52
1070 :
1071 : !Local variables-------------------------------
1072 : integer :: ii,m1,k1,m2,k2
1073 : real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,ff1
1074 :
1075 : !..load the coefficients of the expansion
1076 : data an,m1,k1,m2,k2 /2.5e0_dp, 2, 3, 6, 6/
1077 : data (a1(ii),ii=1,3)/ 2.138969250409e2_dp, 3.539903493971e1_dp,&
1078 : & 1.0e0_dp/
1079 : data (b1(ii),ii=1,4)/ 7.108545512710e2_dp, 9.873746988121e1_dp,&
1080 : & 1.067755522895e0_dp, -1.182798726503e-2_dp/
1081 : data (a2(ii),ii=1,7)/-3.312041011227e-2_dp, 1.315763372315e-1_dp,&
1082 : & -4.820942898296e-1_dp, 5.099038074944e-1_dp,&
1083 : & 5.495613498630e-1_dp, -1.498867562255e0_dp,&
1084 : & 1.0e0_dp/
1085 : data (b2(ii),ii=1,7)/-2.315515517515e-2_dp, 9.198776585252e-2_dp,&
1086 : & -3.835879295548e-1_dp, 5.415026856351e-1_dp,&
1087 : & -3.847241692193e-1_dp, 3.739781456585e-2_dp,&
1088 : & -3.008504449098e-2_dp/
1089 :
1090 : ! *************************************************************************
1091 :
1092 : if (ff .lt. 4.0e0_dp) then
1093 : rn = ff + a1(m1)
1094 : do ii=m1-1,1,-1
1095 : rn = rn*ff + a1(ii)
1096 : end do
1097 : den = b1(k1+1)
1098 : do ii=k1,1,-1
1099 : den = den*ff + b1(ii)
1100 : end do
1101 : ifermi52 = log(ff * rn/den)
1102 :
1103 : else
1104 : ff1 = one/ff**(one/(one + an))
1105 : rn = ff1 + a2(m2)
1106 : do ii=m2-1,1,-1
1107 : rn = rn*ff1 + a2(ii)
1108 : end do
1109 : den = b2(k2+1)
1110 : do ii=k2,1,-1
1111 : den = den*ff1 + b2(ii)
1112 : end do
1113 : ifermi52 = rn/(den*ff1)
1114 : end if
1115 :
1116 : end function ifermi52
1117 : !!***
1118 :
1119 : !!****f* ABINIT/fp12a1
1120 : !! NAME
1121 : !! fp12a1
1122 : !!
1123 : !! FUNCTION
1124 : !!
1125 : !! INPUTS
1126 : !!
1127 : !! OUTPUT
1128 : !!
1129 : !! SOURCE
1130 :
1131 : function fp12a1 (x)
1132 :
1133 : ! Arguments -------------------------------
1134 : real(dp),intent(in) :: x
1135 : real(dp) :: fp12a1
1136 :
1137 : real(dp) :: y
1138 :
1139 : !**********************************************************************
1140 : !* *
1141 : !* Integrale de Fermi d'ordre 1/2 *
1142 : !* Fp12(x) = somme de 0 a l'infini de (dt*t**1/2)/(1+exp(t-x)) *
1143 : !* *
1144 : !**********************************************************************
1145 : !
1146 : !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1147 : !Erreur relative maximum annoncee 5.54 e-5
1148 : !Erreur relative maximum constatee : -5.53e-5 pour eta = 2
1149 : !
1150 : if (x.lt.2._dp) then
1151 : y=exp(x)
1152 : fp12a1=y*(21.8168_dp+y*(13.1693_dp+y))&
1153 : & /(24.6180_dp+y*(23.5546_dp+y*(4.76290_dp+y*0.134481_dp)))
1154 : else
1155 : y=one/(x*x)
1156 : fp12a1=x*sqrt(x)*(0.0473011_dp+y*(0.548433_dp+y))&
1157 : & /(0.0709478_dp+y*(0.737041_dp+y*0.382065_dp))
1158 : end if
1159 : !
1160 : !**********************************************************************
1161 : end function fp12a1
1162 : !!***
1163 :
1164 : !!****f* ABINIT/fp32a1
1165 : !! NAME
1166 : !! fp32a1
1167 : !!
1168 : !! FUNCTION
1169 : !!
1170 : !! INPUTS
1171 : !!
1172 : !! OUTPUT
1173 : !!
1174 : !! SOURCE
1175 :
1176 12255408 : function fp32a1 (x)
1177 :
1178 : !Arguments -------------------------------
1179 : real(dp),intent(in) :: x
1180 : real(dp) :: fp32a1
1181 :
1182 : real(dp) :: y,x2
1183 :
1184 : !
1185 : !**********************************************************************
1186 : !* *
1187 : !* Integrale de Fermi d'ordre 3/2 *
1188 : !* Fp32(x) = somme de 0 a l'infini de (dt*t**3/2)/(1+exp(t-x)) *
1189 : !* *
1190 : !**********************************************************************
1191 : !
1192 : !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1193 : !Erreur relative maximum annoncee 6.54 e-5
1194 : !Erreur relative maximum constatee : -5.84e-5 pour eta = -5
1195 : !
1196 12255408 : if (x.lt.two) then
1197 0 : y=exp(x)
1198 0 : fp32a1=y*(135.863_dp+y*(49.2764_dp+y))/(102.210_dp+y*(55.0312_dp+y*4.23365_dp))
1199 : else
1200 12255408 : x2=x*x
1201 12255408 : y=1._dp/x2
1202 : fp32a1=x2*sqrt(x)*(0.154699_dp+y*(1.20037_dp+y))&
1203 12255408 : & /(0.386765_dp+y*(0.608119_dp-y*0.165665_dp))
1204 : end if
1205 : !
1206 : !**********************************************************************
1207 12255408 : end function fp32a1
1208 : !!***
1209 :
1210 : !!****f* ABINIT/xp12a1
1211 : !! NAME
1212 : !! xp12a1
1213 : !!
1214 : !! FUNCTION
1215 : !!
1216 : !! INPUTS
1217 : !!
1218 : !! OUTPUT
1219 : !!
1220 : !! SOURCE
1221 :
1222 : function xp12a1 (y)
1223 :
1224 : !Arguments -------------------------------
1225 : real(dp) :: xp12a1
1226 : real(dp),intent(in) :: y
1227 :
1228 : real(dp),parameter :: deux=2._dp,deuxs3=deux/3._dp
1229 : real(dp) :: top,den,z
1230 :
1231 : !
1232 : !**********************************************************************
1233 : !* *
1234 : !* Calcul de eta tel que fp12 (eta) = y *
1235 : !* ou fp12 est l'integrale de Fermi d'ordre +1/2 *
1236 : !* *
1237 : !**********************************************************************
1238 : !
1239 : !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1240 : !Erreur relative maximum annoncee sur exp(eta) : 3.02 e-5
1241 : !
1242 : if (y.lt.4._dp) then
1243 : top=44.593646_dp+y*(11.288764_dp+y)
1244 : den=39.519346_dp+y*(-5.7517464_dp+y*0.26594291_dp)
1245 : xp12a1=log(y*top/den)
1246 : else
1247 : z=y**(-deuxs3)
1248 : top=34.873722_dp+z*(-26.922515_dp+z)
1249 : den=26.612832_dp+z*(-20.452930_dp+z*11.808945_dp)
1250 : xp12a1=top/(z*den)
1251 : end if
1252 : !
1253 : !**********************************************************************
1254 : end function xp12a1
1255 : !!***
1256 :
1257 : !!****f* ABINIT/fm12a1
1258 : !! NAME
1259 : !! fm12a1
1260 : !!
1261 : !! FUNCTION
1262 : !!
1263 : !! INPUTS
1264 : !!
1265 : !! OUTPUT
1266 : !!
1267 : !! SOURCE
1268 :
1269 : function fm12a1 (x)
1270 :
1271 : !Arguments -------------------------------
1272 : real(dp),intent(in) :: x
1273 : real(dp) :: fm12a1
1274 :
1275 : real(dp) :: y
1276 :
1277 : !
1278 : !**********************************************************************
1279 : !* *
1280 : !* Integrale de Fermi d'ordre -1/2 *
1281 : !* Fm12(x) = somme de 0 a l'infini de (dt*t**-1/2)/(1+exp(t-x)) *
1282 : !* *
1283 : !**********************************************************************
1284 : !
1285 : !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1286 : !Erreur relative maximum annoncee 4.75 e-5
1287 : !
1288 : if (x.lt.2._dp) then
1289 : y=exp(x)
1290 : fm12a1=y*(23.1456_dp+y*(13.7820_dp+y))&
1291 : & /(13.0586_dp+y*(17.0048_dp+y*(5.07527_dp+y*0.236620_dp)))
1292 : else
1293 : y=1./(x*x)
1294 : fm12a1=sqrt(x)*(0.0153602_dp+y*(0.146815_dp+y))&
1295 : & /(0.00768015_dp+y*(0.0763700_dp+y*0.570485_dp))
1296 : end if
1297 : !
1298 : !**********************************************************************
1299 : end function fm12a1
1300 : !!***
1301 :
1302 : !!****f* ABINIT/fm12a1t
1303 : !! NAME
1304 : !! fm12a1t
1305 : !!
1306 : !! FUNCTION
1307 : !!
1308 : !! INPUTS
1309 : !!
1310 : !! OUTPUT
1311 : !!
1312 : !! SOURCE
1313 :
1314 2182 : subroutine fm12a1t (cktf,rtnewt,tphysel,vtrial,rhor_middx,rhor_mid,nfft)
1315 :
1316 : integer,intent(in) :: nfft
1317 : real(dp),intent(in) :: tphysel,rtnewt,cktf
1318 : real(dp),intent(in) :: vtrial(nfft)
1319 : real(dp),intent(out) :: rhor_middx(nfft),rhor_mid(nfft)
1320 :
1321 : !intrinsic exp,sqrt
1322 : integer :: ifft
1323 : real(dp) :: x,y,sqrtx
1324 :
1325 : !
1326 : !**********************************************************************
1327 : !* *
1328 : !* Integrale de Fermi d'ordre -1/2 *
1329 : !* Fm12(x) = somme de 0 a l'infini de (dt*t**-1/2)/(1+exp(t-x)) *
1330 : !* .... *
1331 : !**********************************************************************
1332 : !
1333 : !H. M. Antia, Astrophys. J. Suppl. 84, 101 (1993)
1334 : !Erreur relative maximum annoncee 4.75 e-5
1335 : !
1336 97127446 : do ifft=1,nfft
1337 97125264 : x=(rtnewt-vtrial(ifft))/tphysel
1338 97127446 : if (x.lt.2._dp) then
1339 9875783 : y=exp(x)
1340 : rhor_middx(ifft)=cktf*y*(23.1456e0_dp+y*(13.7820e0_dp+y))&
1341 9875783 : & /(13.0586e0_dp+y*(17.0048e0_dp+y*(5.07527e0_dp+y*0.236620e0_dp)))
1342 : rhor_mid(ifft)=cktf*y*(21.8168_dp+y*(13.1693_dp+y))&
1343 9875783 : & /(24.6180+y*(23.5546_dp+y*(4.76290_dp+y*0.134481_dp)))
1344 : else
1345 87249481 : y=1._dp/(x*x)
1346 87249481 : sqrtx=sqrt(x)
1347 : rhor_middx(ifft)=cktf*sqrtx*(0.0153602e0_dp+y*(0.146815e0_dp+y))&
1348 87249481 : & /(0.00768015e0_dp+y*(0.0763700e0_dp+y*0.570485e0_dp))
1349 : rhor_mid(ifft)=cktf*x*sqrtx*(0.0473011_dp+y*(0.548433_dp+y))&
1350 87249481 : & /(0.0709478_dp+y*(0.737041_dp+y*0.382065_dp))
1351 : end if
1352 : end do
1353 : !
1354 : !**********************************************************************
1355 2182 : end subroutine fm12a1t
1356 : !!***
1357 :
1358 : end subroutine vtorhotf
1359 : !!***
1360 :
1361 : end module m_vtorhotf
1362 : !!***
|
