1c 2c$Id$ 3c 4 Subroutine xc_ft97(tol_rho, xfac, lxfac, nlxfac, 5 . cfac, lcfac, nlcfac, 6 , rho, delrho, 7 & Amat, Cmat, nq, ipol, Ex, Ec, qwght, 8 & ldew,func) 9 implicit none 10#include "dft2drv.fh" 11 logical ldew ! [input] 12 logical lxfac, nlxfac ! [input] 13 logical lcfac, nlcfac ! [input] 14 double precision func(*) ![input/output] 15 double precision xfac ![input] 16 double precision cfac ![input] 17c 18 integer ipol ! no. of spin states [input] 19 integer nq ! no. of quadrature pts [input] 20 double precision tol_rho! [input]!threshold on density 21 double precision Ec ! Correlation energy [input/output] 22 double precision Ex ! Exchange energy [input/output] 23 double precision rho(nq,ipol*(ipol+1)/2)! Charge Density [input] 24 double precision delrho(nq,3,ipol) ! Charge Density Gradient[input] 25 double precision qwght(nq) ! Quadrature Weights [input] 26 double precision Amat(nq,ipol) !Sampling Matrices for the XC [output] 27 double precision Cmat(nq,*)!Potential & Energy [output] 28 integer n 29 double precision gammaval 30c to hcth 31 double precision rhoa 32 double precision rhob 33 double precision za 34 double precision zb 35 double precision dfdza,dfdzb 36c 37 double precision fc_ft97,dfdrac,dfdgaac,dfdrbc,dfdgbbc 38 double precision fx_ft97,dfdrax,dfdgaax,dfdrbx,dfdgbbx 39c 40 if(ipol.eq.1) then 41 do n=1,nq 42 if(rho(n,1).gt.tol_rho) then 43 gammaval=delrho(n,1,1)*delrho(n,1,1) + 44 & delrho(n,2,1)*delrho(n,2,1) + 45 & delrho(n,3,1)*delrho(n,3,1) 46 call xc_rks_ft97(rho(n,1),gammaval, 47 * fc_ft97,dfdrac,dfdgaac, 48 * fx_ft97,dfdrax,dfdgaax,tol_rho) 49 if(ldew) func(n)=func(n)+fc_ft97*cfac+ 50 * fx_ft97*xfac 51 Ex=Ex+qwght(n)*xfac*fx_ft97 52 Ec=Ec+qwght(n)*cfac*fc_ft97 53 Amat(n,1) = Amat(n,1)+dfdrax*xfac+dfdrac*cfac 54 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dfdgaac*cfac+ 55 + dfdgaax*xfac 56 endif 57 enddo 58 else 59 do n=1,nq 60 if(rho(n,1).gt.tol_rho) then 61 rhoa=rho(n,2) 62 rhob=rho(n,3) 63 za=delrho(n,1,1)*delrho(n,1,1) + 64 & delrho(n,2,1)*delrho(n,2,1) + 65 & delrho(n,3,1)*delrho(n,3,1) 66 zb=delrho(n,1,2)*delrho(n,1,2) + 67 & delrho(n,2,2)*delrho(n,2,2) + 68 & delrho(n,3,2)*delrho(n,3,2) 69 call xc_uks_ft97(tol_rho, 70 ( rhoa,rhob,za,zb, 71 * fc_ft97,dfdrac,dfdrbc,dfdgaac,dfdgbbc, 72 * fx_ft97,dfdrax,dfdrbx,dfdgaax,dfdgbbx) 73 if(ldew) func(n)=func(n)+fc_ft97*cfac+fx_ft97*xfac 74 Ex=Ex+qwght(n)*xfac*fx_ft97 75 Ec=Ec+qwght(n)*cfac*fc_ft97 76 Amat(n,1) = Amat(n,1)+dfdrax*xfac+dfdrac*cfac 77 Amat(n,2) = Amat(n,2)+dfdrbx*xfac+dfdrbc*cfac 78 dfdza=dfdgaax*xfac+dfdgaac*cfac 79 dfdzb=dfdgbbx*xfac+dfdgbbc*cfac 80 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dfdza 81 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + dfdzb 82 endif 83 enddo 84 endif 85 return 86 end 87cDFT functional repository: xc_ft97 fortran77 source 88CDFT repository Quantum Chemistry Group 89C 90C CCLRC Density functional repository Copyright notice 91C This database was prepared as a result of work undertaken by CCLRC. 92C Users may view, print and download the content for personal use only 93C and the content must not be used for any commercial purpose without CCLRC 94C prior written permission 95C 96 97c----------------------------------------------------------------------- 98 subroutine xc_rks_ft97(r,g,fc,dfdrac,dfdgaac,fx,dfdrax,dfdgaax, 99 , tol_rho) 100c 101c This subroutine calculates the Filatov-Thiel 97 102c exchange-correlation functional [1,2] for the closed shell case. 103c This functional is taken to be the correlation functional [1] plus 104c the recommended variant B of the exchange functional [2]. 105c Also the derivatives with respect to the density components and 106c the dot product of the gradients are computed. 107c 108c This implementation includes the LDA exchange part [3] of the 109c exchange functional as well as the gradient correction part. 110c 111c The original code was provide by Prof. Walter Thiel. 112c 113c [1] M. Filatov, and Walter Thiel, 114c "A nonlocal correlation energy density functional from a 115c Coulomb hole model", 116c Int.J.Quant.Chem. 62 (1997) 603-616. 117c 118c [2] M. Filatov, and Walter Thiel, 119c "A new gradient-corrected exchange-correlation 120c density functional", 121c Mol.Phys. 91 (1997) 847-859. 122c 123c [3] P.A.M. Dirac, Proceedings of the Cambridge Philosophical 124c Society, Vol. 26 (1930) 376. 125c 126c Parameters: 127c 128c r the total electron density 129c g the dot product of the total density gradient with itself 130c f On return the functional value 131c dfdra On return the derivative of f with respect to the 132c alpha-electron density 133c dfdgaa On return the derivative of f with respect to the dot 134c product of the alpha-electron density gradient with itself 135c 136 implicit none 137c 138 double precision r, g 139 double precision fc ,dfdrac ,dfdrb ,dfdgaac ,dfdgbb 140 double precision fx,dfdrax, dfdgaax 141c 142 double precision rhoa, rhob, rhoa13, rhob13 143 double precision gama, gamb 144c 145c... Parameters 146c 147 double precision r13,tol_rho 148 parameter (r13=1.0d0/3.0d0) 149c 150 rhoa = 0.5d0*r 151 rhoa13 = rhoa**r13 152 gama = 0.25d0*g 153 rhob = rhoa 154 rhob13 = rhoa13 155 gamb = gama 156 call FT97_ECFUN(rhoa,rhob,rhoa13,rhob13,gama,gamb, 157 + fc,dfdrac,dfdrb,dfdgaac,dfdgbb,tol_rho) 158 call FT97_EXFUN(rhoa,rhob,rhoa13,rhob13, 159 + gama,gamb,fx,.false.,.false.,tol_rho) 160 call FT97_EXGRAD(rhoa,rhoa13,gama,dfdrax,dfdgaax,.false.) 161 end 162c----------------------------------------------------------------------- 163 subroutine xc_uks_ft97(tol_rho,ra,rb,gaa,gbb, 164 , fc,dfdrac,dfdrbc,dfdgaac,dfdgbbc, 165 , fx,dfdrax,dfdrbx,dfdgaax,dfdgbbx) 166c 167c This subroutine calculates the Filatov-Thiel 97 168c exchange-correlation functional [1,2] for the spin polarised case. 169c This functional is taken to be the correlation functional [1] plus 170c the recommended variant B of the exchange functional [2]. 171c Also the derivatives with respect to the density components and 172c the dot product of the gradients are computed. 173c 174c This implementation includes the LDA exchange part [3] of the 175c exchange functional as well as the gradient correction part. 176c 177c The original code was provide by Prof. Walter Thiel. 178c 179c [1] M. Filatov, and Walter Thiel, 180c "A nonlocal correlation energy density functional from a 181c Coulomb hole model", 182c Int.J.Quant.Chem. 62 (1997) 603-616. 183c 184c [2] M. Filatov, and Walter Thiel, 185c "A new gradient-corrected exchange-correlation 186c density functional", 187c Mol.Phys. 91 (1997) 847-859. 188c 189c [3] P.A.M. Dirac, Proceedings of the Cambridge Philosophical 190c Society, Vol. 26 (1930) 376. 191c 192c Parameters: 193c 194c ra the alpha-electron density 195c rb the beta-electron density 196c gaa the dot product of the alpha density gradient with itself 197c gbb the dot product of the beta density gradient with itself 198c the beta density 199c f On return the functional value 200c dfdra On return the derivative of f with respect to ra 201c dfdrb On return the derivative of f with respect to rb 202c dfdgaa On return the derivative of f with respect to gaa 203c dfdgbb On return the derivative of f with respect to gbb 204c 205 implicit none 206c 207 double precision ra, rb, gaa, gbb 208 double precision fc ,dfdrac ,dfdrbc ,dfdgaac ,dfdgbbc 209 double precision fx,dfdrax,dfdrbx,dfdgaax,dfdgbbx 210c 211 double precision rhoa13, rhob13 212c 213 double precision r13,tol_rho 214 parameter (r13=1.0d0/3.0d0) 215c 216 rhoa13=0d0 217 if(ra.gt.tol_rho**2) rhoa13 = ra**r13 218 rhob13=0d0 219 if(rb.gt.tol_rho**2) rhob13 = rb**r13 220 call FT97_ECFUN(ra,rb,rhoa13,rhob13,gaa,gbb, 221 + fc,dfdrac,dfdrbc,dfdgaac,dfdgbbc,tol_rho) 222 call FT97_EXFUN(ra,rb,rhoa13,rhob13, 223 + gaa,gbb,fx,.true.,.false.,tol_rho) 224 if(ra.gt.tol_rho**2) 225 + call FT97_EXGRAD(ra,rhoa13,gaa,dfdrax,dfdgaax,.false.) 226 if(rb.gt.tol_rho**2) 227 + call FT97_EXGRAD(rb,rhob13,gbb,dfdrbx,dfdgbbx,.false.) 228 end 229c----------------------------------------------------------------------- 230c----------------------------------------------------------------------- 231 SUBROUTINE FT97_EXFUN (RHOA,RHOB,RHOA13,RHOB13,AMA,AMB,FT97,UHF, 232 + VARIA,tol_rho) 233C * 234C VALUE OF THE FT97 EXCHANGE xc_ft97 AT A GIVEN POINT IN SPACE. 235C * 236C REFERENCE: M.FILATOV AND W.THIEL, MOL.PHYS. 91, 847 (1997). 237C * 238C ARGUMENT LIST. I=INPUT, O=OUTPUT. 239C RHOA DENSITY rho.alpha (I) 240C RHOB DENSITY rho.beta (I) 241C RHOA13 RHOA**(1.0/3.0), CUBIC ROOT OF rho.alpha (I) 242C RHOB13 RHOB**(1.0/3.0), CUBIC ROOT OF rho.beta (I) 243C AMA NORM OF THE GRADIENT OF RHOA WRT X,Y,Z SQUARED (I) 244C AMB NORM OF THE GRADIENT OF RHOB WRT X,Y,Z SQUARED (I) 245C FT97 VALUE OF EXCHANGE xc_ft97 (O) 246C UHF LOGICAL FLAG (.TRUE. FOR UHF) (I) 247C VARIA LOGICAL FLAG FOR CHOICE OF xc_ft97 (I) 248C =.TRUE. VARIANT A (I) 249C =.FALSE. VARIANT B (RECOMMENDED) (I) 250C * 251 IMPLICIT double precision (A-H,O-Z) 252 LOGICAL UHF,VARIA 253 double precision tol_rho 254C Variants A and B, see eqs.(15) and (16). 255 PARAMETER (BETAA=2.930000D-3) 256 PARAMETER (AX=9.474169D-4, BX=2.501149D03, BETA0=2.913644D-3) 257C BETA : Beta(sigma), eq.(16a) for variant B. 258 IF(VARIA) THEN 259 BETA = BETAA 260 ELSE 261 BETA = BETA0+AX*AMA/(BX*BX+AMA) 262 ENDIF 263C XALPHA : Reduced density gradient xi(sigma), below eq.(11). 264 XALPHA = (AMA/(RHOA*RHOA13)**2) 265C DENOMA : Denominator in f(xi), eqs. (15) and (16). 266 SINHIA = LOG(XALPHA+SQRT(XALPHA*XALPHA+1.0D0)) 267 DENOMA = SQRT(1.0D0+9.D0*(BETA**2)*XALPHA*(SINHIA**2)) 268C CORRAI : Nonlocal correction to exchange energy density ex. 269C CORRAI : See eqs.(9),(11),(15) and (16). 270C CORRAI : Part of integrand in eq.(11) beyond ex(LDA). 271 CORRAI =-RHOA*RHOA13*BETA*XALPHA/DENOMA 272 IF(UHF) THEN 273 IF(RHOB.LT.tol_rho) THEN 274 FT97 = CORRAI 275 ELSE 276 IF(.NOT.VARIA) BETA = BETA0+AX*AMB/(BX*BX+AMB) 277 XBETA = (AMB/(RHOB*RHOB13)**2) 278 SINHIB = LOG(XBETA+SQRT(XBETA*XBETA+1.0D0)) 279 DENOMB = SQRT(1.0D0+9.D0*BETA**2*XBETA*SINHIB**2) 280 CORRBI =-BETA*RHOB*RHOB13*XBETA/DENOMB 281 FT97 = CORRAI + CORRBI 282 ENDIF 283 ELSE 284 FT97 = 2.0D0*CORRAI 285 ENDIF 286 RETURN 287 END 288c----------------------------------------------------------------------- 289 SUBROUTINE FT97_EXGRAD (RHO,RHO13,AM,TERM1,U,VARIA) 290C * 291C DERIVATIVES OF THE FT97 EXCHANGE xc_ft97 WITH RESPECT TO THE 292C DENSITY AT A GIVEN POINT IN SPACE (EXCHANGE POTENTIAL). 293C * 294C REFERENCE: M.FILATOV AND W.THIEL, MOL.PHYS. 91, 847 (1997). 295C * 296C ARGUMENT LIST. I=INPUT, O=OUTPUT. 297C RHO DENSITY rho (I) 298C RHO13 RHO**(1.0/3.0), CUBIC ROOT OF rho (I) 299C AM NORM OF THE GRADIENT OF RHO WRT X,Y,Z SQUARED (I) 300C TERM1 DERIVATIVE d(Ex)/d(rho) (O) 301C U DERIVATIVE d(Ex)/d(grad(rho)^2) (O) 302C VARIA LOGICAL FLAG FOR CHOICE OF xc_ft97 (I) 303C =.TRUE. VARIANT A (I) 304C =.FALSE. VARIANT B (RECOMMENDED) (I) 305C 306 IMPLICIT double precision (A-H,O-Z) 307 LOGICAL VARIA 308C Variants A and B, see eqs.(15) and (16). 309 PARAMETER (BETAA=2.930000D-3) 310 PARAMETER (AX=9.474169D-4, BX=2.501149D03, BETA0=2.913644D-3) 311 PARAMETER (FT=4.0D0/3.0D0) 312C BETA : Beta(sigma), eq.(16a) for variant B. 313C DBDG : Derivative of BETA w.r.t. square of density gradient. 314 IF(VARIA) THEN 315 BETA = BETAA 316 dbdg=0d0 317 ELSE 318 BETA = BETA0+AX*AM/(BX*BX+AM) 319 DBDG = AX*BX*BX/(BX*BX+AM)**2 320 ENDIF 321 XL = 1.0D0/(RHO*RHO13) 322C ZT : Reduced density gradient xi(sigma), below eq.(11). 323 ZT = AM*(XL)**2 324C DENOM : Denominator in f(xi), eqs. (15) and (16). 325 SQZT = SQRT(ZT*ZT+1.D0) 326 SH = LOG(ZT+SQZT) 327 DEN = 1.D0+9.D0*BETA**2*ZT*SH**2 328 DENOM = SQRT(DEN) 329C F : Last factor in eqs.(15) and (16). 330 F = BETA*ZT/DENOM 331C FZ : Derivative of F with respect to ZT (see above). 332 FZ = -4.5D0*BETA**3*ZT*(2.D0*ZT*SH/SQRT(1.D0+ZT**2)+SH**2)/ 333 1 (DEN*DENOM) + BETA/DENOM 334C TERM1 : Derivative d(Ex)/d(rho). 335 TERM1 = FT*RHO13*(2.D0*ZT*FZ - F) 336C U : Derivative d(Ex)/d(grad(rho)^2), variants A and B. 337 IF(VARIA) THEN 338 U =-0.5d0*(FZ+FZ)*XL 339 ELSE 340C EB : Derivative of exchange energy with respect to BETA. 341 EB =-RHO*RHO13*ZT/(DEN*DENOM) 342 U =-(FZ*XL-EB*DBDG) 343 ENDIF 344 RETURN 345 END 346c----------------------------------------------------------------------- 347c----------------------------------------------------------------------- 348 SUBROUTINE FT97_ECFUN (RHOA,RHOB,RHOA13,RHOB13,AMA,AMB,EC,V1,V2, 349 1 V3,V4,tol_rho) 350C * 351C VALUE OF THE CORRELATION xc_ft97 (EC) AND ITS DERIVATIVES 352C (V1,V2,V3,V4) WITH REGARD TO THE DENSITY AT A GIVEN POINT 353C IN SPACE WITH CARTESIAN COORDINATES X,Y,Z. 354C * 355C REFERENCE: 356C M.FILATOV AND W.THIEL, 357C "A nonlocal correlation energy density functional from a Coulomb 358C hole model", INT.J.QUANTUM CHEM. Vol. 62, (1997) 603-616. 359C * 360C ARGUMENT LIST. I=INPUT, O=OUTPUT. 361C RHOA DENSITY rho.alpha (I) 362C RHOB DENSITY rho.beta (I) 363C RHOA13 RHOA**(1.0/3.0), CUBIC ROOT OF rho.alpha (I) 364C RHOB13 RHOB**(1.0/3.0), CUBIC ROOT OF rho.beta (I) 365C AMA NORM OF THE GRADIENT OF RHOA WRT X,Y,Z SQUARED (I) 366C AMB NORM OF THE GRADIENT OF RHOB WRT X,Y,Z SQUARED (I) 367C EC CORRELATION ENERGY (O) 368C V1 DERIVATIVE d(Ec)/d(rho.alpha) (I,O) 369C V2 DERIVATIVE d(Ec)/d(rho.beta) (I,O) 370C V3 DERIVATIVE d(Ec)/d(grad(rho.alpha)^2) (I,O) 371C V4 DERIVATIVE d(Ec)/d(grad(rho.beta)^2) (I,O) 372C UHF LOGICAL FLAG (.TRUE. FOR UHF) (I) 373C * 374 IMPLICIT double precision (A-H,O-Z) 375 EC = 0.D0 376 V1 = 0.D0 377 V2 = 0.D0 378 V3 = 0.D0 379 V4 = 0.D0 380 CALL FT97_VCOUR (RHOA,RHOB,RHOA13,RHOB13,AMA,AMB,EAA,EBB,EAB,EBA, 381 1 DEAADR,DEBBDR,DEABDR,DEBADR, 382 2 DEAADG,DEBBDG,DEABDG,DEBADG,tol_rho) 383C EC IS THE CONTRIBUTION TO THE CORRELATION ENERGY AT (X,Y,Z). 384C EC IS DEFINED IN EQUATION (4) OF THE REFERENCE. 385C Add parallel-spin contribution to the energy. 386 EC = EC+RHOA*EAA+RHOB*EBB 387C Add opposite-spin contribution to the energy. 388 EC = EC+RHOA*EAB+RHOB*EBA 389C Add parallel-spin contribution to the derivatives. 390 V1 = V1 +EAA+RHOA*DEAADR 391 V2 = V2 +EBB+RHOB*DEBBDR 392 V3 = V3 +RHOA*DEAADG 393 V4 = V4 +RHOB*DEBBDG 394C Add opposite-spin contribution to the derivatives. 395 V1 = V1 +EAB+RHOB*DEBADR 396 V2 = V2 +EBA+RHOA*DEABDR 397 V3 = V3 +RHOB*DEBADG 398 V4 = V4 +RHOA*DEABDG 399 RETURN 400 END 401c----------------------------------------------------------------------- 402 SUBROUTINE FT97_VCOUR (RHOA,RHOB,RHOA13,RHOB13,AMA,AMB, 403 1 EAA,EBB,EAB,EBA, 404 2 DEAADR,DEBBDR,DEABDR,DEBADR, 405 3 DEAADG,DEBBDG,DEABDG,DEBADG,tol_rho) 406C * 407C COMPUTATION OF TERMS IN THE CORRELATION xc_ft97. 408C * 409C REFERENCE: 410C M.FILATOV AND W.THIEL, INT.J.QUANTUM CHEM. 62, 603 (1997). 411C * 412C NOTATION FOR ARGUMENTS. 413C RHOA - rho(alpha), density, eq.(2) 414C RHOB - rho(beta ) 415C RHOA13 - rho(alpha)**1/3, cubic root of density 416C RHOB13 - rho(beta )**1/3 417C AMA - grad(rho(alpha))**2, gradient norm 418C AMB - grad(rho(beta ))**2 419C EAA - epsilon(alpha,alpha), correlation energy density, eq.(5) 420C EBB - epsilon(beta ,beta ) 421C EAB - epsilon(alpha,beta ) 422C EBA - epsilon(beta ,alpha) 423C DEDRAA - d(Ec)/d(rho.alpha) 424C DEDGAA - d(Ec)/d(grad(rho.alpha)^2) 425C * 426 IMPLICIT double precision (A-H,O-Z) 427 double precision KSSP0,KSS0,MUAB,MUBA,MUBB,MUAA, 428 * KSSPBIG,big,kss0big 429C NOTATION FOR CONSTANTS. 430C C0 = (1-ln2)/(2*Pi**2), see eq.(9) 431C C1 = 2*C0/3 , see eq.(13,(28),(33) 432C C2 = (3/(4*Pi))**1/3 , see eq.(8) 433C C3 = C2/3 434 PARAMETER (C0=0.1554534543483D-1) 435 PARAMETER (C1=0.2072712724644D-1) 436 PARAMETER (C2=0.6203504908994D00) 437 PARAMETER (C3=0.2067834969665D00) 438C OPTIMIZED PARAMETERS IN GRADIENT-CORRECTED CORRELATION xc_ft97. 439C See formulas for FSSP, eq.(45), and FSS, eq.(46). 440 PARAMETER (A1=1.622118767D0) 441 PARAMETER (A2=0.489958076D0) 442 PARAMETER (A3=1.379021941D0) 443 PARAMETER (A4=4.946281353D0) 444 PARAMETER (A5=3.600612059D0) 445C NUMERICAL CUTOFF FOR MUAA, MUAB, MUBA, MUBB, see eqs.(13) and (33). 446 PARAMETER (CUTOFF=1.0D7) 447 parameter(ksspbig=1.291551074D0 - 0.349064173D0,big=1d4) 448 parameter(KSS0BIG = 1.200801774D0 + 0.859614445D0- 449 - 0.812904345D0) 450C 451C *** DEFINITION OF FORMULA FUNCTIONS. 452C 453C R : Density parameter Rs (Wigner radius), see eq.(8). 454C KSSP0 : See eq.(39) for k(sigma,sigma'). 455 KSSP0(R) = 1.291551074D0 - 0.349064173D0*(1.D0- 456 . EXP(-0.083275880D0*R**0.8D0)) 457C KSS0 : See eq.(40) for k(sigma,sigma). 458 KSS0(R) = 1.200801774D0 + 0.859614445D0*(1.D0- 459 . EXP(-1.089338848*SQRT(R)))- 460 . 0.812904345D0*(1.D0-EXP(-0.655638823D0*R**0.4D0)) 461C 462C FSSP : See eq.(45) for F(sigma,sigma'). 463 FSSP(R,GR) = (1.D0+A1*GR+(A2*GR)**2)* 464 . EXP(-(A2*GR)**2)/SQRT(1.D0+A3*GR/R) 465C USSP1 : Derivative of ln(k(sigma,sigma')) with respect to Rs. 466 USSP1(R) = -4.65098019D-2*EXP(-0.083275880D0*R**0.8D0)*R**0.8D0 467C USSP2 : Derivative of ln(F(sigma,sigma')) with respect to Rs. 468 USSP2(R,GR) = A3*GR/(R+A3*GR) 469C WSSP : Derivative of ln(F(sigma,sigma')) w.r.t. Nabla(Rs)^2 470 WSSP(R,GR) = ((A1+A1)-A1*(A2+A2)**2*GR**2-(A2*(A2+A2))**2*GR**3) 471 . /(1.D0+A1*GR+(A2*GR)**2) - A3/(R+A3*GR) 472C 473C FSS : See eq.(46) for F(sigma,sigma). 474 FSS(R,GR)= (1.D0+(A4*GR)**2)*EXP(-(A4*GR)**2)/ 475 . SQRT(1.D0+A5*GR/R) 476C FACTOR : See eq.(34). 477 FACTOR(R)= EXP(-(R/(0.939016D0*SQRT(R)+1.733170D0*R))**2) 478C USS1 : Derivative of ln(k(sigma,sigma)) with respect to Rs. 479 USS1(R) =-4.26377318D-1*R**0.4D0*EXP(-0.655638823D0*R**0.4D0)+ 480 . 0.93641140924D0*SQRT(R)*EXP(-1.089338848*SQRT(R)) 481C USS2 : Derivative of ln(F(sigma,sigma)) with respect to Rs. 482 USS2(R,GR)= A5*GR/(R+A5*GR) 483C WSS : Derivative of ln(F(sigma,sigma)) w.r.t. Nabla(Rs)^2 484 WSS(R,GR) =-(A4*(A4+A4))**2*GR**3/(1.D0+(A4*GR)**2)-A5/(R+A5*GR) 485C YSS : Derivative of FACTOR (eq.(34)) with respect to Rs. 486 YSS(R) = 0.939016D0*SQRT(R)*R*R/(0.939016D0*SQRT(R)+ 487 . 1.733170D0*R)**3 488C 489C FF : Used in approximate normalization, eq.(15). 490 FF(T) = (3.D0+2.D0*(SQRT(T)+T))/ 491 . (3.D0+6.D0*(SQRT(T)+T)) 492C DFFDMU : Derivative of FF with respect to the argument T. 493 DFFDMU(T) =-(2.D0*(SQRT(T)+T+T))/ 494 . (3.D0*(1.D0+2.D0*(SQRT(T)+T))**2) 495C 496C *** RSA - Rs(alpha), GA - grad(Rs(alpha)) 497C 498 IF(RHOA.GT.tol_rho) THEN 499C Wigner radius Rs(alpha), eq.(8) 500 RSA = C2/RHOA13 501 CFT = C3/(RHOA*RHOA13) 502C GA : Nabla(Rs(alpha)) 503C GA2 : Square of Nabla(Rs(alpha)) 504 GA2 = CFT*CFT*AMA 505C mu(beta,alpha), eq.(13) 506C MUBA = C1*RSA/(KSSP0(RSA)*FSSP(RSA,GA2))**2 507 EBA = 0.D0 508 ECMUC0 = 0.D0 509 DEBADR = 0.D0 510 DEBADG = 0.D0 511 if((ga2*a2)**2.gt.big) then 512 DENOM = 0d0 513 else 514 if(rsa.gt.big) then 515 DENOM = (KSSPBIG*FSSP(RSA,GA2))**2 516 else 517 DENOM = (KSSP0(RSA)*FSSP(RSA,GA2))**2 518 endif 519 endif 520 IF(C1*RSA .LE. DENOM*CUTOFF) THEN 521 MUBA = C1*RSA/DENOM 522 EEI = FT97_EXPEI(-MUBA) 523 EEI1 = MUBA*EEI+1.D0 524C EBA : Correlation energy density, eq.(19) 525 EBA = C0*(EEI+2.D0*FF(MUBA)*EEI1) 526C ECMUC0 : Derivative of EBA w.r.t. MUBA 527 ECMUC0 = C0*(EEI1*(1.D0+2.D0*(DFFDMU(MUBA)+MUBA* 528 . FF(MUBA)))+2.D0*FF(MUBA)*MUBA*EEI) 529C d(Ec)/d(Nabla(rho)^2) = -CFT^2*(mu*dEcdmu)*Wss' 530 DEBADG = -WSSP(RSA,GA2)*ECMUC0*CFT*CFT 531C d(Ec)/d(rho) 532 DEBADR = -ECMUC0*(1.D0-USSP1(RSA)/KSSP0(RSA)- 533 . USSP2(RSA,GA2) - 534 . 8.D0*WSSP(RSA,GA2)*GA2) 535 . / (3.D0*RHOA) 536 ENDIF 537C mu(alpha,alpha), eq.(33) 538C MUAA = C1*RSA/(KSS0(RSA)*FSS(RSA,GA2))**2 539 FA = FACTOR(RSA) 540 EAA = 0.D0 541 ECMUC0 = 0.D0 542 DEAADR = 0.D0 543 DEAADG = 0.D0 544 if((ga2*a2)**2.gt.big) then 545 DENOM = 0d0 546 else 547 if(rsa.gt.big) then 548 DENOM = (KSS0big*FSS(RSA,GA2))**2 549 else 550 DENOM = (KSS0(RSA)*FSS(RSA,GA2))**2 551 endif 552 endif 553 IF(C1*RSA .LE. DENOM*CUTOFF) THEN 554 MUAA = C1*RSA/DENOM 555 EEI = FT97_EXPEI(-MUAA) 556 EEI1 = MUAA*EEI+1.D0 557C EAA : Correlation energy density, eqs.(19) and (31). 558 EAA = FA*C0*(EEI+2.D0*FF(MUAA)*EEI1) 559C ECMUC0 : Derivative of EAA w.r.t. MUAA 560 ECMUC0 = FA*C0*(EEI1*(1.D0+2.D0*(DFFDMU(MUAA)+MUAA* 561 . FF(MUAA)))+2.D0*FF(MUAA)*MUAA*EEI) 562C d(Ec)/d(rho) 563 DEAADR = (-ECMUC0*(1.D0-USS1(RSA)/KSS0(RSA)-USS2(RSA,GA2) - 564 . 8.D0*WSS(RSA,GA2)*GA2) + YSS(RSA)*EAA) 565 . / (3.D0*RHOA) 566C d(Ec)/d(Nabla(rho)^2) = -Cft^2*(mu*dEcdmu*Fa)*Wss 567 DEAADG = -WSS(RSA,GA2)*ECMUC0*CFT*CFT 568 ENDIF 569 ELSE 570 EAA = 0.D0 571 EBA = 0.D0 572 DEAADR = 0.D0 573 DEBADR = 0.D0 574 DEAADG = 0.D0 575 DEBADG = 0.D0 576 ENDIF 577C *** RSB - Rs(beta), GB - grad(Rs(beta)) 578 IF(RHOB.GT.tol_rho) THEN 579C Wigner radius Rs(beta), eq.(8) 580 RSB = C2/RHOB13 581 CFT = C3/(RHOB*RHOB13) 582C GB : Nabla(Rs(beta)) 583C GB2 : Square of Nabla(Rs(beta)) 584 GB2 = CFT*CFT*AMB 585C mu(alpha,beta), eq.(13) 586C MUAB = C1*RSB/(KSSP0(RSB)*FSSP(RSB,GB2))**2 587 EAB = 0.D0 588 ECMUC0 = 0.D0 589 DEABDR = 0.D0 590 DEABDG = 0.D0 591 DENOM = (KSSP0(RSB)*FSSP(RSB,GB2))**2 592 if((gb2*a2)**2.gt.big) then 593 DENOM = 0d0 594 else 595 if(rsb.gt.big) then 596 DENOM = (KSSPBIG*FSSP(RSB,GB2))**2 597 else 598 DENOM = (KSSP0(RSB)*FSSP(RSB,GB2))**2 599 endif 600 endif 601 IF(C1*RSB .LE. DENOM*CUTOFF) THEN 602 MUAB = C1*RSB/DENOM 603 EEI = FT97_EXPEI(-MUAB) 604 EEI1 = MUAB*EEI+1.D0 605C EAB : Correlation energy density, eqs.(19). 606 EAB = C0*(EEI+2.D0*FF(MUAB)*EEI1) 607C ECMUC0 : Derivative of EAB w.r.t. MUAB 608 ECMUC0 = C0*(EEI1*(1.D0+2.D0*(DFFDMU(MUAB)+MUAB* 609 . FF(MUAB)))+2.D0*FF(MUAB)*MUAB*EEI) 610C d(Ec)/d(Nabla(rho)^2) = -CFT^2*(mu*dEcdmu)*Wss' 611 DEABDG = -WSSP(RSB,GB2)*ECMUC0*CFT*CFT 612C d(Ec)/d(rho) 613 DEABDR = -ECMUC0*(1.D0-USSP1(RSB)/KSSP0(RSB)- 614 . USSP2(RSB,GB2) - 615 . 8.D0*WSSP(RSB,GB2)*GB2) 616 . / (3.D0*RHOB) 617 ENDIF 618C mu(beta,beta), eq.(33) 619C MUBB = C1*RSB/(KSS0(RSB)*FSS(RSB,GB2))**2 620 FA = FACTOR(RSB) 621 EBB = 0.D0 622 ECMUC0 = 0.D0 623 DEBBDR = 0.D0 624 DEBBDG = 0.D0 625 if((gb2*a2)**2.gt.big) then 626 DENOM = 0d0 627 else 628 if(rsb.gt.big) then 629 DENOM = (KSS0big*FSS(RSB,GB2))**2 630 else 631 DENOM = (KSS0(RSB)*FSS(RSB,GB2))**2 632 endif 633 endif 634 IF(C1*RSB .LE. DENOM*CUTOFF) THEN 635 MUBB = C1*RSB/DENOM 636 EEI = FT97_EXPEI(-MUBB) 637 EEI1 = MUBB*EEI+1.D0 638C EBB : Correlation energy density, eqs.(19) and (31). 639 EBB = FA*C0*(EEI+2.D0*FF(MUBB)*EEI1) 640C ECMUC0 : Derivative of EBB w.r.t. MUBB 641 ECMUC0 = FA*C0*(EEI1*(1.D0+2.D0*(DFFDMU(MUBB)+MUBB* 642 . FF(MUBB)))+2.D0*FF(MUBB)*MUBB*EEI) 643C d(Ec)/d(rho) 644 DEBBDR = (-ECMUC0*(1.D0-USS1(RSB)/KSS0(RSB) 645 . -USS2(RSB,GB2) - 646 . 8.D0*WSS(RSB,GB2)*GB2) + YSS(RSB)*EBB) 647 . / (3.D0*RHOB) 648C d(Ec)/d(Nabla(rho)^2) = -Cft^2*(mu*dEcdmu*Fa)*Wss 649 DEBBDG = -WSS(RSB,GB2)*ECMUC0*CFT*CFT 650 ENDIF 651 ELSE 652 EAB = 0.D0 653 EBB = 0.D0 654 DEABDR = 0.D0 655 DEBBDR = 0.D0 656 DEABDG = 0.D0 657 DEBBDG = 0.D0 658 ENDIF 659 RETURN 660 END 661c----------------------------------------------------------------------- 662 FUNCTION FT97_EXPEI(X) 663C 664C This function program computes approximate values for the 665C function exp(-x) * Ei(x), where Ei(x) is the exponential 666C integral, and x is real. 667C 668C Author: W. J. Cody 669C 670C Latest modification: January 12, 1988 671C 672 INTEGER INT 673 double precision FT97_EXPEI, X, RESULT 674 INT = 3 675 CALL FT97_CALCEI(X,RESULT,INT) 676 FT97_EXPEI = RESULT 677 RETURN 678 END 679c----------------------------------------------------------------------- 680 SUBROUTINE FT97_CALCEI(ARG,RESULT,INT) 681C 682C This Fortran 77 packet computes the exponential integrals Ei(x), 683C E1(x), and exp(-x)*Ei(x) for real arguments x where 684C 685C integral (from t=-infinity to t=x) (exp(t)/t), x > 0, 686C Ei(x) = 687C -integral (from t=-x to t=infinity) (exp(t)/t), x < 0, 688C 689C and where the first integral is a principal value integral. 690C The packet contains three function type subprograms: EI, EONE, 691C and EXPEI; and one subroutine type subprogram: CALCEI. The 692C calling statements for the primary entries are 693C 694C Y = EI(X), where X .NE. 0, 695C 696C Y = EONE(X), where X .GT. 0, 697C and 698C Y = EXPEI(X), where X .NE. 0, 699C 700C and where the entry points correspond to the functions Ei(x), 701C E1(x), and exp(-x)*Ei(x), respectively. The routine CALCEI 702C is intended for internal packet use only, all computations within 703C the packet being concentrated in this routine. The function 704C subprograms invoke CALCEI with the Fortran statement 705C CALL CALCEI(ARG,RESULT,INT) 706C where the parameter usage is as follows 707C 708C Function Parameters for CALCEI 709C Call ARG RESULT INT 710C 711C EI(X) X .NE. 0 Ei(X) 1 712C EONE(X) X .GT. 0 -Ei(-X) 2 713C EXPEI(X) X .NE. 0 exp(-X)*Ei(X) 3 714C 715C The main computation involves evaluation of rational Chebyshev 716C approximations published in Math. Comp. 22, 641-649 (1968), and 717C Math. Comp. 23, 289-303 (1969) by Cody and Thacher. This 718C transportable program is patterned after the machine-dependent 719C FUNPACK packet NATSEI, but cannot match that version for 720C efficiency or accuracy. This version uses rational functions 721C that theoretically approximate the exponential integrals to 722C at least 18 significant decimal digits. The accuracy achieved 723C depends on the arithmetic system, the compiler, the intrinsic 724C functions, and proper selection of the machine-dependent 725C constants. 726C 727C 728C******************************************************************* 729C******************************************************************* 730C 731C Explanation of machine-dependent constants 732C 733C beta = radix for the floating-point system. 734C minexp = smallest representable power of beta. 735C maxexp = smallest power of beta that overflows. 736C XBIG = largest argument acceptable to EONE; solution to 737C equation: 738C exp(-x)/x * (1 + 1/x) = beta ** minexp. 739C XINF = largest positive machine number; approximately 740C beta ** maxexp 741C XMAX = largest argument acceptable to EI; solution to 742C equation: exp(x)/x * (1 + 1/x) = beta ** maxexp. 743C 744C Approximate values for some important machines are: 745C 746C beta minexp maxexp 747C 748C CRAY-1 (S.P.) 2 -8193 8191 749C Cyber 180/185 750C under NOS (S.P.) 2 -975 1070 751C IEEE (IBM/XT, 752C SUN, etc.) (S.P.) 2 -126 128 753C IEEE (IBM/XT, 754C SUN, etc.) (D.P.) 2 -1022 1024 755C IBM 3033 (D.P.) 16 -65 63 756C VAX D-Format (D.P.) 2 -128 127 757C VAX G-Format (D.P.) 2 -1024 1023 758C 759C XBIG XINF XMAX 760C 761C CRAY-1 (S.P.) 5670.31 5.45E+2465 5686.21 762C Cyber 180/185 763C under NOS (S.P.) 669.31 1.26E+322 748.28 764C IEEE (IBM/XT, 765C SUN, etc.) (S.P.) 82.93 3.40E+38 93.24 766C IEEE (IBM/XT, 767C SUN, etc.) (D.P.) 701.84 1.79D+308 716.35 768C IBM 3033 (D.P.) 175.05 7.23D+75 179.85 769C VAX D-Format (D.P.) 84.30 1.70D+38 92.54 770C VAX G-Format (D.P.) 703.22 8.98D+307 715.66 771C 772C******************************************************************* 773C******************************************************************* 774C 775C Error returns 776C 777C The following table shows the types of error that may be 778C encountered in this routine and the function value supplied 779C in each case. 780C 781C Error Argument Function values for 782C Range EI EXPEI EONE 783C 784C UNDERFLOW (-)X .GT. XBIG 0 - 0 785C OVERFLOW X .GE. XMAX XINF - - 786C ILLEGAL X X = 0 -XINF -XINF XINF 787C ILLEGAL X X .LT. 0 - - USE ABS(X) 788C 789C Intrinsic functions required are: 790C 791C ABS, SQRT, EXP 792C 793C 794C Author: W. J. Cody 795C Mathematics abd Computer Science Division 796C Argonne National Laboratory 797C Argonne, IL 60439 798C 799C Latest modification: September 9, 1988 800C 801C---------------------------------------------------------------------- 802 INTEGER I,INT 803 double precision 804 1 A,ARG,B,C,D,EXP40,E,EI,F,FOUR,FOURTY,FRAC,HALF,ONE,P, 805 2 PLG,PX,P037,P1,P2,q,QLG,QX,Q1,Q2,R,RESULT,S,SIX,SUMP, 806 3 SUMQ,T,THREE,TWELVE,TWO,TWO4,W,X,XBIG,XINF,XMAX,XMX0, 807 4 X0,X01,X02,X11,Y,YSQ,ZERO 808 DIMENSION A(7),B(6),C(9),D(9),E(10),F(10),P(10),q(10),R(10), 809 1 S(9),P1(10),Q1(9),P2(10),Q2(9),PLG(4),QLG(4),PX(10),QX(10) 810C---------------------------------------------------------------------- 811C Mathematical constants 812C EXP40 = exp(40) 813C X0 = zero of Ei 814C X01/X11 + X02 = zero of Ei to extra precision 815C---------------------------------------------------------------------- 816 DATA ZERO,P037,HALF,ONE,TWO/0.0D0,0.037D0,0.5D0,1.0D0,2.0D0/, 817 1 THREE,FOUR,SIX,TWELVE,TWO4/3.0D0,4.0D0,6.0D0,12.D0,24.0D0/, 818 2 FOURTY,EXP40/40.0D0,2.3538526683701998541D17/, 819 3 X01,X11,X02/381.5D0,1024.0D0,-5.1182968633365538008D-5/, 820 4 X0/3.7250741078136663466D-1/ 821C---------------------------------------------------------------------- 822C Machine-dependent constants 823C---------------------------------------------------------------------- 824CS DATA XINF/3.40E+38/,XMAX/93.246E0/,XBIG/82.93E0/ 825 DATA XINF/1.79D+308/,XMAX/716.351D0/,XBIG/701.84D0/ 826C---------------------------------------------------------------------- 827C Coefficients for -1.0 <= X < 0.0 828C---------------------------------------------------------------------- 829 DATA A/1.1669552669734461083368D2, 2.1500672908092918123209D3, 830 1 1.5924175980637303639884D4, 8.9904972007457256553251D4, 831 2 1.5026059476436982420737D5,-1.4815102102575750838086D5, 832 3 5.0196785185439843791020D0/ 833 DATA B/4.0205465640027706061433D1, 7.5043163907103936624165D2, 834 1 8.1258035174768735759855D3, 5.2440529172056355429883D4, 835 2 1.8434070063353677359298D5, 2.5666493484897117319268D5/ 836C---------------------------------------------------------------------- 837C Coefficients for -4.0 <= X < -1.0 838C---------------------------------------------------------------------- 839 DATA C/3.828573121022477169108D-1, 1.107326627786831743809D+1, 840 1 7.246689782858597021199D+1, 1.700632978311516129328D+2, 841 2 1.698106763764238382705D+2, 7.633628843705946890896D+1, 842 3 1.487967702840464066613D+1, 9.999989642347613068437D-1, 843 4 1.737331760720576030932D-8/ 844 DATA D/8.258160008564488034698D-2, 4.344836335509282083360D+0, 845 1 4.662179610356861756812D+1, 1.775728186717289799677D+2, 846 2 2.953136335677908517423D+2, 2.342573504717625153053D+2, 847 3 9.021658450529372642314D+1, 1.587964570758947927903D+1, 848 4 1.000000000000000000000D+0/ 849C---------------------------------------------------------------------- 850C Coefficients for X < -4.0 851C---------------------------------------------------------------------- 852 DATA E/1.3276881505637444622987D+2,3.5846198743996904308695D+4, 853 1 1.7283375773777593926828D+5,2.6181454937205639647381D+5, 854 2 1.7503273087497081314708D+5,5.9346841538837119172356D+4, 855 3 1.0816852399095915622498D+4,1.0611777263550331766871D03, 856 4 5.2199632588522572481039D+1,9.9999999999999999087819D-1/ 857 DATA F/3.9147856245556345627078D+4,2.5989762083608489777411D+5, 858 1 5.5903756210022864003380D+5,5.4616842050691155735758D+5, 859 2 2.7858134710520842139357D+5,7.9231787945279043698718D+4, 860 3 1.2842808586627297365998D+4,1.1635769915320848035459D+3, 861 4 5.4199632588522559414924D+1,1.0D0/ 862C---------------------------------------------------------------------- 863C Coefficients for rational approximation to ln(x/a), |1-x/a| < .1 864C---------------------------------------------------------------------- 865 DATA PLG/-2.4562334077563243311D+01,2.3642701335621505212D+02, 866 1 -5.4989956895857911039D+02,3.5687548468071500413D+02/ 867 DATA QLG/-3.5553900764052419184D+01,1.9400230218539473193D+02, 868 1 -3.3442903192607538956D+02,1.7843774234035750207D+02/ 869C---------------------------------------------------------------------- 870C Coefficients for 0.0 < X < 6.0, 871C ratio of Chebyshev polynomials 872C---------------------------------------------------------------------- 873 DATA P/-1.2963702602474830028590D01,-1.2831220659262000678155D03, 874 1 -1.4287072500197005777376D04,-1.4299841572091610380064D06, 875 2 -3.1398660864247265862050D05,-3.5377809694431133484800D08, 876 3 3.1984354235237738511048D08,-2.5301823984599019348858D10, 877 4 1.2177698136199594677580D10,-2.0829040666802497120940D11/ 878 DATA q/ 7.6886718750000000000000D01,-5.5648470543369082846819D03, 879 1 1.9418469440759880361415D05,-4.2648434812177161405483D06, 880 2 6.4698830956576428587653D07,-7.0108568774215954065376D08, 881 3 5.4229617984472955011862D09,-2.8986272696554495342658D10, 882 4 9.8900934262481749439886D10,-8.9673749185755048616855D10/ 883C---------------------------------------------------------------------- 884C J-fraction coefficients for 6.0 <= X < 12.0 885C---------------------------------------------------------------------- 886 DATA R/-2.645677793077147237806D00,-2.378372882815725244124D00, 887 1 -2.421106956980653511550D01, 1.052976392459015155422D01, 888 2 1.945603779539281810439D01,-3.015761863840593359165D01, 889 3 1.120011024227297451523D01,-3.988850730390541057912D00, 890 4 9.565134591978630774217D00, 9.981193787537396413219D-1/ 891 DATA S/ 1.598517957704779356479D-4, 4.644185932583286942650D00, 892 1 3.697412299772985940785D02,-8.791401054875438925029D00, 893 2 7.608194509086645763123D02, 2.852397548119248700147D01, 894 3 4.731097187816050252967D02,-2.369210235636181001661D02, 895 4 1.249884822712447891440D00/ 896C---------------------------------------------------------------------- 897C J-fraction coefficients for 12.0 <= X < 24.0 898C---------------------------------------------------------------------- 899 DATA P1/-1.647721172463463140042D00,-1.860092121726437582253D01, 900 1 -1.000641913989284829961D01,-2.105740799548040450394D01, 901 2 -9.134835699998742552432D-1,-3.323612579343962284333D01, 902 3 2.495487730402059440626D01, 2.652575818452799819855D01, 903 4 -1.845086232391278674524D00, 9.999933106160568739091D-1/ 904 DATA Q1/ 9.792403599217290296840D01, 6.403800405352415551324D01, 905 1 5.994932325667407355255D01, 2.538819315630708031713D02, 906 2 4.429413178337928401161D01, 1.192832423968601006985D03, 907 3 1.991004470817742470726D02,-1.093556195391091143924D01, 908 4 1.001533852045342697818D00/ 909C---------------------------------------------------------------------- 910C J-fraction coefficients for X .GE. 24.0 911C---------------------------------------------------------------------- 912 DATA P2/ 1.75338801265465972390D02,-2.23127670777632409550D02, 913 1 -1.81949664929868906455D01,-2.79798528624305389340D01, 914 2 -7.63147701620253630855D00,-1.52856623636929636839D01, 915 3 -7.06810977895029358836D00,-5.00006640413131002475D00, 916 4 -3.00000000320981265753D00, 1.00000000000000485503D00/ 917 DATA Q2/ 3.97845977167414720840D04, 3.97277109100414518365D00, 918 1 1.37790390235747998793D02, 1.17179220502086455287D02, 919 2 7.04831847180424675988D01,-1.20187763547154743238D01, 920 3 -7.99243595776339741065D00,-2.99999894040324959612D00, 921 4 1.99999999999048104167D00/ 922C---------------------------------------------------------------------- 923 X = ARG 924 IF (X .EQ. ZERO) THEN 925 EI = -XINF 926 IF (INT .EQ. 2) EI = -EI 927 ELSE IF ((X .LT. ZERO) .OR. (INT .EQ. 2)) THEN 928C---------------------------------------------------------------------- 929C Calculate EI for negative argument or for E1. 930C---------------------------------------------------------------------- 931 Y = ABS(X) 932 IF (Y .LE. ONE) THEN 933 SUMP = A(7) * Y + A(1) 934 SUMQ = Y + B(1) 935 DO 110 I = 2, 6 936 SUMP = SUMP * Y + A(I) 937 SUMQ = SUMQ * Y + B(I) 938 110 CONTINUE 939 EI = LOG(Y) - SUMP / SUMQ 940 IF (INT .EQ. 3) EI = EI * EXP(Y) 941 ELSE IF (Y .LE. FOUR) THEN 942 W = ONE / Y 943 SUMP = C(1) 944 SUMQ = D(1) 945 DO 130 I = 2, 9 946 SUMP = SUMP * W + C(I) 947 SUMQ = SUMQ * W + D(I) 948 130 CONTINUE 949 EI = - SUMP / SUMQ 950 IF (INT .NE. 3) EI = EI * EXP(-Y) 951 ELSE 952 IF ((Y .GT. XBIG) .AND. (INT .LT. 3)) THEN 953 EI = ZERO 954 ELSE 955 W = ONE / Y 956 SUMP = E(1) 957 SUMQ = F(1) 958 DO 150 I = 2, 10 959 SUMP = SUMP * W + E(I) 960 SUMQ = SUMQ * W + F(I) 961 150 CONTINUE 962 EI = -W * (ONE - W * SUMP / SUMQ ) 963 IF (INT .NE. 3) EI = EI * EXP(-Y) 964 END IF 965 END IF 966 IF (INT .EQ. 2) EI = -EI 967 ELSE IF (X .LT. SIX) THEN 968C---------------------------------------------------------------------- 969C To improve conditioning, rational approximations are expressed 970C in terms of Chebyshev polynomials for 0 <= X < 6, and in 971C continued fraction form for larger X. 972C---------------------------------------------------------------------- 973 T = X + X 974 T = T / THREE - TWO 975 PX(1) = ZERO 976 QX(1) = ZERO 977 PX(2) = P(1) 978 QX(2) = q(1) 979 DO 210 I = 2, 9 980 PX(I+1) = T * PX(I) - PX(I-1) + P(I) 981 QX(I+1) = T * QX(I) - QX(I-1) + q(I) 982 210 CONTINUE 983 SUMP = HALF * T * PX(10) - PX(9) + P(10) 984 SUMQ = HALF * T * QX(10) - QX(9) + q(10) 985 FRAC = SUMP / SUMQ 986 XMX0 = (X - X01/X11) - X02 987 IF (ABS(XMX0) .GE. P037) THEN 988 EI = LOG(X/X0) + XMX0 * FRAC 989 IF (INT .EQ. 3) EI = EXP(-X) * EI 990 ELSE 991C---------------------------------------------------------------------- 992C Special approximation to ln(X/X0) for X close to X0 993C---------------------------------------------------------------------- 994 Y = XMX0 / (X + X0) 995 YSQ = Y*Y 996 SUMP = PLG(1) 997 SUMQ = YSQ + QLG(1) 998 DO 220 I = 2, 4 999 SUMP = SUMP*YSQ + PLG(I) 1000 SUMQ = SUMQ*YSQ + QLG(I) 1001 220 CONTINUE 1002 EI = (SUMP / (SUMQ*(X+X0)) + FRAC) * XMX0 1003 IF (INT .EQ. 3) EI = EXP(-X) * EI 1004 END IF 1005 ELSE IF (X .LT. TWELVE) THEN 1006 FRAC = ZERO 1007 DO 230 I = 1, 9 1008 FRAC = S(I) / (R(I) + X + FRAC) 1009 230 CONTINUE 1010 EI = (R(10) + FRAC) / X 1011 IF (INT .NE. 3) EI = EI * EXP(X) 1012 ELSE IF (X .LE. TWO4) THEN 1013 FRAC = ZERO 1014 DO 240 I = 1, 9 1015 FRAC = Q1(I) / (P1(I) + X + FRAC) 1016 240 CONTINUE 1017 EI = (P1(10) + FRAC) / X 1018 IF (INT .NE. 3) EI = EI * EXP(X) 1019 ELSE 1020 IF ((X .GE. XMAX) .AND. (INT .LT. 3)) THEN 1021 EI = XINF 1022 ELSE 1023 Y = ONE / X 1024 FRAC = ZERO 1025 DO 250 I = 1, 9 1026 FRAC = Q2(I) / (P2(I) + X + FRAC) 1027 250 CONTINUE 1028 FRAC = P2(10) + FRAC 1029 EI = Y + Y * Y * FRAC 1030 IF (INT .NE. 3) THEN 1031 IF (X .LE. XMAX-TWO4) THEN 1032 EI = EI * EXP(X) 1033 ELSE 1034C---------------------------------------------------------------------- 1035C Calculation reformulated to avoid premature overflow 1036C---------------------------------------------------------------------- 1037 EI = (EI * EXP(X-FOURTY)) * EXP40 1038 END IF 1039 END IF 1040 END IF 1041 END IF 1042 RESULT = EI 1043 RETURN 1044 END 1045c----------------------------------------------------------------------- 1046