1#include "dft2drv.fh" 2c Bc95 correlation functional 3c META GGA 4C utilizes ingredients: 5c rho - density 6c delrho - gradient of density 7c tau (tauN)- K.S kinetic energy density 8 9 Subroutine xc_bc95(tol_rho, cfac, lcfac, nlcfac, rho, delrho, 10 & nq, ipol, Ec, qwght, ldew, func, 11 & tau, Amat, Cmat, Mmat,ijmswitch) 12 13 14c 15c$Id$ 16c 17c Reference 18c Becke, A. D. J. Chem. Phys. 1996, 104, 1040. 19c 20 implicit none 21c 22c 23c 24c Input and other parameters 25c 26 integer ipol, nq 27 28 double precision cfac 29 logical lcfac, nlcfac 30 31 logical lfac, nlfac 32 double precision fac 33 double precision tol_rho 34c 35 logical ldew 36 double precision func(*) 37c 38c Threshold parameters 39c 40 double precision DTol,F1, F2, F3, F4,COpp 41 Data F1/1.0d0/,F2/2.0d0/, 42 & F3/3.0d0/,F4/4.0d0/ 43c 44c Correlation energy 45c 46 double precision Ec 47c 48c Charge Density 49c 50 double precision rho(nq,ipol*(ipol+1)/2) 51c 52c Charge Density Gradient 53c 54 double precision delrho(nq,3,ipol), gammaval, gam12 55 56c 57c Kinetic Energy Density 58c 59 double precision tau(nq,ipol), tauN 60 61c 62c Quadrature Weights 63c 64 double precision qwght(nq) 65c 66c Sampling Matrices for the XC Potential 67c 68 double precision Amat(nq,ipol), Cmat(nq,*) 69 double precision Mmat(nq,*) 70 71 integer n, ijmswitch 72 73c call to the bc95css subroutine 74 double precision PA,GAA,TA,FA,FPA,FGA,FTA,EUA,EUEGA,ChiA,EUPA 75 &,ChiAP,ChiAG 76 double precision PB,GBB,TB,FB,FPB,FGB,FTB,EUB,EUEGB,ChiB,EUPB 77 &,ChiBP,ChiBG 78c 79 double precision sop 80 double precision Pi, F6, F43, Pi34, F13, 81 &RS,RSP,Zeta,dZdA,dZdB,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ 82 double precision P, EUEG,Denom, DenPA, DenPB, DenGA, DenGB 83 double precision EUEGPA,EUEGPB 84 85 86c 87c ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <====== 88c 89 DTol=tol_rho 90 sop=1.0d0 91 if (ijmswitch.eq.1) then 92C Parameters for PW6B95 Correlation 93 COpp=0.00262d0 94 elseif (ijmswitch.eq.2) then 95C Parameters for PWB6K Correlation 96 COpp=0.00353d0 97 else 98C Parameters for B95 (data statement is unsafe as values may survive 99C between calls) 100 COpp=0.0031d0 101 endif 102 Pi = F4*ATan(F1) 103 F6=6.0d0 104 F43 = F4 / F3 105 Pi34 = F3 / (F4*Pi) 106 F13 = F1 / F3 107 108 do 20 n = 1, nq 109 if (rho(n,1).lt.DTol) goto 20 110 if (ipol.eq.1) then 111c 112c get the density, gradient, and tau for the alpha spin from the total 113c 114 PA = rho(n,1)/F2 115 GAA = ( delrho(n,1,1)*delrho(n,1,1) + 116 & delrho(n,2,1)*delrho(n,2,1) + 117 & delrho(n,3,1)*delrho(n,3,1))/4.0d0 118c In the bc95css subroutine, we use 2*TA as the tau, so we do not divide 119c the tau by 2 here 120 121 TA = tau(n,1) 122! TA=0.0005d0 123 124 Call bc95ss(PA,GAA,TA,FA,FPA,FGA,FTA,EUA, 125 & ChiA,EUPA,ChiAP,ChiAG,ijmswitch) 126 PB = PA 127 GBB = GAA 128 TB = TA 129 FB = FA 130 FPB = FPA 131 FGB = FGA 132 FTB = FTA 133 EUB = EUA 134 ChiB = ChiA 135 EUPB = EUPA 136 ChiBP = ChiAP 137 ChiBG = ChiAG 138 139 Ec = Ec + 2.0d0*FA*qwght(n) !factor of 2 account for both spin 140 if(ldew) func(n)=func(n)+ 2.0d0*FA 141 Amat(n,1)=Amat(n,1)+ FPA 142 Cmat(n,D1_GAA)= Cmat(n,D1_GAA) + FGA 143 Mmat(n,1)= Mmat(n,1) + FTA 144#if 0 145 write (0,'(A,3F20.6)') " Amat Cmat Mmat",FPA,FGA,FTA 146#endif 147 148 149c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted 150 else ! ipol=2 151c 152c ======> SPIN-UNRESTRICTED <====== 153c 154c 155c alpha 156c 157 158 PA = rho(n,2) 159 if (PA.le.DTol) go to 25 160 GAA = delrho(n,1,1)*delrho(n,1,1) + 161 & delrho(n,2,1)*delrho(n,2,1) + 162 & delrho(n,3,1)*delrho(n,3,1) 163c 164c In the bc95css subroutine, we use 2*TA as the tau 165c 166 TA = tau(n,1)*2.0d0 167 168 Call bc95ss(PA,GAA,TA,FA,FPA,FGA,FTA,EUA, 169 & ChiA,EUPA,ChiAP,ChiAG,ijmswitch) 170 Ec = Ec + FA*qwght(n) 171 if(ldew) func(n)=func(n)+ FA 172 Amat(n,1)=Amat(n,1)+ FPA 173 Cmat(n,D1_GAA)= Cmat(n,D1_GAA) + FGA 174c 2*0.5=1.0 for Mmat 175 Mmat(n,1)= Mmat(n,1) + FTA 176#if 0 177 write (0,'(A,3F20.6)') "AAmat Cmat Mmat",FPA,FGA,FTA 178#endif 179c 180c In the bc95css subroutine, we use 2*TA as the tau, 181c 182c 183c Beta 184c 185 25 continue 186 PB = rho(n,3) 187 if(PB.le.DTol) go to 30 188 GBB = delrho(n,1,2)*delrho(n,1,2) + 189 & delrho(n,2,2)*delrho(n,2,2) + 190 & delrho(n,3,2)*delrho(n,3,2) 191 192 TB = tau(n,2)*2.0d0 193 194 Call bc95ss(PB,GBB,TB,FB,FPB,FGB,FTB,EUB, 195 & ChiB,EUPB,ChiBP,ChiBG,ijmswitch) 196 Ec = Ec + FB*qwght(n) 197 if(ldew) func(n)=func(n)+ FB 198 Amat(n,2)= Amat(n,2)+ FPB 199 Cmat(n,D1_GBB)= Cmat(n,D1_GBB) + FGB 200 Mmat(n,2)= Mmat(n,2) + FTB 201#if 0 202 write (0,'(A,3F20.6)') "BAmat Cmat Mmat",FPB,FGB,FTB 203#endif 204 endif 205 30 continue 206 P = rho(n,1) 207 If(PA.gt.DTol.and.PB.gt.DTol) then 208 RS = (Pi34/P) ** F13 209 RSP = -RS/(F3*P) 210 Zeta = (PA-PB)/P 211 dZdA = (F1-Zeta)/P 212 dZdB = (-F1-Zeta)/P 213 Call lsdac(dtol, 214 D RS,Zeta,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ, 215 $ d2LdZZ) 216 EUEG = P*PotLC - EUA - EUB 217 Denom = F1 + COpp*(ChiA+ChiB) 218 Ec = Ec + sop*EUEG*qwght(n)/Denom 219 if(ldew) func(n)=func(n)+ sop*EUEG/Denom 220 DenPA = COpp*ChiAP 221 DenPB = COpp*ChiBP 222 DenGA = COpp*ChiAG 223 DenGB = COpp*ChiBG 224 EUEGPA = PotLC + P*dLdS*RSP + P*dLdZ*dZdA - EUPA 225 EUEGPB = PotLC + P*dLdS*RSP + P*dLdZ*dZdB - EUPB 226 if (ipol.eq.1) then 227 Amat(n,1) = Amat(n,1) + 228 & sop*(EUEGPA/Denom - EUEG*DenPA/Denom**2) 229 Cmat(n,D1_GAA)= Cmat(n,D1_GAA) - sop*(EUEG*DenGA/Denom**2) 230 else 231 Amat(n,1) = Amat(n,1) + 232 & sop*(EUEGPA/Denom - EUEG*DenPA/Denom**2) 233 Amat(n,2) = Amat(n,2) + 234 & sop*(EUEGPB/Denom - EUEG*DenPB/Denom**2) 235 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) - sop*EUEG*DenGA/Denom**2 236 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) - sop*EUEG*DenGB/Denom**2 237 endif 238 endIf 239c write (*,*) "Amat(n,1),Cmat(n,1),Mmat(n,1)",Amat(n,1),Cmat(n,1) 240c & ,Mmat(n,1) 241c stop 24220 continue 243 end 244 245 Subroutine xc_bc95_d2() 246 call errquit(' bc95: d2 not coded ',0,0) 247 return 248 end 249 250 251 252 253 Subroutine bc95ss(PX,GX,TX,F,FP,FG,FT,EUEG,Chi,EUEGP, 254 & ChiP,ChiG,ijmswitch) 255 Implicit none 256C 257C Compute the same-spin part of the bc95 correlation functional for one grid 258C point and one spin-case. 259C 260C 261 integer ijmswitch 262 double precision PX, GX, TX, F, FP, FG, FT, DTol 263 double precision EUEG, Chi, EUEGP, ChiP, ChiG,Css 264 double precision Zero, Pt25, F1, F2, F3, F4, F5, F6, F8, F11 265 double precision Pi, Pi34, F13, F23, F43, F53, F83, F113 266 double precision RS, FDUEG, D, RSP,DUEG, Denom, PotLC 267 double precision E, DenomG, DenomP, DUEGP, DP, DG, DT 268 double precision d2LdSS,d2LdSZ,d2LdZZ,dLdS,dLdZ 269 270 271 272 Data Zero/0.0d0/, Pt25/0.25d0/, F1/1.0d0/, F2/2.0d0/, F3/3.0d0/, 273 $ F4/4.0d0/, F5/5.0d0/, F6/6.0d0/, F8/8.0d0/, F11/11.0d0/ 274C 275 if (ijmswitch.eq.1) then 276C Parameters for PW6B95 Correlation 277 Css=0.03668d0 278 elseif (ijmswitch.eq.2) then 279C Parameters for PWB6K Correlation 280 Css=0.04120d0 281 else 282C Parameters for B95 (data statement is unsafe as values may survive 283C between calls) 284 Css=0.038d0 285 endif 286 DTol =1.0d-6 287 If(PX.le.DTol) then 288 EUEG = Zero 289 Chi = Zero 290 EUEGP = Zero 291 ChiP = Zero 292 ChiG = Zero 293 PX = Zero 294 GX = Zero 295 TX = Zero 296 F = Zero 297 FP = Zero 298 FG = Zero 299 FT = Zero 300 else 301 Pi = F4*ATan(F1) 302 Pi34 = F3 / (F4*Pi) 303 F13 = F1 / F3 304 F23 = F2 / F3 305 F43 = F2 * F23 306 F53 = F5 / F3 307 F83 = F8 / F3 308 F113 = F11 / F3 309 FDUEG = (F3/F5)*(F6*Pi*Pi)**F23 310 RS = (Pi34/PX) ** F13 311 Call lsdac(dtol, 312 D RS,F1,PotLC,dLdS,dLdZ,d2LdSS,d2LdSZ,d2LdZZ) 313 EUEG = PX*PotLC 314 D = TX - Pt25*GX/PX 315 DUEG = FDUEG*PX**F53 316 Chi = GX/PX**F83 317 Denom = F1 + Css*Chi 318 E = D*EUEG/(DUEG*Denom*Denom) 319c write (*,*) "ijmswitch, Css, E= ",ijmswitch, Css, E 320c stop 321 F = E 322c 323 RSP = -RS/(F3*Px) 324 ChiG = F1/PX**F83 325 ChiP = -F83*Chi/PX 326 DenomG = Css*ChiG 327 DenomP = Css*ChiP 328 DUEGP = F53*DUEG/PX 329 DP = Pt25*GX/PX**2 330 DG = -Pt25/PX 331 DT = F1 332 EUEGP = PotLC + PX*dLdS*RSP 333 FP = DP*EUEG/(DUEG*Denom*Denom) + 334 $ D*EUEGP/(DUEG*Denom*Denom) 335 $ - D*EUEG*DUEGP/(DUEG*Denom)**2 - 336 $ F2*D*EUEG*DenomP/(DUEG*Denom*Denom*Denom) 337 FG =DG*EUEG/(DUEG*Denom*Denom) - 338 $ F2*D*EUEG*DenomG/(DUEG*Denom*Denom*Denom) 339 FT =DT*EUEG/(DUEG*Denom*Denom) 340 Endif 341 Return 342 End 343 344 345