1#ifndef SECOND_DERIV 2 Subroutine xc_op(tol_rho, whichf, 3 & fac, lfac, nlfac, rho, delrho, 4 & Amat, Cmat, nq, ipol, Ec, qwght,ldew,func) 5#else 6 Subroutine xc_op_d2(tol_rho, whichf, 7 & fac, lfac, nlfac, rho, delrho, 8 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ec, 9 & qwght,ldew,func) 10#endif 11c 12C$Id$ 13c 14 implicit none 15c 16#include "dft2drv.fh" 17c 18 double precision tol_rho, fac, Ec 19 character*4 whichf 20 integer nq, ipol 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 QABOP,QAB88OP,QABPBOP 48 Parameter (QAB88OP=2.3670D0,QABPBOP=2.3789D0) 49c 50c References: 51c Tsuneda, Suzumura, Hirao, JCP 110, 10664 (1999) 52c Tsuneda, Suzumura, Hirao, JCP 111, 5656 (1999) 53c 54c*************************************************************************** 55c 56 integer n 57 double precision rho13, rho43, gamma, x 58 double precision kalpha,kbeta, rho13a, rho13b,rhoa,rhob 59 double precision banb, hbab, hbabx 60 double precision dhdab,dhdabx,dkadra,dkbdrb,dkadxa,dkbdxb, 61 A dbabdra,dbabdrb,dbabdga,dbabdgb,dkadga,dkbdgb, 62 A dbabdka,dbabdkb 63c 64 hbabx(x) = (1.5214d0*x + 0.5764d0)/ 65 / (x**2*(x**2+1.1284d0*x+0.3183d0)) 66 dhdabx(x) = -(4.5642d0*x**4+5.7391d0*x**3+ 67 + 2.4355*x**2+0.3669d0*x)/ 68 / ((x**4+1.1284d0*x**3+0.3183d0*x**2)**2) 69c 70 if(whichf.eq.'be88') then 71 QABOP=QAB88OP 72 endif 73 if(whichf.eq.'pb96') then 74 QABOP=QABPBOP 75 endif 76 if (ipol.eq.1) then 77c 78c ======> SPIN-RESTRICTED <====== 79c 80 do 10 n = 1, nq 81 if (rho(n,1).lt.tol_rho) goto 10 82c 83c Spin alpha: 84c 85 rhoa=rho(n,1)*0.5d0 86 rho13a = (rhoa)**(1.d0/3.d0) 87 rho43 = rho13a**4 88 gamma = delrho(n,1,1)*delrho(n,1,1) + 89 & delrho(n,2,1)*delrho(n,2,1) + 90 & delrho(n,3,1)*delrho(n,3,1) 91 gamma = 0.25d0 * gamma 92 if (dsqrt(gamma).gt.tol_rho)then 93 x = dsqrt(gamma) / rho43 94 call xc_kop(tol_rho,whichf,x, 95 & kalpha, dkadxa) 96 dkadra = -(4d0/3d0)*x*dkadxa/rhoa 97 dkadga = (dkadxa/rho43)*0.5d0/dsqrt(gamma) 98 else 99 x=0d0 100 call xc_kop(tol_rho,whichf,x, 101 & kalpha, dkadxa) 102 dkadra = 0d0 103 dkadga = 0d0 104 endif 105c 106c 107 108 banb = qabop * rho13a * kalpha *0.5d0 109 110 if(banb.ne.0) then 111 dbabdra = banb*0.5d0* 112 / (1d0/(3d0*rhoa)+dkadra/kalpha) 113 114 dbabdga = banb/kalpha*dkadga*0.5d0 115 116 hbab = hbabx(banb) 117 dhdab = dhdabx(banb) 118 else 119 dbabdra =0d0 120 dbabdga =0d0 121 hbab = 0d0 122 dhdab = 0d0 123 endif 124 125 Ec = Ec - rhoa**2*hbab*qwght(n)*fac 126 if(ldew)func(n) = func(n) - rhoa**2*hbab*fac 127 Amat(n,1) = Amat(n,1) - 128 - (rhoa*hbab + rhoa**2*dhdab*dbabdra)*fac 129 130c 131 if (x.gt.tol_rho) then 132 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) - 133 - rhoa**2*dhdab*dbabdga*fac 134 endif 135c 136 10 continue 137c 138 else 139c 140c ======> SPIN-UNRESTRICTED <====== 141c 142 do 20 n = 1, nq 143 if (rho(n,1).lt.tol_rho) goto 20 144 if (rho(n,2).ge.tol_rho*0.5d0) then 145c 146c Spin alpha: 147c 148 rhoa=rho(n,2) 149 rho13a = rhoa**(1.d0/3.d0) 150 rho43 = rho13a*rhoa 151 gamma = delrho(n,1,1)*delrho(n,1,1) + 152 & delrho(n,2,1)*delrho(n,2,1) + 153 & delrho(n,3,1)*delrho(n,3,1) 154 if (dsqrt(gamma).gt.tol_rho)then 155 x = dsqrt(gamma) / rho43 156 call xc_kop(tol_rho,whichf,x, 157 & kalpha, dkadxa) 158 159 dkadra = -(4d0/3d0)*x*dkadxa/rhoa 160 dkadga = dkadxa*0.5d0/(rho43*dsqrt(gamma)) 161 else 162 x = 0d0 163 endif 164 else 165 rhoa=0d0 166 rho13a=0d0 167 x = 0d0 168 endif 169 if(x.eq.0d0) then 170 call xc_kop(tol_rho,whichf,x, 171 & kalpha, dkadxa) 172 dkadra = 0d0 173 dkadga = 0d0 174 endif 175c 176c Spin beta: 177c 178 if (rho(n,3).ge.tol_rho*0.5d0) then 179c 180 rhob=rho(n,3) 181 rho13b = rhob**(1.d0/3.d0) 182 rho43 = rho13b*rhob 183 gamma = delrho(n,1,2)*delrho(n,1,2) + 184 & delrho(n,2,2)*delrho(n,2,2) + 185 & delrho(n,3,2)*delrho(n,3,2) 186 if (dsqrt(gamma).gt.tol_rho)then 187 x = dsqrt(gamma) / rho43 188 call xc_kop(tol_rho,whichf,x, 189 & kbeta, dkbdxb) 190 191 dkbdrb = -(4d0/3d0)*x*dkbdxb/rhob 192 dkbdgb = dkbdxb*0.5d0/(rho43*dsqrt(gamma)) 193 else 194 x = 0d0 195 endif 196 else 197 if(rho13a.eq.0) goto 20 198 rhob=0d0 199 rho13b=0d0 200 x=0d0 201 endif 202 if(x.eq.0d0) then 203 call xc_kop(tol_rho,whichf,x, 204 & kbeta, dkbdxb) 205 dkbdrb = 0d0 206 dkbdgb= 0d0 207 endif 208 209 banb = qabop*(rho13a*kalpha*rho13b*kbeta)/ 210 / (rho13a*kalpha+rho13b*kbeta) 211 212 if(banb.ne.0) then 213 dbabdra = banb*kbeta*rho13b/ 214 / (rho13a*kalpha+rho13b*kbeta)* 215 / (1d0/(3d0*rhoa)+dkadra/kalpha) 216 dbabdrb = banb*kalpha*rho13a/ 217 / (rho13a*kalpha+rho13b*kbeta)* 218 / (1d0/(3d0*rhob)+dkbdrb/kbeta) 219 220 dbabdga = banb*rho13b*kbeta/ 221 / ((rho13a*kalpha+rho13b*kbeta)*kalpha)* 222 * dkadga 223 dbabdgb = banb*rho13a*kalpha/ 224 / ((rho13a*kalpha+rho13b*kbeta)*kbeta)* 225 * dkbdgb 226 227 hbab = hbabx(banb) 228 dhdab = dhdabx(banb) 229 else 230 dbabdra =0d0 231 dbabdrb =0d0 232 dbabdga =0d0 233 dbabdgb =0d0 234 hbab = 0d0 235 dhdab = 0d0 236 endif 237 238 Ec = Ec - rhoa*rhob*hbab*qwght(n)*fac 239 if (ldew) func(n) = func(n) - rhoa*rhob*hbab*fac 240 Amat(n,1) = Amat(n,1) - 241 - (rhob*hbab + rhoa*rhob*dhdab*dbabdra)*fac 242 Amat(n,2) = Amat(n,2) - 243 - (rhoa*hbab + rhoa*rhob*dhdab*dbabdrb)*fac 244c 245c 246 if (x.gt.tol_rho) then 247 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) - 248 - rhoa*rhob*dhdab*dbabdga*fac 249 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) - 250 - rhoa*rhob*dhdab*dbabdgb*fac 251 endif 252 253c 254c 255 20 continue 256c 257 endif 258c 259 return 260 end 261#ifndef SECOND_DERIV 262#define SECOND_DERIV 263c 264c Compile source again for the 2nd derivative case 265c 266#include "xc_op.F" 267#endif 268