1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2#if !defined(NWAD_PRINT) 3C> \ingroup nwxc 4C> @{ 5C> 6C> \file nwxc_c_op.F 7C> The OP correlation functional 8C> 9C> @} 10#endif 11#endif 12C> 13C> \ingroup nwxc_priv 14C> @{ 15C> 16C> \brief Evaluate the OP correlation functional 17C> 18C> The OP correlation functional [1,2] is a functional designed to 19C> have few optimized parameters (only one in this case) and has 20C> a "progressive" form. 21C> 22C> ### References ### 23C> 24C> [1] T. Tsuneda, T. Suzumura, K. Hirao, 25C> "A new one-parameter progressive Colle-Salvetti-type correlation 26C> functional", J. Chem. Phys. <b>110</b>, 10664-10678 (1999), DOI: 27C> <a href="https://doi.org/10.1063/1.479012"> 28C> 10.1063/1.479012</a>. 29C> 30C> [2] T. Tsuneda, T. Suzumura, K. Hirao, 31C> "A reexamination of exchange energy functionals", 32C> J. Chem. Phys. <b>111</b>, 5656-5667 (1999), DOI: 33C> <a href="https://doi.org/10.1063/1.479954"> 34C> 10.1063/1.479954</a>. 35C> 36#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 37#if defined(NWAD_PRINT) 38 Subroutine nwxc_c_op_p(kop,param,tol_rho,ipol,nq,wght,rho,rgamma, 39 & func) 40#else 41 Subroutine nwxc_c_op(kop,param,tol_rho,ipol,nq,wght,rho,rgamma, 42 & func) 43#endif 44#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 45 Subroutine nwxc_c_op_d2(kop,param,tol_rho,ipol,nq,wght,rho,rgamma, 46 & func) 47#else 48 Subroutine nwxc_c_op_d3(kop,param,tol_rho,ipol,nq,wght,rho,rgamma, 49 & func) 50#endif 51c 52C$Id$ 53c 54#include "nwad.fh" 55c 56 implicit none 57c 58#include "nwxc_param.fh" 59c 60c Input and other parameters 61c 62 external kop !< [Input] Subroutine to evaluate the 63 !< GGA exchange enhancement factor (see Eq.(14) in [1]). 64#if defined(NWAD_PRINT) 65#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 66 type(nwad_dble)::param(*)!< [Input] Parameters of functional 67 type(nwad_dble)::QABOP 68#else 69 double precision param(*)!< [Input] Parameters of functional 70 double precision QABOP 71#endif 72#else 73 double precision param(*)!< [Input] Parameters of functional 74 !< - param(1): \f$ q^{\alpha\beta}_{OP} \f$ see Eq.(27) in [1] and 75 !< Table III in [2]. 76 double precision QABOP 77#endif 78 double precision tol_rho !< [Input] The lower limit on the density 79 integer ipol !< [Input] The number of spin channels 80 integer nq !< [Input] The number of points 81 double precision wght !< [Input] The weight of the functional 82c 83c Charge Density 84c 85 type(nwad_dble)::rho(nq,*) !< [Input] The density 86c 87c Charge Density Gradient 88c 89 type(nwad_dble)::rgamma(nq,*) !< [Input] The norm of the density gradients 90c 91c Sampling Matrices for the XC Potential 92c 93 type(nwad_dble)::func(nq) !< [Output] The value of the functional 94c double precision Amat(nq,*) !< [Output] The derivative wrt rho 95c double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma 96c 97 double precision QAB88OP,QABPBOP 98c Parameter (QAB88OP=2.3670D0,QABPBOP=2.3789D0) 99c 100c References: 101c Tsuneda, Suzumura, Hirao, JCP 110, 10664 (1999) 102c Tsuneda, Suzumura, Hirao, JCP 111, 5656 (1999) 103c 104c*************************************************************************** 105c 106 integer n 107 type(nwad_dble)::rho13, rho43, gamma, x 108 type(nwad_dble)::kalpha,kbeta, rho13a, rho13b,rhoa,rhob 109c type(nwad_dble)::banb, hbab, hbabx 110 type(nwad_dble)::banb, hbab 111 double precision dhdab,dhdabx,dkadra,dkbdrb,dkadxa,dkbdxb, 112 A dbabdra,dbabdrb,dbabdga,dbabdgb,dkadga,dkbdgb, 113 A dbabdka,dbabdkb 114c 115c hbabx(x) = (1.5214d0*x + 0.5764d0)/ 116c / (x**2.0d0*(x**2.0d0+1.1284d0*x+0.3183d0)) 117c dhdabx(x) = -(4.5642d0*x**4+5.7391d0*x**3+ 118c + 2.4355*x**2+0.3669d0*x)/ 119c / ((x**4+1.1284d0*x**3+0.3183d0*x**2)**2) 120c 121c if(whichf.eq.'be88') then 122c QABOP=QAB88OP 123c endif 124c if(whichf.eq.'pb96') then 125c QABOP=QABPBOP 126c endif 127 QABOP = param(1) 128 if (ipol.eq.1) then 129c 130c ======> SPIN-RESTRICTED <====== 131c 132 do 10 n = 1, nq 133 if (rho(n,R_T).lt.tol_rho) goto 10 134c 135c Spin alpha: 136c 137 rhoa=rho(n,R_T)*0.5d0 138 rho13a = (rhoa)**(1.d0/3.d0) 139 rho43 = rho13a**4.0d0 140 gamma = rgamma(n,G_TT) 141c gamma = delrho(n,1,1)*delrho(n,1,1) + 142c & delrho(n,2,1)*delrho(n,2,1) + 143c & delrho(n,3,1)*delrho(n,3,1) 144 gamma = 0.25d0 * gamma 145 if (sqrt(gamma).gt.tol_rho)then 146 x = sqrt(gamma) / rho43 147 call kop(tol_rho,x,kalpha) 148c dkadra = -(4d0/3d0)*x*dkadxa/rhoa 149c dkadga = (dkadxa/rho43)*0.5d0/sqrt(gamma) 150 else 151 x=0d0 152 call kop(tol_rho,x,kalpha) 153c dkadra = 0d0 154c dkadga = 0d0 155 endif 156c 157c 158 159 banb = qabop * rho13a * kalpha *0.5d0 160 161 if(banb.ne.0.0d0) then 162c dbabdra = banb*0.5d0* 163c / (1d0/(3d0*rhoa)+dkadra/kalpha) 164 165c dbabdga = banb/kalpha*dkadga*0.5d0 166 167 hbab = hbabx(banb) 168c dhdab = dhdabx(banb) 169 else 170c dbabdra =0d0 171c dbabdga =0d0 172 hbab = 0d0 173c dhdab = 0d0 174 endif 175 176c Ec = Ec - rhoa**2*hbab*qwght(n)*fac 177 func(n) = func(n) - rhoa**2.0d0*hbab*wght 178c Amat(n,D1_RA) = Amat(n,D1_RA) - 179c - (rhoa*hbab + rhoa**2*dhdab*dbabdra)*wght 180 181c 182c if (x.gt.tol_rho) then 183c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) - 184c - rhoa**2*dhdab*dbabdga*wght 185c endif 186c 187 10 continue 188c 189 else 190c 191c ======> SPIN-UNRESTRICTED <====== 192c 193 do 20 n = 1, nq 194 if (rho(n,R_A)+rho(n,R_B).lt.tol_rho) goto 20 195 if (rho(n,R_A).ge.tol_rho*0.5d0) then 196c 197c Spin alpha: 198c 199 rhoa=rho(n,R_A) 200 rho13a = rhoa**(1.d0/3.d0) 201 rho43 = rho13a*rhoa 202 gamma = rgamma(n,G_AA) 203c gamma = delrho(n,1,1)*delrho(n,1,1) + 204c & delrho(n,2,1)*delrho(n,2,1) + 205c & delrho(n,3,1)*delrho(n,3,1) 206 if (sqrt(gamma).gt.tol_rho)then 207 x = sqrt(gamma) / rho43 208 call kop(tol_rho,x,kalpha) 209c dkadra = -(4d0/3d0)*x*dkadxa/rhoa 210c dkadga = dkadxa*0.5d0/(rho43*dsqrt(gamma)) 211 else 212 x = 0d0 213 endif 214 else 215 rhoa=0d0 216 rho13a=0d0 217 x = 0d0 218 endif 219 if(x.eq.0d0) then 220 call kop(tol_rho,x,kalpha) 221c dkadra = 0d0 222c dkadga = 0d0 223 endif 224c 225c Spin beta: 226c 227 if (rho(n,R_B).ge.tol_rho*0.5d0) then 228c 229 rhob=rho(n,R_B) 230 rho13b = rhob**(1.d0/3.d0) 231 rho43 = rho13b*rhob 232 gamma = rgamma(n,G_BB) 233c gamma = delrho(n,1,2)*delrho(n,1,2) + 234c & delrho(n,2,2)*delrho(n,2,2) + 235c & delrho(n,3,2)*delrho(n,3,2) 236 if (sqrt(gamma).gt.tol_rho)then 237 x = sqrt(gamma) / rho43 238 call kop(tol_rho,x, kbeta) 239c dkbdrb = -(4d0/3d0)*x*dkbdxb/rhob 240c dkbdgb = dkbdxb*0.5d0/(rho43*dsqrt(gamma)) 241 else 242 x = 0d0 243 endif 244 else 245 if(rho13a.eq.0.0d0) goto 20 246 rhob=0d0 247 rho13b=0d0 248 x=0d0 249 endif 250 if(x.eq.0d0) then 251 call kop(tol_rho,x, kbeta) 252c dkbdrb = 0d0 253c dkbdgb= 0d0 254 endif 255 256 banb = qabop*(rho13a*kalpha*rho13b*kbeta)/ 257 / (rho13a*kalpha+rho13b*kbeta) 258 259 if(banb.ne.0.0d0) then 260c dbabdra = banb*kbeta*rho13b/ 261c / (rho13a*kalpha+rho13b*kbeta)* 262c / (1d0/(3d0*rhoa)+dkadra/kalpha) 263c dbabdrb = banb*kalpha*rho13a/ 264c / (rho13a*kalpha+rho13b*kbeta)* 265c / (1d0/(3d0*rhob)+dkbdrb/kbeta) 266 267c dbabdga = banb*rho13b*kbeta/ 268c / ((rho13a*kalpha+rho13b*kbeta)*kalpha)* 269c * dkadga 270c dbabdgb = banb*rho13a*kalpha/ 271c / ((rho13a*kalpha+rho13b*kbeta)*kbeta)* 272c * dkbdgb 273 274 hbab = hbabx(banb) 275c dhdab = dhdabx(banb) 276 else 277c dbabdra =0d0 278c dbabdrb =0d0 279c dbabdga =0d0 280c dbabdgb =0d0 281 hbab = 0d0 282c dhdab = 0d0 283 endif 284 285c Ec = Ec - rhoa*rhob*hbab*qwght(n)*fac 286 func(n) = func(n) - rhoa*rhob*hbab*wght 287c Amat(n,D1_RA) = Amat(n,D1_RA) - 288c - (rhob*hbab + rhoa*rhob*dhdab*dbabdra)*wght 289c Amat(n,D1_RB) = Amat(n,D1_RB) - 290c - (rhoa*hbab + rhoa*rhob*dhdab*dbabdrb)*wght 291c 292c 293c if (x.gt.tol_rho) then 294c Cmat(n,D1_GAA) = Cmat(n,D1_GAA) - 295c - rhoa*rhob*dhdab*dbabdga*wght 296c Cmat(n,D1_GBB) = Cmat(n,D1_GBB) - 297c - rhoa*rhob*dhdab*dbabdgb*wght 298c endif 299 300c 301c 302 20 continue 303c 304 endif 305c 306 return 307c 308 contains 309c 310c The combination of statement functions and derived types with 311c overloaded operators is not properly supported in the Fortran 312c standard (apparently). Therefore the statement functions from 313c the original subroutine had to be transformed into contained 314c functions. 315c 316c WARNING: This code is EXTREMELY EVIL! Although there appears to be 317c only one function there are actually three with the same name, 318c each one returning results of a different data type. The trick is 319c that depending on the data type the the subroutine that contains 320c these functions changes its name and hence these different 321c functions of the same name do not lead to conflicts. The 322c alternative would have been to add a forest of conditional 323c compilation constructs (#ifdef's) to change the function names 324c in the declarations and all places where they are used. That 325c would have been extremely ugly, so we are between a rock and a 326c hard place on this one. 327c 328 function hbabx(x) result(s) 329#include "nwad.fh" 330 implicit none 331 type(nwad_dble), intent(in) :: x 332 type(nwad_dble) :: s 333 s = (1.5214d0*x + 0.5764d0)/ 334 . (x**2.0d0*(x**2.0d0+1.1284d0*x+0.3183d0)) 335 end function 336c 337 end 338C> 339C> \brief The Becke88 exchange GGA enhancement factor 340C> 341#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 342#if defined(NWAD_PRINT) 343 subroutine nwxc_k_becke88_p(tol_rho,x, kalpha) 344#else 345 subroutine nwxc_k_becke88(tol_rho,x, kalpha) 346#endif 347#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 348 subroutine nwxc_k_becke88_d2(tol_rho,x, kalpha) 349#else 350 subroutine nwxc_k_becke88_d3(tol_rho,x, kalpha) 351#endif 352#include "nwad.fh" 353 implicit none 354c 355 double precision tol_rho 356 type(nwad_dble)::x,kalpha 357 double precision dkadxa 358c 359 double precision BETA, C 360 Parameter (BETA = 0.0042D0) 361 type(nwad_dble)::g,gdenom 362 double precision dgdenom,dg 363c double precision arcsinh, darcsinh 364c arcsinh(x)=log(x+dsqrt(1d0+x*x)) 365c darcsinh(x)=1d0/dsqrt(1d0+x*x) 366c 367c 368c Uniform electron gas constant 369c 370 C = 3d0*(0.75d0/acos(-1d0))**(1d0/3d0) 371 372 if (x.gt.0d0)then 373 gdenom = 1d0 + 6d0*BETA*x*asinh(x) 374 g = 2d0*BETA*x*x / gdenom 375c dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x)) 376c dg = g*(2d0/x-dgdenom/gdenom) 377 378 kalpha= C + g 379c dkadxa = dg 380 381 else 382 kalpha= C 383c dkadxa = 0d0 384 endif 385 return 386 end 387C> 388C> \brief The PBE96 exchange GGA enhancement factor 389C> 390#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 391#if defined(NWAD_PRINT) 392 subroutine nwxc_k_pbe96_p(tol_rho,x, kalpha) 393#else 394 subroutine nwxc_k_pbe96(tol_rho,x, kalpha) 395#endif 396#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 397 subroutine nwxc_k_pbe96_d2(tol_rho,x, kalpha) 398#else 399 subroutine nwxc_k_pbe96_d3(tol_rho,x, kalpha) 400#endif 401#include "nwad.fh" 402 implicit none 403c 404 double precision tol_rho 405 type(nwad_dble)::x,kalpha,deno 406 double precision dkadxa 407c 408 double precision pi,um, uk, umk 409 parameter(um=0.21951d0, uk=0.804d0, umk=um/uk) 410 double precision C 411 double precision forty8 412c 413c 414c Uniform electron gas constant 415c 416 pi = acos(-1.d0) 417 C = 3d0*(0.75d0/pi)**(1d0/3d0) 418 419 if (x.gt.0d0)then 420 forty8=1d0/((48d0*pi*pi)**(2d0/3d0)) 421 deno=1d0/(1d0+um*x*x*forty8/uk) 422 kalpha= C * (1d0 + uk - uk *deno) 423c dkadxa = C * (2d0*um*x*deno*deno* 424c * forty8) 425 426 else 427 kalpha= C 428c dkadxa = 0d0 429 endif 430 return 431 end 432C> 433C> \brief The Dirac exchange GGA enhancement factor 434C> 435C> Of course the Dirac functional is the exchange part for LDA 436C> so this subroutine just returns a constant. 437C> 438#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 439#if defined(NWAD_PRINT) 440 subroutine nwxc_k_dirac_p(tol_rho,x, kalpha) 441#else 442 subroutine nwxc_k_dirac(tol_rho,x, kalpha) 443#endif 444#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 445 subroutine nwxc_k_dirac_d2(tol_rho,x, kalpha) 446#else 447 subroutine nwxc_k_dirac_d3(tol_rho,x, kalpha) 448#endif 449#include "nwad.fh" 450 implicit none 451c 452 double precision tol_rho 453 type(nwad_dble)::x,kalpha 454 double precision dkadxa 455c 456 double precision pi 457 double precision C 458c 459c 460c Uniform electron gas constant 461c 462 pi = acos(-1.d0) 463 C = 3d0*(0.75d0/pi)**(1d0/3d0) 464 465 kalpha= C 466c dkadxa = 0d0 467 468 return 469 end 470#ifndef NWAD_PRINT 471#define NWAD_PRINT 472c 473c Compile source again for the 2nd derivative case 474c 475#include "nwxc_c_op.F" 476#endif 477#ifndef SECOND_DERIV 478#define SECOND_DERIV 479c 480c Compile source again for the 2nd derivative case 481c 482#include "nwxc_c_op.F" 483#endif 484#ifndef THIRD_DERIV 485#define THIRD_DERIV 486c 487c Compile source again for the 3rd derivative case 488c 489#include "nwxc_c_op.F" 490#endif 491#undef NWAD_PRINT 492C> 493C> @} 494