1c 2c ----------------------------------------------------------------------- 3c Uniform electron gas exchange functional for the erfc(r)/r interaction 4c as implemented in the following paper: 5c "A well-tempered density functional theory of electrons in molecules" 6c Ester Livshits & Roi Baer, Phys. Chem. Chem. Phys., 9, 2932 (2007) 7c The other relevant publication is: 8c R. Baer, D. Neuhauser, Phys. Rev. Lett., 94, 043002 (2005) 9c ----------------------------------------------------------------------- 10c 11 subroutine bnl2007_x(n2ft3d,ispin,dn,x_parameter,gamma,xcp,xce) 12 implicit none 13 integer n2ft3d 14 integer ispin 15 real*8 dn(n2ft3d,2) 16 real*8 x_parameter,gamma 17 18 real*8 xcp(n2ft3d,2) 19 real*8 xce(n2ft3d) 20 21c 22c **** local variables **** 23 double precision dncut 24 parameter (dncut=1.0d-30) 25 26 integer n 27 double precision rhoA1,rhoB1,rhoTotal 28 double precision fA, fB, fpA, fpB, fppA, fppB 29 30 double precision TEpsX 31 double precision TEpsXprime 32 external TEpsX 33 external TEpsXprime 34 35c ----------------------------------------------------------------------- 36c Calculate the first derivatives 37c ----------------------------------------------------------------------- 38c 39c **** spin-restricted **** 40 if (ispin.eq.1) then 41 do n = 1,n2ft3d 42 rhoA1 = dn(n,1) + dn(n,ispin) 43 if (rhoA1.gt.dncut) then 44 fA = TEpsX(rhoA1,gamma) 45 fpA = TEpsXprime(rhoA1,gamma) 46 xcp(n,1) = xcp(n,1) + x_parameter*(fpA*rhoA1+fA) 47 xce(n) = xce(n) + x_parameter*fA 48 end if 49 end do 50 51c **** spin-unrestricted **** 52 else 53 do n = 1,n2ft3d 54 rhoTotal = dn(n,1) + dn(n,ispin) 55 if (rhoTotal.gt.dncut) then 56 rhoA1 = dn(n,1) + dn(n,1) 57 rhoB1 = dn(n,ispin) + dn(n,ispin) 58 fA = TEpsX(rhoA1,gamma) 59 fB = TEpsX(rhoB1,gamma) 60 fpA = TEpsXprime(rhoA1,gamma) 61 fpB = TEpsXprime(rhoB1,gamma) 62 xcp(n,1) = xcp(n,1) + x_parameter*(fpA*rhoA1+fA) 63 xcp(n,2) = xcp(n,2) + x_parameter*(fpB*rhoB1+fB) 64 xce(n) = xce(n) 65 > + x_parameter*( fA*dn(n,1) 66 > + fB*dn(n,ispin) )/(rhoTotal) 67 end if 68 end do 69 end if 70 71 return 72 end 73 74c 75c --------------------------------------------------------------------------------------- 76c Evaluates the actual function 77c --------------------------------------------------------------------------------------- 78c 79 double precision function THqBNL(q) 80 implicit none 81 double precision q 82 83 double precision Pi 84 parameter (Pi = 3.141592653589793d0) 85 double precision TwoSqrtPi,OneOverQ,q2 86 double precision util_erf 87 external util_erf 88 89 OneOverQ = 1.0d0/q 90 TwoSqrtPi = 2.0d0*dsqrt(Pi) 91 q2 = q**2.0d0 92 93 if (q .lt. 1.0d-15) then 94 THqBNL=1.d0 95 return 96 end if 97 98 if (q .lt. 0.1d0) then 99 THqBNL=1.0d0-q*2.0d0/3.0d0*(TwoSqrtPi-q+q*(q2-2.0d0)) 100 return 101 end if 102 103 THqBNL=1.0d0-q*2.0d0/3.0d0*(TwoSqrtPi*util_erf(OneOverQ)-q+ 104 $ q*(q2-2.0d0)*(1.0d0-exp(-OneOverQ*OneOverQ))) 105 106 return 107 end 108c 109c --------------------------------------------------------------------------------------- 110c Calculate the first derivative of the function 111c --------------------------------------------------------------------------------------- 112c 113 double precision function THqBNLPrime(q) 114 implicit none 115 double precision q 116 117 double precision Pi 118 parameter (Pi = 3.141592653589793d0) 119 double precision OneOverQ,q2,q3 120 double precision util_erf 121 external util_erf 122 123 OneOverQ = 1.0d0/q 124 q2 = q**2.0d0 125 q3 = q**3.0d0 126 127 if (q .lt. 0.1d0) then 128 THqBNLPrime = -4.0d0/3.0d0*(dsqrt(Pi)+2.0d0*q3-3.0d0*q) 129 return 130 end if 131 132 THqBNLPrime = 4.0d0/3.0d0*(q*(exp(-OneOverQ*OneOverQ)*(2.0d0*q2 133 $ -1.0d0)+(3.0d0-2.0d0*q2))-dsqrt(Pi)*util_erf(OneOverQ)) 134 135 return 136 end 137 138 139c 140c --------------------------------------------------------------------------------------- 141c Calculate the local Fermi vector for the provided density 142c --------------------------------------------------------------------------------------- 143c 144 double precision function TFermiK(den) 145 implicit none 146 double precision den 147 148 double precision Pi,F13 149 parameter (Pi = 3.141592653589793d0) 150 parameter (F13 = 1.0d0/3.0d0) 151 152 TFermiK = (3.d0*Pi*Pi*den)**F13 153 return 154 end 155c 156c --------------------------------------------------------------------------------------- 157c Calculate the first derivative of the local Fermi vector (it depends on the density) 158c --------------------------------------------------------------------------------------- 159c 160 double precision function TFermiKPrime(den) 161 implicit none 162 double precision den 163 164 double precision Pi,F23 165 parameter (Pi = 3.141592653589793d0) 166 parameter (F23 = 2.0d0/3.0d0) 167 168 169 TFermiKPrime = (Pi/(3.0d0*den))**F23 170 return 171 end 172 173c 174c --------------------------------------------------------------------------------------- 175c Calculate the first derivative of q (q=gamma/kf) (it implicitly depends on the density) 176c --------------------------------------------------------------------------------------- 177 double precision function TQPrime(gamma,kF) 178 implicit none 179 double precision kF, gamma 180 181 182 TQPrime = -gamma/(kF*kF) 183 return 184 end 185 186 187c 188c --------------------------------------------------------------------------------------- 189c Calculate the function EpsX at the given density value and gamma 190c --------------------------------------------------------------------------------------- 191c 192 double precision function TEpsX(Rho,gamma) 193 implicit none 194 double precision Rho,gamma 195 196 double precision Pi 197 parameter (Pi = 3.141592653589793d0) 198 double precision kF,Cs 199 double precision THqBNL 200 double precision TFermiK 201 external THqBNL 202 external TFermiK 203 204 if (Rho.le.0.0d0) then 205 TEpsX = 0.0d0 206 return 207 end if 208 209 kF = TFermiK(Rho) 210 Cs = -3.0D0/(4.0d0*Pi) 211 TEpsX = Cs * kF * THqBNL(gamma/kF) 212 return 213 end 214 215 216 217c 218c --------------------------------------------------------------------------------------- 219c Calculate the first derivative of TEpsX 220c --------------------------------------------------------------------------------------- 221c 222 double precision function TEpsXprime(Rho,gamma) 223 implicit none 224 double precision Rho,gamma 225 226 double precision Pi 227 parameter (Pi = 3.141592653589793d0) 228 double precision Cs,kF,CsPrime 229 230 double precision THqBNL 231 double precision THqBNLPrime 232 double precision TQPrime 233 double precision TFermiK 234 double precision TFermiKPrime 235 external THqBNL 236 external THqBNLPrime 237 external TQPrime 238 external TFermiK 239 external TFermiKPrime 240 241 242 kF = TFermiK(Rho) 243 CsPrime = -3.0D0/(4.0d0*Pi) 244 Cs = CsPrime*kF 245 246 if (Rho.le.0d0) then 247 TEpsXprime = 0.0d0 248 return 249 end if 250 251 TEpsXprime = TFermiKPrime(Rho)*(CsPrime*THqBNL(gamma/kF)+ 252 $ TQPrime(gamma,kF)*THqBNLPrime(gamma/kF)*Cs) 253 return 254 end 255 256c $Id$ 257