1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 2#if !defined(NWAD_PRINT) 3C> \ingroup nwxc 4C> @{ 5C> 6C> \file nwxc_x_camlsd.F 7C> The CAM-LSD exchange functional 8C> 9C> @} 10#endif 11#endif 12C> \ingroup nwxc_priv 13C> @{ 14C> 15C> \brief Evaluate the CAM-LSD exchange functional 16C> 17C> Evaluate the CAM-LSD functional [1,2]. This routine is 18C> also used to implement CAM-B3LYP. 19C> 20C> ### References ### 21C> 22C> [1] T. Yanai, D.P. Tew, N.C. Handy, 23C> "A new hybrid exchange-correlation functional using the Coulomb-attenuating 24C> method (CAM-B3LYP)", 25C> Chem. Phys. Lett. <b>393</b>, 51-57 (2004), DOI: 26C> <a href="https://doi.org/10.1016/j.cplett.2004.06.011"> 27C> 10.1016/j.cplett.2004.06.011</a>. 28C> 29C> [2] A.D. Becke, 30C> "Density-functional exchange-energy approximation with correct 31C> asymptotic behavior", 32C> Phys. Rev. A <b>38</b>, 3098-3100 (1998), DOI: 33C> <a href="https://doi.org/10.1103/PhysRevA.38.3098"> 34C> 10.1103/PhysRevA.38.3098</a>. 35C> 36c 37c Modified to handle second derivatives while reusing code 38c 39c BGJ - 8/98 40c 41#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 42#if defined(NWAD_PRINT) 43 Subroutine nwxc_x_camlsd_p(param, tol_rho, ipol, nq, wght, rho, 44 + func) 45#else 46 Subroutine nwxc_x_camlsd(param, tol_rho, ipol, nq, wght, rho, 47 + func) 48#endif 49#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV) 50c For locations of 2nd derivatives of functionals in array 51 Subroutine nwxc_x_camlsd_d2(param, tol_rho, ipol, nq, wght, rho, 52 + func) 53#else 54 Subroutine nwxc_x_camlsd_d3(param, tol_rho, ipol, nq, wght, rho, 55 + func) 56#endif 57c 58C$Id$ 59c 60#include "nwad.fh" 61 Implicit none 62c 63#include "nwxc_param.fh" 64c 65#if defined(NWAD_PRINT) 66#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 67 type(nwad_dble)::param(*)!< [Input] Parameters of functional 68#else 69 double precision param(*)!< [Input] Parameters of functional 70#endif 71#else 72 double precision param(*)!< [Input] Parameters of functional 73 !< - param(1): \f$ \alpha_{CAM} \f$ 74 !< - param(2): \f$ \beta_{CAM} \f$ 75 !< - param(3): \f$ \omega_{CAM} \f$ 76#endif 77 double precision tol_rho !< [Input] The lower limit on the density 78 integer nq !< [Input] The number of points 79 integer ipol !< [Input] The number of spin channels 80 double precision wght !< [Input] The weight of the functional 81c 82c Charge Density 83c 84 type(nwad_dble)::rho(nq,*) !< [Input] The density 85c 86c The Exchange Energy Functional 87c 88 type(nwad_dble)::func(nq) !< [Output] The value of the functional 89c 90c Partial First Derivatives of the Exchange Energy Functional 91c 92c double precision Amat(nq,*) !< [Output] 1st order partial derivatives 93c 94#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 95c 96c Partial Second Derivatives of the Exchange Energy Functional 97c 98c double precision Amat2(nq,*) !< [Output] 2nd order partial derivatives 99#endif 100#if defined(THIRD_DERIV) 101c 102c Partial Third Order Derivatives of the Exchange Energy Functional 103c 104c double precision Amat3(nq,*) !< [Output] 3rd order partial derivatives 105#endif 106c 107c Compute the partial derivatives of the exchange functional of Dirac. 108c 109 double precision Atmp,Ctmp,A2tmp,C2tmp,C3tmp 110 type(nwad_dble):: Etmp 111 double precision A3tmp, C4tmp, C5tmp, C6tmp 112 double precision rhom23 113 double precision P1, P2, P3, P4 114c 115c P1 = -(3/PI)**(1/3) 116c P2 = -(3/4)*(3/PI)**(1/3) 117c P3 = -(6/PI)**(1/3) 118c P4 = -(3/4)*(6/PI)**(1/3) 119c 120 Parameter (P1 = -0.9847450218426959D+00) 121 Parameter (P2 = -0.7385587663820219D+00) 122 Parameter (P3 = -0.1240700981798799D+01) 123 Parameter (P4 = -0.9305257363490993D+00) 124 double precision one_third,two_ninth 125 type(nwad_dble):: rho13, rho32, rho33 126 Parameter (one_third = 1.d0/3.d0) 127 Parameter (two_ninth = 2.d0/9.d0) 128 integer n 129c 130 if (ipol.eq.1)then 131c 132c ======> SPIN-RESTRICTED <====== 133c 134 do 10 n = 1, nq 135 if (rho(n,R_T).gt.tol_rho)then 136 rho13=rho(n,R_T)**one_third 137 Etmp = rho(n,R_T)*rho13*P2*wght 138c Atmp = rho13*P1*wght 139c Ctmp = 0.d0 140#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 141c A2tmp = (rho13/rho(n,R_T))*2.0d0*one_third*P1*wght 142c C2tmp = 0.d0 143c C3tmp = 0.d0 144#endif 145#if defined(THIRD_DERIV) 146c rhom23 = rho13/rho(n,R_T) 147c A3tmp = (rhom23/rho(n,R_T))*(-4.0d0)*two_ninth*P1*wght 148c C4tmp = 0.0d0 149c C5tmp = 0.0d0 150c C6tmp = 0.0d0 151#endif 152#if defined(THIRD_DERIV) 153 call nwxc_x_att_d3(param,tol_rho,rho(n,R_T),ipol, 154 & Etmp) 155c 156c Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp 157c 158c Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + A3tmp 159#elif defined(SECOND_DERIV) 160 161 call nwxc_x_att_d2(param,tol_rho,rho(n,R_T),ipol, 162 & Etmp) 163c Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp 164#else 165#if defined(NWAD_PRINT) 166 call nwxc_x_att_p(param,tol_rho,rho(n,R_T),ipol, 167 & Etmp) 168#else 169 call nwxc_x_att(param,tol_rho,rho(n,R_T),ipol, 170 & Etmp) 171#endif 172#endif 173 func(n) = func(n) + Etmp 174c Amat(n,D1_RA) = Amat(n,D1_RA) + Atmp 175 endif 176 10 continue 177c 178 else 179c 180c ======> SPIN-UNRESTRICTED <====== 181c 182 do 20 n = 1,nq 183 if (rho(n,R_A).gt.0.5d0*tol_rho) then 184 rho32=rho(n,R_A)**one_third 185 else 186 rho32=0.0d0 187 endif 188 if (rho(n,R_B).gt.0.5d0*tol_rho) then 189 rho33=rho(n,R_B)**one_third 190 else 191 rho33=0.0d0 192 endif 193c 194 if (rho(n,R_A).gt.0.5d0*tol_rho) then 195 Etmp = rho32*rho(n,R_A)*P4*wght 196 else 197 Etmp = 0.0d0 198 endif 199c Atmp = P3*rho32*wght 200c Ctmp = 0.d0 201#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 202c A2tmp = 0.d0 203c C2tmp = 0.d0 204c C3tmp = 0.d0 205c if (rho(n,R_A).gt.0.5d0*tol_rho) then 206c A2tmp = one_third*P3*rho32/rho(n,R_A)*wght 207c end if 208#endif 209#if defined(THIRD_DERIV) 210c A3tmp = 0.0d0 211c C4tmp = 0.0d0 212c C5tmp = 0.0d0 213c C6tmp = 0.0d0 214c 215c if (rho(n,R_A).gt.0.5d0*tol_rho) then 216c A3tmp = -two_ninth*P3*rho32/(rho(n,R_A)**2)*wght 217c endif 218#endif 219#if defined(THIRD_DERIV) 220 if (rho(n,R_A).gt.0.5d0*tol_rho) then 221 call nwxc_x_att_d3(param,tol_rho,rho(n,R_A),ipol, 222 & Etmp) 223 endif 224c 225c Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp 226c 227c Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) + A3tmp 228#elif defined(SECOND_DERIV) 229 if (rho(n,R_A).gt.0.5d0*tol_rho) then 230 call nwxc_x_att_d2(param,tol_rho,rho(n,R_A),ipol, 231 & Etmp) 232 end if 233c Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp 234#else 235#if defined(NWAD_PRINT) 236 if (rho(n,R_A).gt.0.5d0*tol_rho) then 237 call nwxc_x_att_p(param,tol_rho,rho(n,R_A),ipol, 238 & Etmp) 239 endif 240#else 241 if (rho(n,R_A).gt.0.5d0*tol_rho) then 242 call nwxc_x_att(param,tol_rho,rho(n,R_A),ipol, 243 & Etmp) 244 endif 245#endif 246#endif 247 func(n) = func(n) + Etmp 248c Amat(n,D1_RA) = Amat(n,D1_RA) + Atmp 249c 250c Beta spin channel 251c 252 if (rho(n,R_B).gt.0.5d0*tol_rho) then 253 Etmp = rho33*rho(n,R_B)*P4*wght 254 else 255 Etmp = 0.0d0 256 endif 257c Atmp = P3*rho33*wght 258c Ctmp = 0.d0 259#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 260c A2tmp = 0.d0 261c C2tmp = 0.d0 262c C3tmp = 0.d0 263c if (rho(n,R_B).gt.0.5d0*tol_rho) then 264c A2tmp = one_third*P3*rho33/rho(n,R_B)*wght 265c end if 266#endif 267#if defined(THIRD_DERIV) 268c A3tmp = 0.0d0 269c C4tmp = 0.0d0 270c C5tmp = 0.0d0 271c C6tmp = 0.0d0 272c 273c if (rho(n,R_B).gt.0.5d0*tol_rho) then 274c A3tmp = -two_ninth*P3*rho33/(rho(n,R_B)**2)*wght 275c endif 276#endif 277#if defined(THIRD_DERIV) 278 if (rho(n,R_B).gt.0.5d0*tol_rho) then 279 call nwxc_x_att_d3(param,tol_rho,rho(n,R_B),ipol, 280 & Etmp) 281 endif 282c 283c Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp 284c 285c Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) + A3tmp 286#elif defined(SECOND_DERIV) 287 if (rho(n,R_B).gt.0.5d0*tol_rho) then 288 call nwxc_x_att_d2(param,tol_rho,rho(n,R_B),ipol, 289 & Etmp) 290 end if 291c Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp 292#else 293#if defined(NWAD_PRINT) 294 if (rho(n,R_B).gt.0.5d0*tol_rho) then 295 call nwxc_x_att_p(param,tol_rho,rho(n,R_B),ipol, 296 & Etmp) 297 endif 298#else 299 if (rho(n,R_B).gt.0.5d0*tol_rho) then 300 call nwxc_x_att(param,tol_rho,rho(n,R_B),ipol, 301 & Etmp) 302 endif 303#endif 304#endif 305 func(n) = func(n) + Etmp 306c Amat(n,D1_RB) = Amat(n,D1_RB) + Atmp 307c 308c func(n) = func(n) + ( rho32*rho(n,R_A) + 309c & rho33*rho(n,R_B) )*P4*wght 310 20 continue 311c 312 endif 313c 314 return 315 end 316#ifndef NWAD_PRINT 317#define NWAD_PRINT 318c 319c Compile source again for the 2nd derivative case 320c 321#include "nwxc_x_camlsd.F" 322#endif 323#ifndef SECOND_DERIV 324#define SECOND_DERIV 325c 326c Compile source again for the 2nd derivative case 327c 328#include "nwxc_x_camlsd.F" 329#endif 330#ifndef THIRD_DERIV 331#define THIRD_DERIV 332c 333c Compile source again for the 3rd derivative case 334c 335#include "nwxc_x_camlsd.F" 336#endif 337#undef NWAD_PRINT 338C> @} 339