1c 2c$Id$ 3c 4#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 5C> \ingroup nwxc 6C> @{ 7C> 8C> \file nwxc_c_ft97.F 9C> The Filatov, Thiel correlation functional 10C> 11C> @} 12#endif 13C> 14C> \ingroup nwxc_priv 15C> @{ 16C> 17C> \brief Evaluate the Filatov-Thiel correlation functional 18C> 19C> This subroutine calculates the Filatov-Thiel 97 correlation 20C> functional [1,2]. 21C> Also the derivatives with respect to the density components and 22C> the dot product of the gradients are computed. 23C> 24C> This implementation includes the LDA exchange part [3] of the 25C> exchange functional as well as the gradient correction part. 26C> 27C> The original code was provided by Prof. Walter Thiel. 28C> 29C> ### References ### 30C> 31C> [1] M. Filatov, W. Thiel, 32C> "A nonlocal correlation energy density functional from a 33C> Coulomb hole model", Int. J. Quant. Chem. <b>62</b> (1997) 34C> 603-616, DOI: 35C> <a href="https://doi.org/10.1002/(SICI)1097-461X(1997)62:6<603::AID-QUA4>3.0.CO;2-%23"> 36C> 10.1002/(SICI)1097-461X(1997)62:6<603::AID-QUA4>3.0.CO;2-#</a>. 37C> 38C> [2] M. Filatov, W. Thiel, 39C> "A new gradient-corrected exchange-correlation 40C> density functional", Mol. Phys. <b>91</b> (1997) 847-859, DOI: 41C> <a href="https://doi.org/10.1080/002689797170950"> 42C> 10.1080/002689797170950</a>. 43C> 44C> [3] P.A.M. Dirac, "Note on exchange phenomena in the Thomas atom", 45C> Math. Proc. Cambridge Philos. Soc. <b>26</b> (1930) 376-385, DOI: 46C> <a href="https://doi.org/10.1017/S0305004100016108"> 47C> 10.1017/S0305004100016108</a>. 48C> 49#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 50#if defined(NWAD_PRINT) 51 Subroutine nwxc_c_ft97_p(tol_rho, ipol, nq, wght, rho, rgamma, 52 & func) 53#else 54 Subroutine nwxc_c_ft97(tol_rho, ipol, nq, wght, rho, rgamma, 55 & func) 56#endif 57#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 58 Subroutine nwxc_c_ft97_d2(tol_rho, ipol, nq, wght, rho, rgamma, 59 & func) 60#else 61 Subroutine nwxc_c_ft97_d3(tol_rho, ipol, nq, wght, rho, rgamma, 62 & func) 63#endif 64#include "nwad.fh" 65 implicit none 66c 67#include "intf_nwxc_rks_c_ft97.fh" 68#include "intf_nwxc_uks_c_ft97.fh" 69c 70#include "nwxc_param.fh" 71c 72c Input and other parameters 73c 74 double precision tol_rho !< [Input] The lower limit on the density 75 integer ipol !< [Input] The number of spin channels 76 integer nq !< [Input] The number of points 77 double precision wght !< [Input] The weight of the functional 78c 79c Charge Density 80c 81 type(nwad_dble)::rho(nq,*) !< [Input] The density 82c 83c Charge Density Gradient 84c 85 type(nwad_dble)::rgamma(nq,*) !< [Input] The norm of the density gradients 86c 87c Sampling Matrices for the XC Potential 88c 89 type(nwad_dble)::func(nq) !< [Output] The value of the functional 90c double precision Amat(nq,*) !< [Output] The derivative wrt rho 91c double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 92 integer n 93 type(nwad_dble)::gammaval 94c to hcth 95 type(nwad_dble)::rhoa 96 type(nwad_dble)::rhob 97 type(nwad_dble)::za 98 type(nwad_dble)::zb 99 double precision dfdza,dfdzb 100c 101 type(nwad_dble)::fc_ft97 102 double precision dfdrac,dfdgaac,dfdrbc,dfdgbbc 103c 104 if(ipol.eq.1) then 105 do n=1,nq 106 if(rho(n,R_T).gt.tol_rho) then 107 gammaval=rgamma(n,G_TT) 108c gammaval=delrho(n,1,1)*delrho(n,1,1) + 109c & delrho(n,2,1)*delrho(n,2,1) + 110c & delrho(n,3,1)*delrho(n,3,1) 111#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 112#if defined(NWAD_PRINT) 113 call nwxc_rks_c_ft97_p(rho(n,R_T),gammaval, 114 * fc_ft97,dfdrac,dfdgaac,tol_rho) 115#else 116 call nwxc_rks_c_ft97(rho(n,R_T),gammaval, 117 * fc_ft97,dfdrac,dfdgaac,tol_rho) 118#endif 119#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 120 call nwxc_rks_c_ft97_d2(rho(n,R_T),gammaval, 121 * fc_ft97,dfdrac,dfdgaac,tol_rho) 122#else 123 call nwxc_rks_c_ft97_d3(rho(n,R_T),gammaval, 124 * fc_ft97,dfdrac,dfdgaac,tol_rho) 125#endif 126 func(n)=func(n)+fc_ft97*wght 127c Amat(n,D1_RA) = Amat(n,D1_RA)+dfdrac*wght 128c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dfdgaac*wght 129 endif 130 enddo 131 else 132 do n=1,nq 133 rhoa=0.0d0 134 rhob=0.0d0 135 za =0.0d0 136 zb =0.0d0 137 if(rho(n,R_A)+rho(n,R_B).gt.tol_rho) then 138 if (rho(n,R_A).gt.0.5d0*tol_rho) then 139 rhoa=rho(n,R_A) 140 za=rgamma(n,G_AA) 141 endif 142 if (rho(n,R_B).gt.0.5d0*tol_rho) then 143 rhob=rho(n,R_B) 144 zb=rgamma(n,G_BB) 145 endif 146c za=delrho(n,1,1)*delrho(n,1,1) + 147c & delrho(n,2,1)*delrho(n,2,1) + 148c & delrho(n,3,1)*delrho(n,3,1) 149c zb=delrho(n,1,2)*delrho(n,1,2) + 150c & delrho(n,2,2)*delrho(n,2,2) + 151c & delrho(n,3,2)*delrho(n,3,2) 152#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 153#if defined(NWAD_PRINT) 154 call nwxc_uks_c_ft97_p(tol_rho, 155 ( rhoa,rhob,za,zb, 156 * fc_ft97,dfdrac,dfdrbc,dfdgaac,dfdgbbc) 157#else 158 call nwxc_uks_c_ft97(tol_rho, 159 ( rhoa,rhob,za,zb, 160 * fc_ft97,dfdrac,dfdrbc,dfdgaac,dfdgbbc) 161#endif 162#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 163 call nwxc_uks_c_ft97_d2(tol_rho, 164 ( rhoa,rhob,za,zb, 165 * fc_ft97,dfdrac,dfdrbc,dfdgaac,dfdgbbc) 166#else 167 call nwxc_uks_c_ft97_d3(tol_rho, 168 ( rhoa,rhob,za,zb, 169 * fc_ft97,dfdrac,dfdrbc,dfdgaac,dfdgbbc) 170#endif 171 func(n)=func(n)+fc_ft97*wght 172c Amat(n,D1_RA) = Amat(n,D1_RA)+dfdrac*wght 173c Amat(n,D1_RB) = Amat(n,D1_RB)+dfdrbc*wght 174c dfdza=dfdgaac*wght 175c dfdzb=dfdgbbc*wght 176c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dfdza 177c Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + dfdzb 178 endif 179 enddo 180 endif 181 return 182 end 183cDFT functional repository: xc_ft97 fortran77 source 184CDFT repository Quantum Chemistry Group 185C 186C CCLRC Density functional repository Copyright notice 187C This database was prepared as a result of work undertaken by CCLRC. 188C Users may view, print and download the content for personal use only 189C and the content must not be used for any commercial purpose without CCLRC 190C prior written permission 191C 192 193c----------------------------------------------------------------------- 194#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 195#if defined(NWAD_PRINT) 196 subroutine nwxc_rks_c_ft97_p(r,g,fc,dfdrac,dfdgaac,tol_rho) 197#else 198 subroutine nwxc_rks_c_ft97(r,g,fc,dfdrac,dfdgaac,tol_rho) 199#endif 200#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 201 subroutine nwxc_rks_c_ft97_d2(r,g,fc,dfdrac,dfdgaac,tol_rho) 202#else 203 subroutine nwxc_rks_c_ft97_d3(r,g,fc,dfdrac,dfdgaac,tol_rho) 204#endif 205c 206c This subroutine calculates the Filatov-Thiel 97 207c exchange-correlation functional [1,2] for the closed shell case. 208c This functional is taken to be the correlation functional [1] plus 209c the recommended variant B of the exchange functional [2]. 210c Also the derivatives with respect to the density components and 211c the dot product of the gradients are computed. 212c 213c This implementation includes the LDA exchange part [3] of the 214c exchange functional as well as the gradient correction part. 215c 216c The original code was provide by Prof. Walter Thiel. 217c 218c [1] M. Filatov, and Walter Thiel, 219c "A nonlocal correlation energy density functional from a 220c Coulomb hole model", 221c Int.J.Quant.Chem. 62 (1997) 603-616. 222c 223c [2] M. Filatov, and Walter Thiel, 224c "A new gradient-corrected exchange-correlation 225c density functional", 226c Mol.Phys. 91 (1997) 847-859. 227c 228c [3] P.A.M. Dirac, Proceedings of the Cambridge Philosophical 229c Society, Vol. 26 (1930) 376. 230c 231c Parameters: 232c 233c r the total electron density 234c g the dot product of the total density gradient with itself 235c f On return the functional value 236c dfdra On return the derivative of f with respect to the 237c alpha-electron density 238c dfdgaa On return the derivative of f with respect to the dot 239c product of the alpha-electron density gradient with itself 240c 241#include "nwad.fh" 242 implicit none 243c 244#include "intf_nwxc_ft97_ecfun.fh" 245c 246 type(nwad_dble)::r, g 247 type(nwad_dble)::fc 248 double precision dfdrac ,dfdrb ,dfdgaac ,dfdgbb 249c 250 type(nwad_dble)::rhoa, rhob, rhoa13, rhob13 251 type(nwad_dble)::gama, gamb 252c 253c... Parameters 254c 255 double precision r13,tol_rho 256 parameter (r13=1.0d0/3.0d0) 257c 258 rhoa = 0.5d0*r 259 rhoa13 = rhoa**r13 260 gama = 0.25d0*g 261 rhob = rhoa 262 rhob13 = rhoa13 263 gamb = gama 264#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 265#if defined(NWAD_PRINT) 266 call nwxc_FT97_ECFUN_p(rhoa,rhob,rhoa13,rhob13,gama,gamb, 267 + fc,dfdrac,dfdrb,dfdgaac,dfdgbb,tol_rho) 268#else 269 call nwxc_FT97_ECFUN(rhoa,rhob,rhoa13,rhob13,gama,gamb, 270 + fc,dfdrac,dfdrb,dfdgaac,dfdgbb,tol_rho) 271#endif 272#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 273 call nwxc_FT97_ECFUN_d2(rhoa,rhob,rhoa13,rhob13,gama,gamb, 274 + fc,dfdrac,dfdrb,dfdgaac,dfdgbb,tol_rho) 275#else 276 call nwxc_FT97_ECFUN_d3(rhoa,rhob,rhoa13,rhob13,gama,gamb, 277 + fc,dfdrac,dfdrb,dfdgaac,dfdgbb,tol_rho) 278#endif 279c call FT97_EXFUN(rhoa,rhob,rhoa13,rhob13, 280c + gama,gamb,fx,.false.,.false.,tol_rho) 281c call FT97_EXGRAD(rhoa,rhoa13,gama,dfdrax,dfdgaax,.false.) 282 end 283c----------------------------------------------------------------------- 284#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 285#if defined(NWAD_PRINT) 286 subroutine nwxc_uks_c_ft97_p(tol_rho,ra,rb,gaa,gbb, 287 , fc,dfdrac,dfdrbc,dfdgaac,dfdgbbc) 288#else 289 subroutine nwxc_uks_c_ft97(tol_rho,ra,rb,gaa,gbb, 290 , fc,dfdrac,dfdrbc,dfdgaac,dfdgbbc) 291#endif 292#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 293 subroutine nwxc_uks_c_ft97_d2(tol_rho,ra,rb,gaa,gbb, 294 , fc,dfdrac,dfdrbc,dfdgaac,dfdgbbc) 295#else 296 subroutine nwxc_uks_c_ft97_d3(tol_rho,ra,rb,gaa,gbb, 297 , fc,dfdrac,dfdrbc,dfdgaac,dfdgbbc) 298#endif 299c 300c This subroutine calculates the Filatov-Thiel 97 301c exchange-correlation functional [1,2] for the spin polarised case. 302c This functional is taken to be the correlation functional [1] plus 303c the recommended variant B of the exchange functional [2]. 304c Also the derivatives with respect to the density components and 305c the dot product of the gradients are computed. 306c 307c This implementation includes the LDA exchange part [3] of the 308c exchange functional as well as the gradient correction part. 309c 310c The original code was provide by Prof. Walter Thiel. 311c 312c [1] M. Filatov, and Walter Thiel, 313c "A nonlocal correlation energy density functional from a 314c Coulomb hole model", 315c Int.J.Quant.Chem. 62 (1997) 603-616. 316c 317c [2] M. Filatov, and Walter Thiel, 318c "A new gradient-corrected exchange-correlation 319c density functional", 320c Mol.Phys. 91 (1997) 847-859. 321c 322c [3] P.A.M. Dirac, Proceedings of the Cambridge Philosophical 323c Society, Vol. 26 (1930) 376. 324c 325c Parameters: 326c 327c ra the alpha-electron density 328c rb the beta-electron density 329c gaa the dot product of the alpha density gradient with itself 330c gbb the dot product of the beta density gradient with itself 331c the beta density 332c f On return the functional value 333c dfdra On return the derivative of f with respect to ra 334c dfdrb On return the derivative of f with respect to rb 335c dfdgaa On return the derivative of f with respect to gaa 336c dfdgbb On return the derivative of f with respect to gbb 337c 338#include "nwad.fh" 339 implicit none 340c 341#include "intf_nwxc_ft97_ecfun.fh" 342c 343 type(nwad_dble)::ra, rb, gaa, gbb 344 type(nwad_dble)::fc 345 type(nwad_dble)::fx 346 double precision dfdrac ,dfdrbc ,dfdgaac ,dfdgbbc 347 double precision dfdrax,dfdrbx,dfdgaax,dfdgbbx 348c 349 type(nwad_dble)::rhoa13, rhob13 350c 351 double precision r13,tol_rho 352 parameter (r13=1.0d0/3.0d0) 353c 354 rhoa13=0d0 355 if(ra.gt.tol_rho**2) rhoa13 = ra**r13 356 rhob13=0d0 357 if(rb.gt.tol_rho**2) rhob13 = rb**r13 358#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 359#if defined(NWAD_PRINT) 360 call nwxc_FT97_ECFUN_p(ra,rb,rhoa13,rhob13,gaa,gbb, 361 + fc,dfdrac,dfdrbc,dfdgaac,dfdgbbc,tol_rho) 362#else 363 call nwxc_FT97_ECFUN(ra,rb,rhoa13,rhob13,gaa,gbb, 364 + fc,dfdrac,dfdrbc,dfdgaac,dfdgbbc,tol_rho) 365#endif 366#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 367 call nwxc_FT97_ECFUN_d2(ra,rb,rhoa13,rhob13,gaa,gbb, 368 + fc,dfdrac,dfdrbc,dfdgaac,dfdgbbc,tol_rho) 369#else 370 call nwxc_FT97_ECFUN_d3(ra,rb,rhoa13,rhob13,gaa,gbb, 371 + fc,dfdrac,dfdrbc,dfdgaac,dfdgbbc,tol_rho) 372#endif 373c call FT97_EXFUN(ra,rb,rhoa13,rhob13, 374c + gaa,gbb,fx,.true.,.false.,tol_rho) 375c if(ra.gt.tol_rho**2) 376c + call FT97_EXGRAD(ra,rhoa13,gaa,dfdrax,dfdgaax,.false.) 377c if(rb.gt.tol_rho**2) 378c + call FT97_EXGRAD(rb,rhob13,gbb,dfdrbx,dfdgbbx,.false.) 379 end 380c----------------------------------------------------------------------- 381c----------------------------------------------------------------------- 382#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 383#if defined(NWAD_PRINT) 384 SUBROUTINE nwxc_FT97_ECFUN_p(RHOA,RHOB,RHOA13,RHOB13,AMA,AMB,EC, 385 1 V1,V2,V3,V4,tol_rho) 386#else 387 SUBROUTINE nwxc_FT97_ECFUN(RHOA,RHOB,RHOA13,RHOB13,AMA,AMB,EC, 388 1 V1,V2,V3,V4,tol_rho) 389#endif 390#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 391 SUBROUTINE nwxc_FT97_ECFUN_d2(RHOA,RHOB,RHOA13,RHOB13,AMA,AMB,EC, 392 1 V1,V2,V3,V4,tol_rho) 393#else 394 SUBROUTINE nwxc_FT97_ECFUN_d3(RHOA,RHOB,RHOA13,RHOB13,AMA,AMB,EC, 395 1 V1,V2,V3,V4,tol_rho) 396#endif 397C * 398C VALUE OF THE CORRELATION xc_ft97 (EC) AND ITS DERIVATIVES 399C (V1,V2,V3,V4) WITH REGARD TO THE DENSITY AT A GIVEN POINT 400C IN SPACE WITH CARTESIAN COORDINATES X,Y,Z. 401C * 402C REFERENCE: 403C M.FILATOV AND W.THIEL, 404C "A nonlocal correlation energy density functional from a Coulomb 405C hole model", INT.J.QUANTUM CHEM. Vol. 62, (1997) 603-616. 406C * 407C ARGUMENT LIST. I=INPUT, O=OUTPUT. 408C RHOA DENSITY rho.alpha (I) 409C RHOB DENSITY rho.beta (I) 410 411C RHOA13 RHOA**(1.0/3.0), CUBIC ROOT OF rho.alpha (I) 412C RHOB13 RHOB**(1.0/3.0), CUBIC ROOT OF rho.beta (I) 413C AMA NORM OF THE GRADIENT OF RHOA WRT X,Y,Z SQUARED (I) 414C AMB NORM OF THE GRADIENT OF RHOB WRT X,Y,Z SQUARED (I) 415C EC CORRELATION ENERGY (O) 416C V1 DERIVATIVE d(Ec)/d(rho.alpha) (I,O) 417C V2 DERIVATIVE d(Ec)/d(rho.beta) (I,O) 418C V3 DERIVATIVE d(Ec)/d(grad(rho.alpha)^2) (I,O) 419C V4 DERIVATIVE d(Ec)/d(grad(rho.beta)^2) (I,O) 420C UHF LOGICAL FLAG (.TRUE. FOR UHF) (I) 421C * 422#include "nwad.fh" 423 IMPLICIT double precision (A-H,O-Z) 424c 425#include "intf_nwxc_ft97_vcour.fh" 426c 427 type(nwad_dble)::RHOA,RHOB,RHOA13,RHOB13,AMA,AMB,EAA,EBB, 428 1 EAB,EBA,EC 429 EC = 0.D0 430 V1 = 0.D0 431 V2 = 0.D0 432 V3 = 0.D0 433 V4 = 0.D0 434#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 435#if defined(NWAD_PRINT) 436 CALL nwxc_FT97_VCOUR_p(RHOA,RHOB,RHOA13,RHOB13,AMA,AMB,EAA,EBB, 437 1 EAB,EBA, 438 2 DEAADR,DEBBDR,DEABDR,DEBADR, 439 3 DEAADG,DEBBDG,DEABDG,DEBADG,tol_rho) 440#else 441 CALL nwxc_FT97_VCOUR(RHOA,RHOB,RHOA13,RHOB13,AMA,AMB,EAA,EBB, 442 1 EAB,EBA, 443 2 DEAADR,DEBBDR,DEABDR,DEBADR, 444 3 DEAADG,DEBBDG,DEABDG,DEBADG,tol_rho) 445#endif 446#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 447 CALL nwxc_FT97_VCOUR_d2(RHOA,RHOB,RHOA13,RHOB13,AMA,AMB,EAA,EBB, 448 1 EAB,EBA, 449 2 DEAADR,DEBBDR,DEABDR,DEBADR, 450 3 DEAADG,DEBBDG,DEABDG,DEBADG,tol_rho) 451#else 452 CALL nwxc_FT97_VCOUR_d3(RHOA,RHOB,RHOA13,RHOB13,AMA,AMB,EAA,EBB, 453 1 EAB,EBA, 454 2 DEAADR,DEBBDR,DEABDR,DEBADR, 455 3 DEAADG,DEBBDG,DEABDG,DEBADG,tol_rho) 456#endif 457C EC IS THE CONTRIBUTION TO THE CORRELATION ENERGY AT (X,Y,Z). 458C EC IS DEFINED IN EQUATION (4) OF THE REFERENCE. 459C Add parallel-spin contribution to the energy. 460 EC = EC+RHOA*EAA+RHOB*EBB 461C Add opposite-spin contribution to the energy. 462 EC = EC+RHOA*EAB+RHOB*EBA 463C Add parallel-spin contribution to the derivatives. 464c V1 = V1 +EAA+RHOA*DEAADR 465c V2 = V2 +EBB+RHOB*DEBBDR 466c V3 = V3 +RHOA*DEAADG 467c V4 = V4 +RHOB*DEBBDG 468C Add opposite-spin contribution to the derivatives. 469c V1 = V1 +EAB+RHOB*DEBADR 470c V2 = V2 +EBA+RHOA*DEABDR 471c V3 = V3 +RHOB*DEBADG 472c V4 = V4 +RHOA*DEABDG 473 RETURN 474 END 475c----------------------------------------------------------------------- 476#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 477#if defined(NWAD_PRINT) 478 SUBROUTINE nwxc_FT97_VCOUR_p(RHOA,RHOB,RHOA13,RHOB13,AMA,AMB, 479 1 EAA,EBB,EAB,EBA, 480 2 DEAADR,DEBBDR,DEABDR,DEBADR, 481 3 DEAADG,DEBBDG,DEABDG,DEBADG,tol_rho) 482#else 483 SUBROUTINE nwxc_FT97_VCOUR (RHOA,RHOB,RHOA13,RHOB13,AMA,AMB, 484 1 EAA,EBB,EAB,EBA, 485 2 DEAADR,DEBBDR,DEABDR,DEBADR, 486 3 DEAADG,DEBBDG,DEABDG,DEBADG,tol_rho) 487#endif 488#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 489 SUBROUTINE nwxc_FT97_VCOUR_d2(RHOA,RHOB,RHOA13,RHOB13,AMA,AMB, 490 1 EAA,EBB,EAB,EBA, 491 2 DEAADR,DEBBDR,DEABDR,DEBADR, 492 3 DEAADG,DEBBDG,DEABDG,DEBADG,tol_rho) 493#else 494 SUBROUTINE nwxc_FT97_VCOUR_d3(RHOA,RHOB,RHOA13,RHOB13,AMA,AMB, 495 1 EAA,EBB,EAB,EBA, 496 2 DEAADR,DEBBDR,DEABDR,DEBADR, 497 3 DEAADG,DEBBDG,DEABDG,DEBADG,tol_rho) 498#endif 499C * 500C COMPUTATION OF TERMS IN THE CORRELATION xc_ft97. 501C * 502C REFERENCE: 503C M.FILATOV AND W.THIEL, INT.J.QUANTUM CHEM. 62, 603 (1997). 504C * 505C NOTATION FOR ARGUMENTS. 506C RHOA - rho(alpha), density, eq.(2) 507C RHOB - rho(beta ) 508C RHOA13 - rho(alpha)**1/3, cubic root of density 509C RHOB13 - rho(beta )**1/3 510C AMA - grad(rho(alpha))**2, gradient norm 511C AMB - grad(rho(beta ))**2 512C EAA - epsilon(alpha,alpha), correlation energy density, eq.(5) 513C EBB - epsilon(beta ,beta ) 514C EAB - epsilon(alpha,beta ) 515C EBA - epsilon(beta ,alpha) 516C DEDRAA - d(Ec)/d(rho.alpha) 517C DEDGAA - d(Ec)/d(grad(rho.alpha)^2) 518C * 519#include "nwad.fh" 520 IMPLICIT double precision (A-H,O-Z) 521c 522#include "intf_nwxc_ft97_expei.fh" 523c 524 type(nwad_dble)::RHOA,RHOB,RHOA13,RHOB13,AMA,AMB,EAA,EBB,EAB,EBA 525 type(nwad_dble)::EEI,MUAB,MUBA,MUBB,MUAA,GA2,GB2,CFT,RSA,RSB 526 type(nwad_dble)::R,GR,T,EEI1,DENOM,FA 527c type(nwad_dble)::ECMUC0 528c type(nwad_dble)::KSSP0,KSS0,FSSP,FSS,FACTOR,FF 529c type(nwad_dble)::USSP1,USSP2,WSSP,USS1,USS2,WSS,YSS,DFFDMU 530 double precision KSSPBIG,big,kss0big 531c 532c Nwxc_FT97_EXPEI is already declared in the corresponding 533c interface block. 534c 535#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 536c type(nwad_dble)::nwxc_FT97_EXPEI 537#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 538c type(nwad_dble)::nwxc_FT97_EXPEI_d2 539#else 540c type(nwad_dble)::nwxc_FT97_EXPEI_d3 541#endif 542C NOTATION FOR CONSTANTS. 543C C0 = (1-ln2)/(2*Pi**2), see eq.(9) 544C C1 = 2*C0/3 , see eq.(13,(28),(33) 545C C2 = (3/(4*Pi))**1/3 , see eq.(8) 546C C3 = C2/3 547 PARAMETER (C0=0.1554534543483D-1) 548 PARAMETER (C1=0.2072712724644D-1) 549 PARAMETER (C2=0.6203504908994D00) 550 PARAMETER (C3=0.2067834969665D00) 551C OPTIMIZED PARAMETERS IN GRADIENT-CORRECTED CORRELATION xc_ft97. 552C See formulas for FSSP, eq.(45), and FSS, eq.(46). 553 PARAMETER (A1=1.622118767D0) 554 PARAMETER (A2=0.489958076D0) 555 PARAMETER (A3=1.379021941D0) 556 PARAMETER (A4=4.946281353D0) 557 PARAMETER (A5=3.600612059D0) 558C NUMERICAL CUTOFF FOR MUAA, MUAB, MUBA, MUBB, see eqs.(13) and (33). 559 PARAMETER (CUTOFF=1.0D7) 560 parameter(ksspbig=1.291551074D0 - 0.349064173D0,big=1d4) 561 parameter(KSS0BIG = 1.200801774D0 + 0.859614445D0- 562 - 0.812904345D0) 563C 564C *** DEFINITION OF FORMULA FUNCTIONS. 565C 566C R : Density parameter Rs (Wigner radius), see eq.(8). 567C KSSP0 : See eq.(39) for k(sigma,sigma'). 568chvd KSSP0(R) = 1.291551074D0 - 0.349064173D0*(1.D0- 569chvd . EXP(-0.083275880D0*R**0.8D0)) 570C KSS0 : See eq.(40) for k(sigma,sigma). 571chvd KSS0(R) = 1.200801774D0 + 0.859614445D0*(1.D0- 572chvd . EXP(-1.089338848d0*SQRT(R)))- 573chvd . 0.812904345D0*(1.D0-EXP(-0.655638823D0*R**0.4D0)) 574C 575C FSSP : See eq.(45) for F(sigma,sigma'). 576chvd FSSP(R,GR) = (1.D0+A1*GR+(A2*GR)**2.0d0)* 577chvd . EXP(-(A2*GR)**2.0d0)/SQRT(1.D0+A3*GR/R) 578C USSP1 : Derivative of ln(k(sigma,sigma')) with respect to Rs. 579c USSP1(R) = -4.65098019D-2*EXP(-0.083275880D0*R**0.8D0)*R**0.8D0 580C USSP2 : Derivative of ln(F(sigma,sigma')) with respect to Rs. 581c USSP2(R,GR) = A3*GR/(R+A3*GR) 582C WSSP : Derivative of ln(F(sigma,sigma')) w.r.t. Nabla(Rs)^2 583c WSSP(R,GR) = ((A1+A1)-A1*(A2+A2)**2.0d0*GR**2.0d0 584c . -(A2*(A2+A2))**2.0d0*GR**3.0d0) 585c . /(1.D0+A1*GR+(A2*GR)**2.0d0) - A3/(R+A3*GR) 586C 587C FSS : See eq.(46) for F(sigma,sigma). 588chvd FSS(R,GR)= (1.D0+(A4*GR)**2.0d0)*EXP(-(A4*GR)**2.0d0)/ 589chvd . SQRT(1.D0+A5*GR/R) 590C FACTOR : See eq.(34). 591chvd FACTOR(R)= EXP(-(R/(0.939016D0*SQRT(R)+1.733170D0*R))**2.0d0) 592C USS1 : Derivative of ln(k(sigma,sigma)) with respect to Rs. 593c USS1(R) =-4.26377318D-1*R**0.4D0*EXP(-0.655638823D0*R**0.4D0)+ 594c . 0.93641140924D0*SQRT(R)*EXP(-1.089338848d0*SQRT(R)) 595C USS2 : Derivative of ln(F(sigma,sigma)) with respect to Rs. 596c USS2(R,GR)= A5*GR/(R+A5*GR) 597C WSS : Derivative of ln(F(sigma,sigma)) w.r.t. Nabla(Rs)^2 598c WSS(R,GR) =-(A4*(A4+A4))**2.0d0*GR**3.0d0/(1.D0+(A4*GR)**2.0d0) 599c . -A5/(R+A5*GR) 600C YSS : Derivative of FACTOR (eq.(34)) with respect to Rs. 601c YSS(R) = 0.939016D0*SQRT(R)*R*R/(0.939016D0*SQRT(R)+ 602c . 1.733170D0*R)**3.0d0 603C 604C FF : Used in approximate normalization, eq.(15). 605chvd FF(T) = (3.D0+2.D0*(SQRT(T)+T))/ 606chvd . (3.D0+6.D0*(SQRT(T)+T)) 607C DFFDMU : Derivative of FF with respect to the argument T. 608c DFFDMU(T) =-(2.D0*(SQRT(T)+T+T))/ 609c . (3.D0*(1.D0+2.D0*(SQRT(T)+T))**2.0d0) 610C 611C *** RSA - Rs(alpha), GA - grad(Rs(alpha)) 612C 613 IF(RHOA.GT.tol_rho) THEN 614C Wigner radius Rs(alpha), eq.(8) 615 RSA = C2/RHOA13 616 CFT = C3/(RHOA*RHOA13) 617C GA : Nabla(Rs(alpha)) 618C GA2 : Square of Nabla(Rs(alpha)) 619 GA2 = CFT*CFT*AMA 620C mu(beta,alpha), eq.(13) 621C MUBA = C1*RSA/(KSSP0(RSA)*FSSP(RSA,GA2))**2 622 EBA = 0.D0 623c ECMUC0 = 0.D0 624 DEBADR = 0.D0 625 DEBADG = 0.D0 626 if((ga2*a2)**2.0d0.gt.big) then 627 DENOM = 0d0 628 else 629 if(rsa.gt.big) then 630 DENOM = (KSSPBIG*FSSP(RSA,GA2))**2.0d0 631 else 632 DENOM = (KSSP0(RSA)*FSSP(RSA,GA2))**2.0d0 633 endif 634 endif 635 IF(C1*RSA .LE. DENOM*CUTOFF) THEN 636 MUBA = C1*RSA/DENOM 637#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 638#if defined(NWAD_PRINT) 639 EEI = nwxc_FT97_EXPEI_p(-MUBA) 640#else 641 EEI = nwxc_FT97_EXPEI(-MUBA) 642#endif 643#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 644 EEI = nwxc_FT97_EXPEI_d2(-MUBA) 645#else 646 EEI = nwxc_FT97_EXPEI_d3(-MUBA) 647#endif 648 EEI1 = MUBA*EEI+1.D0 649C EBA : Correlation energy density, eq.(19) 650 EBA = C0*(EEI+2.D0*FF(MUBA)*EEI1) 651C ECMUC0 : Derivative of EBA w.r.t. MUBA 652c ECMUC0 = C0*(EEI1*(1.D0+2.D0*(DFFDMU(MUBA)+MUBA* 653c . FF(MUBA)))+2.D0*FF(MUBA)*MUBA*EEI) 654C d(Ec)/d(Nabla(rho)^2) = -CFT^2*(mu*dEcdmu)*Wss' 655c DEBADG = -WSSP(RSA,GA2)*ECMUC0*CFT*CFT 656C d(Ec)/d(rho) 657c DEBADR = -ECMUC0*(1.D0-USSP1(RSA)/KSSP0(RSA)- 658c . USSP2(RSA,GA2) - 659c . 8.D0*WSSP(RSA,GA2)*GA2) 660c . / (3.D0*RHOA) 661 ENDIF 662C mu(alpha,alpha), eq.(33) 663C MUAA = C1*RSA/(KSS0(RSA)*FSS(RSA,GA2))**2 664 FA = FACTOR(RSA) 665 EAA = 0.D0 666c ECMUC0 = 0.D0 667c DEAADR = 0.D0 668c DEAADG = 0.D0 669 if((ga2*a2)**2.0d0.gt.big) then 670 DENOM = 0d0 671 else 672 if(rsa.gt.big) then 673 DENOM = (KSS0big*FSS(RSA,GA2))**2.0d0 674 else 675 DENOM = (KSS0(RSA)*FSS(RSA,GA2))**2.0d0 676 endif 677 endif 678 IF(C1*RSA .LE. DENOM*CUTOFF) THEN 679 MUAA = C1*RSA/DENOM 680#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 681#if defined(NWAD_PRINT) 682 EEI = nwxc_FT97_EXPEI_p(-MUAA) 683#else 684 EEI = nwxc_FT97_EXPEI(-MUAA) 685#endif 686#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 687 EEI = nwxc_FT97_EXPEI_d2(-MUAA) 688#else 689 EEI = nwxc_FT97_EXPEI_d3(-MUAA) 690#endif 691 EEI1 = MUAA*EEI+1.D0 692C EAA : Correlation energy density, eqs.(19) and (31). 693 EAA = FA*C0*(EEI+2.D0*FF(MUAA)*EEI1) 694C ECMUC0 : Derivative of EAA w.r.t. MUAA 695c ECMUC0 = FA*C0*(EEI1*(1.D0+2.D0*(DFFDMU(MUAA)+MUAA* 696c . FF(MUAA)))+2.D0*FF(MUAA)*MUAA*EEI) 697C d(Ec)/d(rho) 698c DEAADR = (-ECMUC0*(1.D0-USS1(RSA)/KSS0(RSA)-USS2(RSA,GA2) - 699c . 8.D0*WSS(RSA,GA2)*GA2) + YSS(RSA)*EAA) 700c . / (3.D0*RHOA) 701C d(Ec)/d(Nabla(rho)^2) = -Cft^2*(mu*dEcdmu*Fa)*Wss 702c DEAADG = -WSS(RSA,GA2)*ECMUC0*CFT*CFT 703 ENDIF 704 ELSE 705 EAA = 0.D0 706 EBA = 0.D0 707c DEAADR = 0.D0 708c DEBADR = 0.D0 709c DEAADG = 0.D0 710c DEBADG = 0.D0 711 ENDIF 712C *** RSB - Rs(beta), GB - grad(Rs(beta)) 713 IF(RHOB.GT.tol_rho) THEN 714C Wigner radius Rs(beta), eq.(8) 715 RSB = C2/RHOB13 716 CFT = C3/(RHOB*RHOB13) 717C GB : Nabla(Rs(beta)) 718C GB2 : Square of Nabla(Rs(beta)) 719 GB2 = CFT*CFT*AMB 720C mu(alpha,beta), eq.(13) 721C MUAB = C1*RSB/(KSSP0(RSB)*FSSP(RSB,GB2))**2 722 EAB = 0.D0 723c ECMUC0 = 0.D0 724 DEABDR = 0.D0 725 DEABDG = 0.D0 726 DENOM = (KSSP0(RSB)*FSSP(RSB,GB2))**2.0d0 727 if((gb2*a2)**2.0d0.gt.big) then 728 DENOM = 0d0 729 else 730 if(rsb.gt.big) then 731 DENOM = (KSSPBIG*FSSP(RSB,GB2))**2.0d0 732 else 733 DENOM = (KSSP0(RSB)*FSSP(RSB,GB2))**2.0d0 734 endif 735 endif 736 IF(C1*RSB .LE. DENOM*CUTOFF) THEN 737 MUAB = C1*RSB/DENOM 738#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 739#if defined(NWAD_PRINT) 740 EEI = nwxc_FT97_EXPEI_p(-MUAB) 741#else 742 EEI = nwxc_FT97_EXPEI(-MUAB) 743#endif 744#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 745 EEI = nwxc_FT97_EXPEI_d2(-MUAB) 746#else 747 EEI = nwxc_FT97_EXPEI_d3(-MUAB) 748#endif 749 EEI1 = MUAB*EEI+1.D0 750C EAB : Correlation energy density, eqs.(19). 751 EAB = C0*(EEI+2.D0*FF(MUAB)*EEI1) 752C ECMUC0 : Derivative of EAB w.r.t. MUAB 753c ECMUC0 = C0*(EEI1*(1.D0+2.D0*(DFFDMU(MUAB)+MUAB* 754c . FF(MUAB)))+2.D0*FF(MUAB)*MUAB*EEI) 755C d(Ec)/d(Nabla(rho)^2) = -CFT^2*(mu*dEcdmu)*Wss' 756c DEABDG = -WSSP(RSB,GB2)*ECMUC0*CFT*CFT 757C d(Ec)/d(rho) 758c DEABDR = -ECMUC0*(1.D0-USSP1(RSB)/KSSP0(RSB)- 759c . USSP2(RSB,GB2) - 760c . 8.D0*WSSP(RSB,GB2)*GB2) 761c . / (3.D0*RHOB) 762 ENDIF 763C mu(beta,beta), eq.(33) 764C MUBB = C1*RSB/(KSS0(RSB)*FSS(RSB,GB2))**2 765 FA = FACTOR(RSB) 766 EBB = 0.D0 767c ECMUC0 = 0.D0 768c DEBBDR = 0.D0 769c DEBBDG = 0.D0 770 if((gb2*a2)**2.0d0.gt.big) then 771 DENOM = 0d0 772 else 773 if(rsb.gt.big) then 774 DENOM = (KSS0big*FSS(RSB,GB2))**2.0d0 775 else 776 DENOM = (KSS0(RSB)*FSS(RSB,GB2))**2.0d0 777 endif 778 endif 779 IF(C1*RSB .LE. DENOM*CUTOFF) THEN 780 MUBB = C1*RSB/DENOM 781#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 782#if defined(NWAD_PRINT) 783 EEI = nwxc_FT97_EXPEI_p(-MUBB) 784#else 785 EEI = nwxc_FT97_EXPEI(-MUBB) 786#endif 787#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 788 EEI = nwxc_FT97_EXPEI_d2(-MUBB) 789#else 790 EEI = nwxc_FT97_EXPEI_d3(-MUBB) 791#endif 792 EEI1 = MUBB*EEI+1.D0 793C EBB : Correlation energy density, eqs.(19) and (31). 794 EBB = FA*C0*(EEI+2.D0*FF(MUBB)*EEI1) 795C ECMUC0 : Derivative of EBB w.r.t. MUBB 796c ECMUC0 = FA*C0*(EEI1*(1.D0+2.D0*(DFFDMU(MUBB)+MUBB* 797c . FF(MUBB)))+2.D0*FF(MUBB)*MUBB*EEI) 798C d(Ec)/d(rho) 799c DEBBDR = (-ECMUC0*(1.D0-USS1(RSB)/KSS0(RSB) 800c . -USS2(RSB,GB2) - 801c . 8.D0*WSS(RSB,GB2)*GB2) + YSS(RSB)*EBB) 802c . / (3.D0*RHOB) 803C d(Ec)/d(Nabla(rho)^2) = -Cft^2*(mu*dEcdmu*Fa)*Wss 804c DEBBDG = -WSS(RSB,GB2)*ECMUC0*CFT*CFT 805 ENDIF 806 ELSE 807 EAB = 0.D0 808 EBB = 0.D0 809c DEABDR = 0.D0 810c DEBBDR = 0.D0 811c DEABDG = 0.D0 812c DEBBDG = 0.D0 813 ENDIF 814 RETURN 815c 816 CONTAINS 817c 818c The combination of statement functions and derived types with 819c overloaded operators is not properly supported in the Fortran 820c standard (apparently). Therefore the statement functions from 821c the original subroutine had to be transformed into contained 822c functions. 823c 824c WARNING: This code is EXTREMELY EVIL! Although there appears to be 825c only one function there are actually three with the same name, 826c each one returning results of a different data type. The trick is 827c that depending on the data type the the subroutine that contains 828c these functions changes its name and hence these different 829c functions of the same name do not lead to conflicts. The 830c alternative would have been to add a forest of conditional 831c compilation constructs (#ifdef's) to change the function names 832c in the declarations and all places where they are used. That 833c would have been extremely ugly, so we are between a rock and a 834c hard place on this one. 835c 836 function KSS0(R) result(s) 837#include "nwad.fh" 838 implicit none 839 type(nwad_dble), intent(in) :: R 840 type(nwad_dble) :: S 841 S = 1.200801774D0 + 0.859614445D0*(1.D0- 842 . EXP(-1.089338848d0*SQRT(R)))- 843 . 0.812904345D0*(1.D0-EXP(-0.655638823D0*R**0.4D0)) 844 end function 845c 846 function KSSP0(R) result(s) 847#include "nwad.fh" 848 implicit none 849 type(nwad_dble), intent(in) :: R 850 type(nwad_dble) :: S 851 S = 1.291551074D0 - 0.349064173D0*(1.D0- 852 . EXP(-0.083275880D0*R**0.8D0)) 853 end function 854c 855 function FSSP(R,GR) result(s) 856#include "nwad.fh" 857 implicit none 858 type(nwad_dble), intent(in) :: R 859 type(nwad_dble), intent(in) :: GR 860 type(nwad_dble) :: S 861 S = (1.D0+A1*GR+(A2*GR)**2.0d0)* 862 . EXP(-(A2*GR)**2.0d0)/SQRT(1.D0+A3*GR/R) 863 end function 864c 865 function FSS(R,GR) result(s) 866#include "nwad.fh" 867 implicit none 868 type(nwad_dble), intent(in) :: R 869 type(nwad_dble), intent(in) :: GR 870 type(nwad_dble) :: S 871 S = (1.D0+(A4*GR)**2.0d0)*EXP(-(A4*GR)**2.0d0)/ 872 . SQRT(1.D0+A5*GR/R) 873 end function 874c 875 function FACTOR(R) result(s) 876#include "nwad.fh" 877 implicit none 878 type(nwad_dble), intent(in) :: R 879 type(nwad_dble) :: S 880 S = EXP(-(R/(0.939016D0*SQRT(R)+1.733170D0*R))**2.0d0) 881 end function 882c 883 function FF(T) result(s) 884#include "nwad.fh" 885 implicit none 886 type(nwad_dble), intent(in) :: T 887 type(nwad_dble) :: S 888 S = (3.D0+2.D0*(SQRT(T)+T))/ 889 . (3.D0+6.D0*(SQRT(T)+T)) 890 end function 891c 892 END 893c----------------------------------------------------------------------- 894#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 895#if defined(NWAD_PRINT) 896 FUNCTION nwxc_FT97_EXPEI_p(X) 897#else 898 FUNCTION nwxc_FT97_EXPEI(X) 899#endif 900#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 901 FUNCTION nwxc_FT97_EXPEI_d2(X) 902#else 903 FUNCTION nwxc_FT97_EXPEI_d3(X) 904#endif 905C 906C This function program computes approximate values for the 907C function exp(-x) * Ei(x), where Ei(x) is the exponential 908C integral, and x is real. 909C 910C Author: W. J. Cody 911C 912C Latest modification: January 12, 1988 913C 914#include "nwad.fh" 915c 916#include "intf_nwxc_ft97_calcei.fh" 917c 918 INTEGER INT 919 type(nwad_dble),intent(in):: X 920 type(nwad_dble):: nwxc_FT97_EXPEI, RESULT 921 type(nwad_dble):: nwxc_FT97_EXPEI_p 922 type(nwad_dble):: nwxc_FT97_EXPEI_d2 923 type(nwad_dble):: nwxc_FT97_EXPEI_d3 924 INT = 3 925#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 926#if defined(NWAD_PRINT) 927 CALL nwxc_FT97_CALCEI_p(X,RESULT,INT) 928 nwxc_FT97_EXPEI_p = RESULT 929#else 930 CALL nwxc_FT97_CALCEI(X,RESULT,INT) 931 nwxc_FT97_EXPEI = RESULT 932#endif 933#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 934 CALL nwxc_FT97_CALCEI_d2(X,RESULT,INT) 935 nwxc_FT97_EXPEI_d2 = RESULT 936#else 937 CALL nwxc_FT97_CALCEI_d3(X,RESULT,INT) 938 nwxc_FT97_EXPEI_d3 = RESULT 939#endif 940 RETURN 941 END 942c----------------------------------------------------------------------- 943#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 944#if defined(NWAD_PRINT) 945 SUBROUTINE nwxc_FT97_CALCEI_p(ARG,RESULT,INT) 946#else 947 SUBROUTINE nwxc_FT97_CALCEI(ARG,RESULT,INT) 948#endif 949#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 950 SUBROUTINE nwxc_FT97_CALCEI_d2(ARG,RESULT,INT) 951#else 952 SUBROUTINE nwxc_FT97_CALCEI_d3(ARG,RESULT,INT) 953#endif 954C 955C This Fortran 77 packet computes the exponential integrals Ei(x), 956C E1(x), and exp(-x)*Ei(x) for real arguments x where 957C 958C integral (from t=-infinity to t=x) (exp(t)/t), x > 0, 959C Ei(x) = 960C -integral (from t=-x to t=infinity) (exp(t)/t), x < 0, 961C 962C and where the first integral is a principal value integral. 963C The packet contains three function type subprograms: EI, EONE, 964C and EXPEI; and one subroutine type subprogram: CALCEI. The 965C calling statements for the primary entries are 966C 967C Y = EI(X), where X .NE. 0, 968C 969C Y = EONE(X), where X .GT. 0, 970C and 971C Y = EXPEI(X), where X .NE. 0, 972C 973C and where the entry points correspond to the functions Ei(x), 974C E1(x), and exp(-x)*Ei(x), respectively. The routine CALCEI 975C is intended for internal packet use only, all computations within 976C the packet being concentrated in this routine. The function 977C subprograms invoke CALCEI with the Fortran statement 978C CALL CALCEI(ARG,RESULT,INT) 979C where the parameter usage is as follows 980C 981C Function Parameters for CALCEI 982C Call ARG RESULT INT 983C 984C EI(X) X .NE. 0 Ei(X) 1 985C EONE(X) X .GT. 0 -Ei(-X) 2 986C EXPEI(X) X .NE. 0 exp(-X)*Ei(X) 3 987C 988C The main computation involves evaluation of rational Chebyshev 989C approximations published in Math. Comp. 22, 641-649 (1968), and 990C Math. Comp. 23, 289-303 (1969) by Cody and Thacher. This 991C transportable program is patterned after the machine-dependent 992C FUNPACK packet NATSEI, but cannot match that version for 993C efficiency or accuracy. This version uses rational functions 994C that theoretically approximate the exponential integrals to 995C at least 18 significant decimal digits. The accuracy achieved 996C depends on the arithmetic system, the compiler, the intrinsic 997C functions, and proper selection of the machine-dependent 998C constants. 999C 1000C 1001C******************************************************************* 1002C******************************************************************* 1003C 1004C Explanation of machine-dependent constants 1005C 1006C beta = radix for the floating-point system. 1007C minexp = smallest representable power of beta. 1008C maxexp = smallest power of beta that overflows. 1009C XBIG = largest argument acceptable to EONE; solution to 1010C equation: 1011C exp(-x)/x * (1 + 1/x) = beta ** minexp. 1012C XINF = largest positive machine number; approximately 1013C beta ** maxexp 1014C XMAX = largest argument acceptable to EI; solution to 1015C equation: exp(x)/x * (1 + 1/x) = beta ** maxexp. 1016C 1017C Approximate values for some important machines are: 1018C 1019C beta minexp maxexp 1020C 1021C CRAY-1 (S.P.) 2 -8193 8191 1022C Cyber 180/185 1023C under NOS (S.P.) 2 -975 1070 1024C IEEE (IBM/XT, 1025C SUN, etc.) (S.P.) 2 -126 128 1026C IEEE (IBM/XT, 1027C SUN, etc.) (D.P.) 2 -1022 1024 1028C IBM 3033 (D.P.) 16 -65 63 1029C VAX D-Format (D.P.) 2 -128 127 1030C VAX G-Format (D.P.) 2 -1024 1023 1031C 1032C XBIG XINF XMAX 1033C 1034C CRAY-1 (S.P.) 5670.31 5.45E+2465 5686.21 1035C Cyber 180/185 1036C under NOS (S.P.) 669.31 1.26E+322 748.28 1037C IEEE (IBM/XT, 1038C SUN, etc.) (S.P.) 82.93 3.40E+38 93.24 1039C IEEE (IBM/XT, 1040C SUN, etc.) (D.P.) 701.84 1.79D+308 716.35 1041C IBM 3033 (D.P.) 175.05 7.23D+75 179.85 1042C VAX D-Format (D.P.) 84.30 1.70D+38 92.54 1043C VAX G-Format (D.P.) 703.22 8.98D+307 715.66 1044C 1045C******************************************************************* 1046C******************************************************************* 1047C 1048C Error returns 1049C 1050C The following table shows the types of error that may be 1051C encountered in this routine and the function value supplied 1052C in each case. 1053C 1054C Error Argument Function values for 1055C Range EI EXPEI EONE 1056C 1057C UNDERFLOW (-)X .GT. XBIG 0 - 0 1058C OVERFLOW X .GE. XMAX XINF - - 1059C ILLEGAL X X = 0 -XINF -XINF XINF 1060C ILLEGAL X X .LT. 0 - - USE ABS(X) 1061C 1062C Intrinsic functions required are: 1063C 1064C ABS, SQRT, EXP 1065C 1066C 1067C Author: W. J. Cody 1068C Mathematics abd Computer Science Division 1069C Argonne National Laboratory 1070C Argonne, IL 60439 1071C 1072C Latest modification: September 9, 1988 1073C 1074C---------------------------------------------------------------------- 1075c 1076c The approach followed here to approximate the derivatives of the 1077c integral is to differentiate the Chebychev approximation to the 1078c the integral. Strictly speaking this is not correct but given the 1079c implementation of EXPEI we have this the only feasible approach. 1080c 1081#include "nwad.fh" 1082 implicit none 1083 INTEGER,intent(in):: INT 1084 INTEGER I 1085 type(nwad_dble),intent(in)::arg 1086 type(nwad_dble),intent(out)::result 1087 type(nwad_dble):: 1088 1 EI,FRAC, 1089 2 PX,QX,SUMP, 1090 3 SUMQ,T,W,X,XMX0, 1091 4 Y,YSQ 1092 double precision 1093 1 A,B,C,D,EXP40,E,F,FOUR,FOURTY,HALF,ONE,P, 1094 2 PLG,P037,P1,P2,q,QLG,Q1,Q2,R,S,SIX, 1095 3 THREE,TWELVE,TWO,TWO4,XBIG,XINF,XMAX, 1096 4 X0,X01,X02,X11,ZERO 1097 DIMENSION A(7),B(6),C(9),D(9),E(10),F(10),P(10),q(10),R(10), 1098 1 S(9),P1(10),Q1(9),P2(10),Q2(9),PLG(4),QLG(4),PX(10),QX(10) 1099C---------------------------------------------------------------------- 1100C Mathematical constants 1101C EXP40 = exp(40) 1102C X0 = zero of Ei 1103C X01/X11 + X02 = zero of Ei to extra precision 1104C---------------------------------------------------------------------- 1105 DATA ZERO,P037,HALF,ONE,TWO/0.0D0,0.037D0,0.5D0,1.0D0,2.0D0/, 1106 1 THREE,FOUR,SIX,TWELVE,TWO4/3.0D0,4.0D0,6.0D0,12.D0,24.0D0/, 1107 2 FOURTY,EXP40/40.0D0,2.3538526683701998541D17/, 1108 3 X01,X11,X02/381.5D0,1024.0D0,-5.1182968633365538008D-5/, 1109 4 X0/3.7250741078136663466D-1/ 1110C---------------------------------------------------------------------- 1111C Machine-dependent constants 1112C---------------------------------------------------------------------- 1113CS DATA XINF/3.40E+38/,XMAX/93.246E0/,XBIG/82.93E0/ 1114 DATA XINF/1.79D+308/,XMAX/716.351D0/,XBIG/701.84D0/ 1115C---------------------------------------------------------------------- 1116C Coefficients for -1.0 <= X < 0.0 1117C---------------------------------------------------------------------- 1118 DATA A/1.1669552669734461083368D2, 2.1500672908092918123209D3, 1119 1 1.5924175980637303639884D4, 8.9904972007457256553251D4, 1120 2 1.5026059476436982420737D5,-1.4815102102575750838086D5, 1121 3 5.0196785185439843791020D0/ 1122 DATA B/4.0205465640027706061433D1, 7.5043163907103936624165D2, 1123 1 8.1258035174768735759855D3, 5.2440529172056355429883D4, 1124 2 1.8434070063353677359298D5, 2.5666493484897117319268D5/ 1125C---------------------------------------------------------------------- 1126C Coefficients for -4.0 <= X < -1.0 1127C---------------------------------------------------------------------- 1128 DATA C/3.828573121022477169108D-1, 1.107326627786831743809D+1, 1129 1 7.246689782858597021199D+1, 1.700632978311516129328D+2, 1130 2 1.698106763764238382705D+2, 7.633628843705946890896D+1, 1131 3 1.487967702840464066613D+1, 9.999989642347613068437D-1, 1132 4 1.737331760720576030932D-8/ 1133 DATA D/8.258160008564488034698D-2, 4.344836335509282083360D+0, 1134 1 4.662179610356861756812D+1, 1.775728186717289799677D+2, 1135 2 2.953136335677908517423D+2, 2.342573504717625153053D+2, 1136 3 9.021658450529372642314D+1, 1.587964570758947927903D+1, 1137 4 1.000000000000000000000D+0/ 1138C---------------------------------------------------------------------- 1139C Coefficients for X < -4.0 1140C---------------------------------------------------------------------- 1141 DATA E/1.3276881505637444622987D+2,3.5846198743996904308695D+4, 1142 1 1.7283375773777593926828D+5,2.6181454937205639647381D+5, 1143 2 1.7503273087497081314708D+5,5.9346841538837119172356D+4, 1144 3 1.0816852399095915622498D+4,1.0611777263550331766871D03, 1145 4 5.2199632588522572481039D+1,9.9999999999999999087819D-1/ 1146 DATA F/3.9147856245556345627078D+4,2.5989762083608489777411D+5, 1147 1 5.5903756210022864003380D+5,5.4616842050691155735758D+5, 1148 2 2.7858134710520842139357D+5,7.9231787945279043698718D+4, 1149 3 1.2842808586627297365998D+4,1.1635769915320848035459D+3, 1150 4 5.4199632588522559414924D+1,1.0D0/ 1151C---------------------------------------------------------------------- 1152C Coefficients for rational approximation to ln(x/a), |1-x/a| < .1 1153C---------------------------------------------------------------------- 1154 DATA PLG/-2.4562334077563243311D+01,2.3642701335621505212D+02, 1155 1 -5.4989956895857911039D+02,3.5687548468071500413D+02/ 1156 DATA QLG/-3.5553900764052419184D+01,1.9400230218539473193D+02, 1157 1 -3.3442903192607538956D+02,1.7843774234035750207D+02/ 1158C---------------------------------------------------------------------- 1159C Coefficients for 0.0 < X < 6.0, 1160C ratio of Chebyshev polynomials 1161C---------------------------------------------------------------------- 1162 DATA P/-1.2963702602474830028590D01,-1.2831220659262000678155D03, 1163 1 -1.4287072500197005777376D04,-1.4299841572091610380064D06, 1164 2 -3.1398660864247265862050D05,-3.5377809694431133484800D08, 1165 3 3.1984354235237738511048D08,-2.5301823984599019348858D10, 1166 4 1.2177698136199594677580D10,-2.0829040666802497120940D11/ 1167 DATA q/ 7.6886718750000000000000D01,-5.5648470543369082846819D03, 1168 1 1.9418469440759880361415D05,-4.2648434812177161405483D06, 1169 2 6.4698830956576428587653D07,-7.0108568774215954065376D08, 1170 3 5.4229617984472955011862D09,-2.8986272696554495342658D10, 1171 4 9.8900934262481749439886D10,-8.9673749185755048616855D10/ 1172C---------------------------------------------------------------------- 1173C J-fraction coefficients for 6.0 <= X < 12.0 1174C---------------------------------------------------------------------- 1175 DATA R/-2.645677793077147237806D00,-2.378372882815725244124D00, 1176 1 -2.421106956980653511550D01, 1.052976392459015155422D01, 1177 2 1.945603779539281810439D01,-3.015761863840593359165D01, 1178 3 1.120011024227297451523D01,-3.988850730390541057912D00, 1179 4 9.565134591978630774217D00, 9.981193787537396413219D-1/ 1180 DATA S/ 1.598517957704779356479D-4, 4.644185932583286942650D00, 1181 1 3.697412299772985940785D02,-8.791401054875438925029D00, 1182 2 7.608194509086645763123D02, 2.852397548119248700147D01, 1183 3 4.731097187816050252967D02,-2.369210235636181001661D02, 1184 4 1.249884822712447891440D00/ 1185C---------------------------------------------------------------------- 1186C J-fraction coefficients for 12.0 <= X < 24.0 1187C---------------------------------------------------------------------- 1188 DATA P1/-1.647721172463463140042D00,-1.860092121726437582253D01, 1189 1 -1.000641913989284829961D01,-2.105740799548040450394D01, 1190 2 -9.134835699998742552432D-1,-3.323612579343962284333D01, 1191 3 2.495487730402059440626D01, 2.652575818452799819855D01, 1192 4 -1.845086232391278674524D00, 9.999933106160568739091D-1/ 1193 DATA Q1/ 9.792403599217290296840D01, 6.403800405352415551324D01, 1194 1 5.994932325667407355255D01, 2.538819315630708031713D02, 1195 2 4.429413178337928401161D01, 1.192832423968601006985D03, 1196 3 1.991004470817742470726D02,-1.093556195391091143924D01, 1197 4 1.001533852045342697818D00/ 1198C---------------------------------------------------------------------- 1199C J-fraction coefficients for X .GE. 24.0 1200C---------------------------------------------------------------------- 1201 DATA P2/ 1.75338801265465972390D02,-2.23127670777632409550D02, 1202 1 -1.81949664929868906455D01,-2.79798528624305389340D01, 1203 2 -7.63147701620253630855D00,-1.52856623636929636839D01, 1204 3 -7.06810977895029358836D00,-5.00006640413131002475D00, 1205 4 -3.00000000320981265753D00, 1.00000000000000485503D00/ 1206 DATA Q2/ 3.97845977167414720840D04, 3.97277109100414518365D00, 1207 1 1.37790390235747998793D02, 1.17179220502086455287D02, 1208 2 7.04831847180424675988D01,-1.20187763547154743238D01, 1209 3 -7.99243595776339741065D00,-2.99999894040324959612D00, 1210 4 1.99999999999048104167D00/ 1211C---------------------------------------------------------------------- 1212 X = ARG 1213 IF (X .EQ. ZERO) THEN 1214 EI = -XINF 1215 IF (INT .EQ. 2) EI = -EI 1216 ELSE IF ((X .LT. ZERO) .OR. (INT .EQ. 2)) THEN 1217C---------------------------------------------------------------------- 1218C Calculate EI for negative argument or for E1. 1219C---------------------------------------------------------------------- 1220 Y = ABS(X) 1221 IF (Y .LE. ONE) THEN 1222 SUMP = A(7) * Y + A(1) 1223 SUMQ = Y + B(1) 1224 DO 110 I = 2, 6 1225 SUMP = SUMP * Y + A(I) 1226 SUMQ = SUMQ * Y + B(I) 1227 110 CONTINUE 1228 EI = LOG(Y) - SUMP / SUMQ 1229 IF (INT .EQ. 3) EI = EI * EXP(Y) 1230 ELSE IF (Y .LE. FOUR) THEN 1231 W = ONE / Y 1232 SUMP = C(1) 1233 SUMQ = D(1) 1234 DO 130 I = 2, 9 1235 SUMP = SUMP * W + C(I) 1236 SUMQ = SUMQ * W + D(I) 1237 130 CONTINUE 1238 EI = - SUMP / SUMQ 1239 IF (INT .NE. 3) EI = EI * EXP(-Y) 1240 ELSE 1241 IF ((Y .GT. XBIG) .AND. (INT .LT. 3)) THEN 1242 EI = ZERO 1243 ELSE 1244 W = ONE / Y 1245 SUMP = E(1) 1246 SUMQ = F(1) 1247 DO 150 I = 2, 10 1248 SUMP = SUMP * W + E(I) 1249 SUMQ = SUMQ * W + F(I) 1250 150 CONTINUE 1251 EI = -W * (ONE - W * SUMP / SUMQ ) 1252 IF (INT .NE. 3) EI = EI * EXP(-Y) 1253 END IF 1254 END IF 1255 IF (INT .EQ. 2) EI = -EI 1256 ELSE IF (X .LT. SIX) THEN 1257C---------------------------------------------------------------------- 1258C To improve conditioning, rational approximations are expressed 1259C in terms of Chebyshev polynomials for 0 <= X < 6, and in 1260C continued fraction form for larger X. 1261C---------------------------------------------------------------------- 1262 T = X + X 1263 T = T / THREE - TWO 1264 PX(1) = ZERO 1265 QX(1) = ZERO 1266 PX(2) = P(1) 1267 QX(2) = q(1) 1268 DO 210 I = 2, 9 1269 PX(I+1) = T * PX(I) - PX(I-1) + P(I) 1270 QX(I+1) = T * QX(I) - QX(I-1) + q(I) 1271 210 CONTINUE 1272 SUMP = HALF * T * PX(10) - PX(9) + P(10) 1273 SUMQ = HALF * T * QX(10) - QX(9) + q(10) 1274 FRAC = SUMP / SUMQ 1275 XMX0 = (X - X01/X11) - X02 1276 IF (ABS(XMX0) .GE. P037) THEN 1277 EI = LOG(X/X0) + XMX0 * FRAC 1278 IF (INT .EQ. 3) EI = EXP(-X) * EI 1279 ELSE 1280C---------------------------------------------------------------------- 1281C Special approximation to ln(X/X0) for X close to X0 1282C---------------------------------------------------------------------- 1283 Y = XMX0 / (X + X0) 1284 YSQ = Y*Y 1285 SUMP = PLG(1) 1286 SUMQ = YSQ + QLG(1) 1287 DO 220 I = 2, 4 1288 SUMP = SUMP*YSQ + PLG(I) 1289 SUMQ = SUMQ*YSQ + QLG(I) 1290 220 CONTINUE 1291 EI = (SUMP / (SUMQ*(X+X0)) + FRAC) * XMX0 1292 IF (INT .EQ. 3) EI = EXP(-X) * EI 1293 END IF 1294 ELSE IF (X .LT. TWELVE) THEN 1295 FRAC = ZERO 1296 DO 230 I = 1, 9 1297 FRAC = S(I) / (R(I) + X + FRAC) 1298 230 CONTINUE 1299 EI = (R(10) + FRAC) / X 1300 IF (INT .NE. 3) EI = EI * EXP(X) 1301 ELSE IF (X .LE. TWO4) THEN 1302 FRAC = ZERO 1303 DO 240 I = 1, 9 1304 FRAC = Q1(I) / (P1(I) + X + FRAC) 1305 240 CONTINUE 1306 EI = (P1(10) + FRAC) / X 1307 IF (INT .NE. 3) EI = EI * EXP(X) 1308 ELSE 1309 IF ((X .GE. XMAX) .AND. (INT .LT. 3)) THEN 1310 EI = XINF 1311 ELSE 1312 Y = ONE / X 1313 FRAC = ZERO 1314 DO 250 I = 1, 9 1315 FRAC = Q2(I) / (P2(I) + X + FRAC) 1316 250 CONTINUE 1317 FRAC = P2(10) + FRAC 1318 EI = Y + Y * Y * FRAC 1319 IF (INT .NE. 3) THEN 1320 IF (X .LE. XMAX-TWO4) THEN 1321 EI = EI * EXP(X) 1322 ELSE 1323C---------------------------------------------------------------------- 1324C Calculation reformulated to avoid premature overflow 1325C---------------------------------------------------------------------- 1326 EI = (EI * EXP(X-FOURTY)) * EXP40 1327 END IF 1328 END IF 1329 END IF 1330 END IF 1331 RESULT = EI 1332 RETURN 1333 END 1334c----------------------------------------------------------------------- 1335#ifndef NWAD_PRINT 1336#define NWAD_PRINT 1337c 1338c Compile source again for Maxima 1339c 1340#include "nwxc_c_ft97.F" 1341#endif 1342#ifndef SECOND_DERIV 1343#define SECOND_DERIV 1344c 1345c Compile source again for the 2nd derivative case 1346c 1347#include "nwxc_c_ft97.F" 1348#endif 1349#ifndef THIRD_DERIV 1350#define THIRD_DERIV 1351c 1352c Compile source again for the 3rd derivative case 1353c 1354#include "nwxc_c_ft97.F" 1355#endif 1356#undef NWAD_PRINT 1357C> @} 1358