1c SSB-D exchange functional part 1 2c (the one that depends on s) 3c 4c References: 5c [a] J.P. Perdew, K. Burke, and M. Ernzerhof, PRL 77, 3865 (1996). 6c [b] M. Swart, M. Sola, and F.M. Bickelhaupt, JCP 131, 094103 (2009). 7c 8#ifndef SECOND_DERIV 9 Subroutine xc_ssbD_1(tol_rho, fac, lfac, nlfac, rho, delrho, 10 & Amat, Cmat, nq, ipol, Ex, qwght,ldew,func) 11#else 12 Subroutine xc_ssbD_1_d2(tol_rho, fac, lfac, nlfac, rho, delrho, 13 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 14 & qwght,ldew,func) 15#endif 16c 17c$Id$ 18c 19 implicit none 20c 21#include "dft2drv.fh" 22c 23 double precision fac, Ex 24 integer nq, ipol 25 logical lfac, nlfac,ldew 26 double precision func(*) ! value of the functional [output] 27c 28c Charge Density & Its Cube Root 29c 30 double precision rho(nq,ipol*(ipol+1)/2) 31c 32c Charge Density Gradient 33c 34 double precision delrho(nq,3,ipol) 35c 36c Quadrature Weights 37c 38 double precision qwght(nq) 39c 40c Sampling Matrices for the XC Potential & Energy 41c 42 double precision amat(nq,ipol), cmat(nq,*) 43#ifdef SECOND_DERIV 44 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 45#endif 46c 47 double precision tol_rho, pi 48 double precision rA, rB, rC, rD, rE, rU 49 double precision C, Cs 50 double precision F43, F13 51#ifdef SECOND_DERIV 52 double precision F73 53#endif 54 parameter (rA=1.079966d0, rB=0.197465d0, rC=0.272729d0) 55 parameter (rE=5.873645d0, rU=-0.749940d0) 56 parameter (rD=rB*(1.0d0-rU)) 57c 58 parameter (F43=4.d0/3.d0, F13=1.d0/3.d0) 59#ifdef SECOND_DERIV 60 parameter (F73=7.d0/3.d0) 61#endif 62c 63 integer n 64 double precision rrho, rho43, rho13, gamma, gam12, s, d1s(2), 65 & d, g, gp, d1g(2), facSSB 66#ifdef SECOND_DERIV 67 double precision rhom23, d2s(3), gpp, d2g(3), gssb2 68#endif 69 double precision gssb0,gssb1 70 gssb0(s)= rB*s*s/(1d0+rC*s*s) 71 + - rD*s*s/(1d0+rE*s**4) 72 gssb1(s)= 2d0*rB*s/(1d0+rC*s*s)**2 + 73 + (2d0*rD*rE*s**5 - 2d0*rD*s)/(1d0+rE*s**4)**2 74#ifdef SECOND_DERIV 75 gssb2(s)= 8d0*rB/(1d0+rC*s*s)**3 - 6d0*rB/(1d0+rC*s*s)**2 76 + + 36d0*rD/(1d0+rE*s**4)**2 - 32d0*rD/(1d0+rE*s**4)**3 77 - - 6d0*rD/(1d0+rE*s**4) 78#endif 79c 80 pi = acos(-1.d0) 81 C = -3d0/(4d0*pi)*(3d0*pi*pi)**F13 82 Cs = 0.5d0/(3d0*pi*pi)**F13 83 Cs = Cs * C ! account for including C in rho43 84c 85 if (ipol.eq.1 )then 86c 87c ======> SPIN-RESTRICTED <====== 88c 89#ifdef IFCV81 90CDEC$ NOSWP 91#endif 92 do 10 n = 1, nq 93 if (rho(n,1).lt.tol_rho) goto 10 94 rho43 = C*rho(n,1)**F43 95 rrho = 1d0/rho(n,1) 96 rho13 = F43*rho43*rrho 97#ifdef SECOND_DERIV 98 rhom23 = F13*rho13*rrho 99#endif 100c 101 gamma = delrho(n,1,1)*delrho(n,1,1) + 102 & delrho(n,2,1)*delrho(n,2,1) + 103 & delrho(n,3,1)*delrho(n,3,1) 104 gam12 = dsqrt(gamma) 105 if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 10 106c 107 108 s = Cs*gam12/rho43 109 d1s(1) = -F43*s*rrho 110 d1s(2) = 0.5d0*s/gamma 111c 112c Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1 113c 114 g=gssb0(s) 115 gp=gssb1(s) 116c 117 d1g(1) = gp*d1s(1) 118 d1g(2) = gp*d1s(2) 119 Ex = Ex + rho43*g*qwght(n)*fac 120 if(ldew)func(n) = func(n) + rho43*g*fac 121 Amat(n,1) = Amat(n,1) + (rho13*g+rho43*d1g(1))*fac 122 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 2d0*rho43*d1g(2)*fac 123#ifdef SECOND_DERIV 124 d2s(1) = -F73*d1s(1)*rrho 125 d2s(2) = -F43*d1s(2)*rrho 126 d2s(3) = -0.5d0*d1s(2)/gamma 127 gpp=gssb2(s) 128 d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1) 129 d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2) 130 d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2) 131 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 132 & +(rhom23*g 133 & + 2.d0*rho13*d1g(1) 134 & + rho43*d2g(1))*fac*2d0 135 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 136 & +(rho13*d1g(2) 137 & + rho43*d2g(2))*fac*4d0 138 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 139 & + rho43*d2g(3)*fac*8d0 140#endif 141 10 continue 142c 143 else 144c 145c ======> SPIN-UNRESTRICTED <====== 146c 147#ifdef IFCV81 148CDEC$ NOSWP 149#endif 150 do 20 n = 1, nq 151 if (rho(n,1).lt.tol_rho) goto 20 152c 153c Alpha 154c 155 if (rho(n,2).lt.tol_rho) goto 25 156 rho43 = C*(2d0*rho(n,2))**F43 157 rrho = 0.5d0/rho(n,2) 158 rho13 = F43*rho43*rrho 159#ifdef SECOND_DERIV 160 rhom23 = F13*rho13*rrho 161#endif 162 gamma = delrho(n,1,1)*delrho(n,1,1) + 163 & delrho(n,2,1)*delrho(n,2,1) + 164 & delrho(n,3,1)*delrho(n,3,1) 165 gam12 = 2d0*dsqrt(gamma) 166 if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 25 167c 168 s = Cs*gam12/rho43 169 d1s(1) = -F43*s*rrho 170 d1s(2) = 0.5d0*s/gamma 171c 172c Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1 173c 174 175 g=gssb0(s) 176 gp=gssb1(s) 177c 178 d1g(1) = gp*d1s(1) 179 d1g(2) = gp*d1s(2) 180 Ex = Ex + rho43*g*qwght(n)*fac*0.5d0 181 if(ldew)func(n) = func(n) + rho43*g*fac*0.5d0 182 Amat(n,1) = Amat(n,1) + (rho13*g+rho43*d1g(1))*fac 183 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 0.5d0*rho43*d1g(2)*fac 184#ifdef SECOND_DERIV 185 d2s(1) = -F73*d1s(1)*rrho 186 d2s(2) = -F43*d1s(2)*rrho 187 d2s(3) = -0.5d0*d1s(2)/gamma 188 gpp=gssb2(s) 189 d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1) 190 d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2) 191 d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2) 192 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 193 & +(rhom23*g 194 & + 2.d0*rho13*d1g(1) 195 & + rho43*d2g(1))*fac*2d0 196 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 197 & +(rho13*d1g(2) 198 & + rho43*d2g(2))*fac 199 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 200 & + rho43*d2g(3)*fac*0.5d0 201#endif 202c 203c Beta 204c 205 25 continue 206 if (rho(n,3).lt.tol_rho) goto 20 207 rho43 = C*(2d0*rho(n,3))**F43 208 rrho = 0.5d0/rho(n,3) 209 rho13 = F43*rho43*rrho 210#ifdef SECOND_DERIV 211 rhom23 = F13*rho13*rrho 212#endif 213 gamma = delrho(n,1,2)*delrho(n,1,2) + 214 & delrho(n,2,2)*delrho(n,2,2) + 215 & delrho(n,3,2)*delrho(n,3,2) 216 gam12 = 2d0*dsqrt(gamma) 217 if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 20 218c 219 s = Cs*gam12/rho43 220 d1s(1) = -F43*s*rrho 221 d1s(2) = 0.5d0*s/gamma 222c 223c Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1 224c 225 g=gssb0(s) 226 gp=gssb1(s) 227c 228 d1g(1) = gp*d1s(1) 229 d1g(2) = gp*d1s(2) 230 Ex = Ex + rho43*g*qwght(n)*fac*0.5d0 231 if(ldew)func(n) = func(n) + rho43*g*fac*0.5d0 232 Amat(n,2) = Amat(n,2) + (rho13*g+rho43*d1g(1))*fac 233 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + 0.5d0*rho43*d1g(2)*fac 234#ifdef SECOND_DERIV 235 d2s(1) = -F73*d1s(1)*rrho 236 d2s(2) = -F43*d1s(2)*rrho 237 d2s(3) = -0.5d0*d1s(2)/gamma 238 gpp=gssb2(s) 239 d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1) 240 d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2) 241 d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2) 242 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 243 & +(rhom23*g 244 & + 2.d0*rho13*d1g(1) 245 & + rho43*d2g(1))*fac*2d0 246 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) 247 & +(rho13*d1g(2) 248 & + rho43*d2g(2))*fac 249 Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) 250 & + rho43*d2g(3)*fac*0.5d0 251#endif 252c 253 20 continue 254 endif 255c 256 return 257 end 258#ifndef SECOND_DERIV 259#define SECOND_DERIV 260c 261c Compile source again for the 2nd derivative case 262c 263#include "xc_ssbD_1.F" 264#endif 265