1#if defined(FUJITSU_VPP) 2!ocl scalar 3#endif 4#ifndef SECOND_DERIV 5 Subroutine xc_s12x(tol_rho, fac, lfac, nlfac, rho, delrho, 6 & Amat, Cmat, nq, ipol, Ex, qwght,ldew,func,is12x) 7#else 8 Subroutine xc_s12x_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,*) 39c 40#ifdef SECOND_DERIV 41c 42c Second Derivatives of the Exchange Energy Functional 43c 44 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 45#endif 46c 47 double precision rB, rC, rD, rA, rK, rE, rH, rG2, rH2 48 double precision ums, vms, dudx, dvdx, d2udx2, d2vdx2 49c 50c References: 51c 52c Swart, Chem. Phys. Lett. (2013), DOI:10.1016/j.cplett.2013.06.045. 53c 54c*************************************************************************** 55c 56 integer n 57 double precision C, rho13, rho43, gamma, x, g, gdenom, dg 58 double precision dgdenom, t, x2, x3, x4, g1, g2 59 double precision g1h1, g2h1, g1h2, g2h2, g1h3, g3h1 60 double precision hdenom, dhdenom, d2hdenom, PI, rM 61 double precision gc4, dgc4, d2gc4 62 parameter (rM=60.770665d0) 63 parameter (PI = 3.1415926535897932385d0) 64#ifdef SECOND_DERIV 65 double precision rhom23, d2g, d2gdenom, d2g1, d2g2, d2g3 66#endif 67c 68 69 if (is12x.eq.1) then 70c 71cswar1 1.03842032 0.00403198 0.00104596 0.75700000 0.00000000 72cswar2 0.00594635 1.17755954 0.84432515 0.00000000 0.00000000 73cswar3 1.00000000 0.00000000 0.00000000 0.00000000 1.00000000 0.00000000 74c 75 rA = 1.03842032d0 76 rK = 0.757d0 77 rC = 0.00403198d0 78 rD = 0.00104596d0 79 rE = 0.00594635d0 80 elseif (is12x.eq.2) then 81c 82cswar1 1.02543951 0.00761554 0.00211063 0.75700000 0.25000000 83cswar2 0.00604672 1.07735222 0.37705816 0.00000000 0.00000000 84cswar3 1.00000000 0.00000000 0.00000000 0.00000000 1.00000000 0.00000000 85c 86 rA = 1.02543951d0 87 rK = 0.757d0 88 rC = 0.00761554d0 89 rD = 0.00211063d0 90 rE = 0.00604672d0 91 else 92 stop 'error in xc_s12x.F' 93 endif 94 rB = 1d0 + rK - rA 95c 96c Uniform electron gas constant 97c 98 C = -(1.5d0)*(0.75d0/acos(-1d0))**(1d0/3d0) 99c 100 if (ipol.eq.1) then 101c 102c ======> SPIN-RESTRICTED <====== 103c 104 do 10 n = 1, nq 105 if (rho(n,1).lt.tol_rho) goto 10 106c 107c Spin alpha: 108c 109 rho13 = (0.5d0*rho(n,1))**(1.d0/3.d0) 110 rho43 = rho13**4 111 gamma = delrho(n,1,1)*delrho(n,1,1) + 112 & delrho(n,2,1)*delrho(n,2,1) + 113 & delrho(n,3,1)*delrho(n,3,1) 114 if (dsqrt(gamma).gt.tol_rho)then 115 gamma = 0.25d0 * gamma 116 x = dsqrt(gamma) / rho43 117 x2 = x*x 118 else 119 x = 0d0 120 x2 = 0d0 121 endif 122c 123 gdenom = 1d0 + rC*x2 + rD*x2*x2 124 hdenom = 1d0 + rE*x2 125 ums = 1d0 - 1d0 / gdenom 126 vms = 1d0 - 1d0 / hdenom 127 g = C*rB*ums*vms 128c 129 dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2) 130 dvdx = 2d0*rE*x/(hdenom**2) 131 dg = C*rB*(dudx*vms + ums*dvdx) 132c 133cmswart if (lfac) then 134cmswart Ex = Ex + rA*2d0*rho43*C*qwght(n)*fac 135cmswart if(ldew)func(n) = func(n) + rA*2.d0*rho43*C*fac 136cmswart Amat(n,1) = Amat(n,1) + rA*(4d0/3d0)*rho13*C*fac 137cmswart endif 138c 139 if (nlfac) then 140 Ex = Ex + 2d0*rho43*g*qwght(n)*fac 141 if(ldew)func(n) = func(n) + 2.d0*rho43*g*fac 142 Amat(n,1) = Amat(n,1) + (4d0/3d0)*rho13*(g-x*dg)*fac 143 endif 144c 145 if (x.gt.tol_rho) then 146 t = 0.5d0 * dg / sqrt(gamma) * fac 147 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + t 148 endif 149c 150#ifdef SECOND_DERIV 151cmswart if(lfac) g = g + rA*C ! Add local contribution back to g 152 rhom23 = rho13 / (0.5d0*rho(n,1)) 153 154 d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2 155 & -20d0*rD*rD*x2*x2*x2)/(gdenom**3) 156 d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3) 157 d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2) 158c 159 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 160 & + (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac 161 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 162 & - (4d0/3d0)*(rhom23**2/rho(n,1))*d2g*fac 163 if (x.gt.tol_rho) then 164 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 165 & - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac 166 endif 167#endif 168c 169 10 continue 170c 171 else 172c 173c ======> SPIN-UNRESTRICTED <====== 174c 175 do 20 n = 1, nq 176 if (rho(n,1).lt.tol_rho) goto 20 177 if (rho(n,2).lt.tol_rho) goto 25 178c 179c Spin alpha: 180c 181 rho13 = rho(n,2)**(1.d0/3.d0) 182 rho43 = rho13*rho(n,2) 183 gamma = delrho(n,1,1)*delrho(n,1,1) + 184 & delrho(n,2,1)*delrho(n,2,1) + 185 & delrho(n,3,1)*delrho(n,3,1) 186 if (dsqrt(gamma).gt.tol_rho)then 187 x = dsqrt(gamma) / rho43 188 x2 = x*x 189 else 190 x = 0d0 191 x2 = 0d0 192 endif 193c 194 gdenom = 1d0 + rC*x2 + rD*x2*x2 195 hdenom = 1d0 + rE*x2 196 ums = 1d0 - 1d0 / gdenom 197 vms = 1d0 - 1d0 / hdenom 198 g = C*rB*ums*vms 199c 200 dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2) 201 dvdx = 2d0*rE*x/(hdenom**2) 202 dg = C*rB*(dudx*vms + ums*dvdx) 203c 204cmswart if (lfac) then 205cmswart Ex = Ex + rA*rho43*C*qwght(n)*fac 206cmswart if (ldew)func(n) = func(n) + rA*rho43*C*fac 207cmswart Amat(n,1) = Amat(n,1) + rA*(4d0/3d0)*rho13*C*fac 208cmswart endif 209c 210 if (nlfac) then 211 Ex = Ex + rho43*g*qwght(n)*fac 212 if (ldew)func(n) = func(n) + rho43*g*fac 213 Amat(n,1) = Amat(n,1) + (4d0/3d0)*rho13*(g-x*dg)*fac 214 endif 215c 216 if (x.gt.tol_rho) then 217 t = dg / sqrt(gamma) * fac 218 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + t * 0.5d0 219 endif 220c 221#ifdef SECOND_DERIV 222cmswart if (lfac) g = g + rA*C ! Add local contribution back to g 223 rhom23 = rho13 / rho(n,2) 224 225 d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2 226 & -20d0*rD*rD*x2*x2*x2)/(gdenom**3) 227 d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3) 228 d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2) 229c 230 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 231 & + (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac 232 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 233 & - (2d0/3d0)*(rhom23**2/rho(n,2))*d2g*fac 234 if (x.gt.tol_rho) then 235 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 236 & - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac 237 endif 238#endif 239c 240 25 continue 241c 242c Spin beta: 243c 244 if (rho(n,3).lt.tol_rho) goto 20 245c 246 rho13 = rho(n,3)**(1.d0/3.d0) 247 rho43 = rho13*rho(n,3) 248 gamma = delrho(n,1,2)*delrho(n,1,2) + 249 & delrho(n,2,2)*delrho(n,2,2) + 250 & delrho(n,3,2)*delrho(n,3,2) 251 if (dsqrt(gamma).gt.tol_rho)then 252 x = dsqrt(gamma) / rho43 253 x2 = x*x 254 else 255 x = 0d0 256 x2 = 0d0 257 endif 258c 259 gdenom = 1d0 + rC*x2 + rD*x2*x2 260 hdenom = 1d0 + rE*x2 261 ums = 1d0 - 1d0 / gdenom 262 vms = 1d0 - 1d0 / hdenom 263 g = C*rB*ums*vms 264c 265 dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2) 266 dvdx = 2d0*rE*x/(hdenom**2) 267 dg = C*rB*(dudx*vms + ums*dvdx) 268c 269cmswart if (lfac) then 270cmswart Ex = Ex + rA*rho43*C*qwght(n)*fac 271cmswart if (ldew)func(n) = func(n) + rA*rho43*C*fac 272cmswart Amat(n,2) = Amat(n,2) + rA*(4d0/3d0)*rho13*C*fac 273cmswart endif 274c 275 if (nlfac) then 276 Ex = Ex + rho43*g*qwght(n)*fac 277 if (ldew)func(n) = func(n) +rho43*g*fac 278 Amat(n,2) = Amat(n,2) + (4d0/3d0)*rho13*(g-x*dg)*fac 279 endif 280c 281 if (x.gt.tol_rho) then 282 t = dg / sqrt(gamma) * fac 283 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + t * 0.5d0 284 endif 285c 286#ifdef SECOND_DERIV 287cmswart if(lfac) g = g + rA*C ! Add local contribution back to g 288 rhom23 = rho13 / rho(n,3) 289 290 d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2 291 & -20d0*rD*rD*x2*x2*x2)/(gdenom**3) 292 d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3) 293 d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2) 294 295c 296 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 297 & + (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac 298 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) 299 & - (2d0/3d0)*(rhom23**2/rho(n,3))*d2g*fac 300 if (x.gt.tol_rho) then 301 Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) 302 & - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac 303 endif 304#endif 305c 306 20 continue 307c 308 endif 309c 310 return 311 end 312#ifndef SECOND_DERIV 313#define SECOND_DERIV 314c 315c Compile source again for the 2nd derivative case 316c 317#include "xc_s12x.F" 318#endif 319