Branch data 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 : : !!***
|
