1#if defined(FUJITSU_VPP) 2!ocl scalar 3#endif 4#ifndef SECOND_DERIV 5 Subroutine xc_cams12x(tol_rho, fac, lfac, nlfac, rho, delrho, 6 & Amat, Cmat, nq, ipol, Ex, qwght,ldew,func,is12x) 7#else 8 Subroutine xc_cams12x_d2(tol_rho, fac, lfac, nlfac, rho, delrho, 9 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 10 & qwght,ldew,func,is12x) 11#endif 12c 13C$Id$ 14c 15 implicit none 16c 17#include "dft2drv.fh" 18c 19 double precision tol_rho, fac, Ex 20 integer nq, ipol, is12x 21 logical lfac, nlfac,ldew 22 double precision func(*) ! value of the functional [output] 23c 24c Charge Density 25c 26 double precision rho(nq,ipol*(ipol+1)/2) 27c 28c Charge Density Gradient 29c 30 double precision delrho(nq,3,ipol) 31c 32c Quadrature Weights 33c 34 double precision qwght(nq) 35c 36c Sampling Matrices for the XC Potential 37c 38 double precision Amat(nq,ipol), Cmat(nq,*) 39 double precision Atmp, Ctmp, Etmp 40 double precision A2tmp, C2tmp, C3tmp 41c 42#ifdef SECOND_DERIV 43c 44c Second Derivatives of the Exchange Energy Functional 45c 46 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 47#endif 48c 49 double precision rB, rC, rD, rA, rK, rE, rH, rG2, rH2 50 double precision ums, vms, dudx, dvdx, d2udx2, d2vdx2 51c 52c References: 53c 54c Swart, Chem. Phys. Lett. (2013) DOI:10.1016/j.cplett.2013.06.045. 55c 56c*************************************************************************** 57c 58 integer n 59 double precision C, rho13, rho43, gamma, x, g, gdenom, dg 60 double precision dgdenom, t, x2, x3, x4, g1, g2 61 double precision g1h1, g2h1, g1h2, g2h2, g1h3, g3h1 62 double precision hdenom, dhdenom, d2hdenom, PI, rM 63 double precision gc4, dgc4, d2gc4 64 parameter (rM=60.770665d0) 65 parameter (PI = 3.1415926535897932385d0) 66#ifdef SECOND_DERIV 67 double precision rhom23, d2g, d2gdenom, d2g1, d2g2, d2g3 68#endif 69c 70 if (is12x.eq.1) then 71c 72cswar1 1.03323556 0.00417251 0.00115216 0.75700000 0.00000000 73cswar2 0.00706184 1.20250451 0.86124355 0.00000000 0.34485046 74cswar3 1.00000000 0.00000000 0.00000000 0.00000000 1.00000000 1.52420731 75c 76 rA = 1.03323556d0 77 rK = 0.757d0 78 rC = 0.00417251d0 79 rD = 0.00115216d0 80 rE = 0.00706184d0 81 elseif (is12x.eq.2) then 82c 83cswar1 1.02149642 0.00825905 0.00235804 0.75700000 0.25000000 84cswar2 0.00654977 1.08034183 0.37999939 0.00000000 0.35897845 85cswar3 1.00000000 0.00000000 0.00000000 0.00000000 1.00000000 0.48516891 86c 87 rA = 1.02149642d0 88 rK = 0.757d0 89 rC = 0.00825905d0 90 rD = 0.00235804d0 91 rE = 0.00654977d0 92 else 93 stop 'error in xc_cams12x.F' 94 endif 95 rB = 1d0 + rK - rA 96c 97c Uniform electron gas constant 98c 99 C = -(1.5d0)*(0.75d0/acos(-1d0))**(1d0/3d0) 100c 101 if (ipol.eq.1) then 102c 103c ======> SPIN-RESTRICTED <====== 104c 105 do 10 n = 1, nq 106 if (rho(n,1).lt.tol_rho) goto 10 107c 108c Spin alpha: 109c 110 rho13 = (0.5d0*rho(n,1))**(1.d0/3.d0) 111 rho43 = rho13**4 112c 113 Etmp = 0.d0 114 Atmp = 0.d0 115 Ctmp = 0.d0 116 if (lfac) then 117 Etmp = rA * 2d0*rho43*C*fac 118 Atmp = rA * (4d0/3d0)*rho13*C*fac 119 endif 120c 121 gamma = delrho(n,1,1)*delrho(n,1,1) + 122 & delrho(n,2,1)*delrho(n,2,1) + 123 & delrho(n,3,1)*delrho(n,3,1) 124 if (dsqrt(gamma).gt.tol_rho)then 125 gamma = 0.25d0 * gamma 126 x = dsqrt(gamma) / rho43 127 x2 = x*x 128 else 129 x = 0d0 130 x2 = 0d0 131 endif 132c 133 gdenom = 1d0 + rC*x2 + rD*x2*x2 134 hdenom = 1d0 + rE*x2 135 ums = 1d0 - 1d0 / gdenom 136 vms = 1d0 - 1d0 / hdenom 137 g = C*rB*ums*vms 138c 139 dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2) 140 dvdx = 2d0*rE*x/(hdenom**2) 141 dg = C*rB*(dudx*vms + ums*dvdx) 142c 143 if (nlfac) then 144 Etmp = Etmp + 2d0*rho43*g*fac 145 Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac 146 endif 147c 148 if (x.gt.tol_rho) then 149 Ctmp = 0.5d0 * dg / sqrt(gamma) * fac 150 endif 151c 152#ifdef SECOND_DERIV 153c 154c Add local contribution back to g 155c 156 if(lfac) g = g + rA * C 157c 158 rhom23 = rho13 / (0.5d0*rho(n,1)) 159 160 d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2 161 & -20d0*rD*rD*x2*x2*x2)/(gdenom**3) 162 d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3) 163 d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2) 164c 165 A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac 166 C2tmp = - (4d0/3d0)*(rhom23**2/rho(n,1))*d2g*fac 167 if (x.gt.tol_rho) then 168 C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac 169 else 170 C3tmp = 0d0 171 endif 172c 173 call xc_att_xc_d2(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp, 174 & C2tmp,C3tmp) 175 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp 176 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp 177 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp 178#else 179 call xc_att_xc(rho(n,1),ipol,Etmp,Atmp,Ctmp) 180#endif 181 Ex = Ex + qwght(n)*Etmp 182 if (ldew) func(n) = func(n) + Etmp 183 Amat(n,1) = Amat(n,1) + Atmp 184 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp 185c 186 10 continue 187c 188 else 189c 190c ======> SPIN-UNRESTRICTED <====== 191c 192 do 20 n = 1, nq 193 if (rho(n,1).lt.tol_rho) goto 20 194 if (rho(n,2).lt.tol_rho) goto 25 195c 196c Spin alpha: 197c 198 rho13 = rho(n,2)**(1.d0/3.d0) 199 rho43 = rho13*rho(n,2) 200c 201 Etmp = 0.d0 202 Atmp = 0.d0 203 Ctmp = 0.d0 204 if (lfac) then 205 Etmp = rA * rho43*C*fac 206 Atmp = rA * (4d0/3d0)*rho13*C*fac 207 endif 208c 209 gamma = delrho(n,1,1)*delrho(n,1,1) + 210 & delrho(n,2,1)*delrho(n,2,1) + 211 & delrho(n,3,1)*delrho(n,3,1) 212 if (dsqrt(gamma).gt.tol_rho)then 213 x = dsqrt(gamma) / rho43 214 x2 = x*x 215 else 216 x = 0d0 217 x2 = 0d0 218 endif 219c 220 gdenom = 1d0 + rC*x2 + rD*x2*x2 221 hdenom = 1d0 + rE*x2 222 ums = 1d0 - 1d0 / gdenom 223 vms = 1d0 - 1d0 / hdenom 224 g = C*rB*ums*vms 225c 226 dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2) 227 dvdx = 2d0*rE*x/(hdenom**2) 228 dg = C*rB*(dudx*vms + ums*dvdx) 229c 230 if (nlfac) then 231 Etmp = Etmp + rho43*g*fac 232 Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac 233 endif 234c 235 if (x.gt.tol_rho) then 236 t = dg / sqrt(gamma) * fac 237 Ctmp = t * 0.5d0 238 endif 239c 240#ifdef SECOND_DERIV 241c 242c Add local contribution back to g 243c 244 if (lfac) g = g + rA * C 245c 246 rhom23 = rho13 / rho(n,2) 247c 248 d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2 249 & -20d0*rD*rD*x2*x2*x2)/(gdenom**3) 250 d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3) 251 d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2) 252c 253 A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac 254 C2tmp = - (2d0/3d0)*(rhom23**2/rho(n,2))*d2g*fac 255 if (x.gt.tol_rho) then 256 C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac 257 else 258 C3tmp = 0d0 259 endif 260 261 call xc_att_xc_d2(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp, 262 & C2tmp,C3tmp) 263 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp 264 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp 265 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp 266#else 267 call xc_att_xc(rho(n,2),ipol,Etmp,Atmp,Ctmp) 268#endif 269 Ex = Ex + qwght(n)*Etmp 270 if (ldew) func(n) = func(n) + Etmp 271 Amat(n,1) = Amat(n,1) + Atmp 272 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp 273c 274 25 continue 275c 276c Spin beta: 277c 278 if (rho(n,3).lt.tol_rho) goto 20 279c 280 rho13 = rho(n,3)**(1.d0/3.d0) 281 rho43 = rho13*rho(n,3) 282c 283 Etmp = 0.d0 284 Atmp = 0.d0 285 Ctmp = 0.d0 286 if (lfac) then 287 Etmp = rA * rho43*C*fac 288 Atmp = rA * (4d0/3d0)*rho13*C*fac 289 endif 290c 291 gamma = delrho(n,1,2)*delrho(n,1,2) + 292 & delrho(n,2,2)*delrho(n,2,2) + 293 & delrho(n,3,2)*delrho(n,3,2) 294 if (dsqrt(gamma).gt.tol_rho)then 295 x = dsqrt(gamma) / rho43 296 x2 = x*x 297 else 298 x = 0d0 299 x2 = 0d0 300 endif 301c 302 gdenom = 1d0 + rC*x2 + rD*x2*x2 303 hdenom = 1d0 + rE*x2 304 ums = 1d0 - 1d0 / gdenom 305 vms = 1d0 - 1d0 / hdenom 306 g = C*rB*ums*vms 307c 308 dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2) 309 dvdx = 2d0*rE*x/(hdenom**2) 310 dg = C*rB*(dudx*vms + ums*dvdx) 311c 312 if (nlfac) then 313 Etmp = Etmp + rho43*g*fac 314 Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac 315 endif 316c 317 if (x.gt.tol_rho) then 318 t = dg / sqrt(gamma) * fac 319 Ctmp = t * 0.5d0 320 endif 321c 322#ifdef SECOND_DERIV 323c 324c Add local contribution back to g 325c 326 if (lfac) g = g + rA * C 327c 328 rhom23 = rho13 / rho(n,3) 329c 330 d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2 331 & -20d0*rD*rD*x2*x2*x2)/(gdenom**3) 332 d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3) 333 d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2) 334 335c 336 A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac 337 C2tmp = -(2d0/3d0)*(rhom23**2/rho(n,3))*d2g*fac 338 if (x.gt.tol_rho) then 339 C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac 340 else 341 C3tmp = 0d0 342 endif 343c 344 call xc_att_xc_d2(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp, 345 & C2tmp,C3tmp) 346 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp 347 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + C2tmp 348 Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + C3tmp 349#else 350 call xc_att_xc(rho(n,3),ipol,Etmp,Atmp,Ctmp) 351#endif 352 Ex = Ex + qwght(n)*Etmp 353 if (ldew) func(n) = func(n) + Etmp 354 Amat(n,2) = Amat(n,2) + Atmp 355 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + Ctmp 356c 357 20 continue 358c 359 endif 360c 361 return 362 end 363#ifndef SECOND_DERIV 364#define SECOND_DERIV 365c 366c Compile source again for the 2nd derivative case 367c 368#include "xc_cams12x.F" 369#endif 370