1!--------------------------------------------------------------------------------------------------! 2! DFTB+: general package for performing fast atomistic simulations ! 3! Copyright (C) 2006 - 2019 DFTB+ developers group ! 4! ! 5! See the LICENSE file for terms of usage and distribution. ! 6!--------------------------------------------------------------------------------------------------! 7 8!> Contains functions to calculate the short-ranged part of gamma, and the distance beyond which it 9!> becomes negligible. 10module dftbp_shortgamma 11 use dftbp_accuracy 12 use dftbp_message 13 implicit none 14 15 private 16 17 public :: expGamma, expGammaDamped, expGammaPrime, expGammaDampedPrime 18 public :: expGammaCutoff 19 20 21 !> Used to return runtime diagnostics 22 character(len=100) :: error_string 23 24contains 25 26 27 !> Determines the cut off where the short range part goes to zero (a small constant) 28 function expGammaCutoff(U2, U1, minValue) 29 30 !> Hubbard U value in a.u. 31 real(dp), intent(in) :: U2 32 33 !> Hubbard U value in a.u. 34 real(dp), intent(in) :: U1 35 36 !> value below which the short range contribution is considered negligible. If not set this 37 !> comes from a constant in the precision module. 38 real(dp), intent(in), optional :: minValue 39 40 !> Returned cut off distance 41 real(dp) :: expGammaCutoff 42 43 real(dp) :: cutValue 44 real(dp) :: rab, MaxGamma, MinGamma, lowerGamma, gamma 45 real(dp) :: cut, MaxCutOff, MinCutOff 46 47 expGammaCutoff = 0.0_dp 48 49 if (present(minValue)) then 50 cutValue = minValue 51 else 52 cutValue = minShortGamma 53 end if 54 55 if (cutValue < tolShortGamma) then 5699000 format ('Failure in determining short-range cut-off,', & 57 & ' -ve cutoff negative :',f12.6) 58 write(error_string, 99000) cutValue 59 call error(error_string) 60 else if (U1 < minHubTol .or. U2 < minHubTol) then 6199010 format ('Failure in short-range gamma, U too small :',f12.6,f12.6) 62 write(error_string, 99010) U1, U2 63 call error(error_string) 64 end if 65 66 rab = 1.0_dp 67 do while(expGamma(rab,U2,U1) > cutValue) 68 rab = 2.0_dp*rab 69 end do 70 if (rab < 2.0_dp) then 7199020 format ('Failure in short-range gamma cut-off : ', & 72 & 'requested tolerance too large : ',f10.6) 73 write(error_string, 99020) cutValue 74 call error(error_string) 75 end if 76 ! bisection search for the value where the contribution drops below cutValue 77 MinCutOff = rab + 0.1_dp 78 MaxCutOff = 0.5_dp * rab - 0.1_dp 79 maxGamma = expGamma(MaxCutOff,U2,U1) 80 minGamma = expGamma(MinCutOff,U2,U1) 81 lowerGamma = expGamma(MinCutOff,U2,U1) 82 cut = MaxCutOff + 0.1_dp 83 gamma = expGamma(cut,U2,U1) 84 do While ((gamma-lowerGamma) > tolShortGamma) 85 MaxCutOff = 0.5_dp*(cut + MinCutOff) 86 if ((maxGamma >= minGamma) .eqv. (cutValue >= & 87 & expGamma(MaxCutOff,U2,U1))) then 88 MinCutOff = MaxCutOff 89 lowerGamma = expGamma(MinCutOff,U2,U1) 90 else 91 cut = MaxCutOff 92 gamma = expGamma(cut,U2,U1) 93 end if 94 end do 95 expGammaCutoff = MinCutOff 96 97 end function expGammaCutoff 98 99 100 !> Determines the value of the short range contribution to gamma with the exponential form 101 function expGamma(rab,Ua,Ub) 102 103 !> separation of sites a and b 104 real(dp), intent(in) :: rab 105 106 !> Hubbard U for site a 107 real(dp), intent(in) :: Ua 108 109 !> Hubbard U for site b 110 real(dp), intent(in) :: Ub 111 112 !> contribution 113 real(dp) :: expGamma 114 115 real(dp) :: tauA, tauB, tauMean 116 117 if (rab < 0.0_dp) then 11899030 format ('Failure in short-range gamma, r_ab negative :',f12.6) 119 write(error_string, 99030) rab 120 call error(error_string) 121 else if (Ua < MinHubTol) then 12299040 format ('Failure in short-range gamma, U too small :',f12.6) 123 write(error_string, 99040) Ua 124 call error(error_string) 125 else if (Ub < MinHubTol) then 12699050 format ('Failure in short-range gamma, U too small : ',f12.6) 127 write(error_string, 99050) Ub 128 call error(error_string) 129 end if 130 131 ! 16/5 * U, see review papers / theses 132 tauA = 3.2_dp*Ua 133 tauB = 3.2_dp*Ub 134 if (rab < tolSameDist) then 135 ! on-site case with R~0 136 if (abs(Ua - Ub) < MinHubDiff) then 137 ! same Hubbard U values, onsite , NOTE SIGN CHANGE! 138 expGamma = -0.5_dp*(Ua + Ub) 139 else 140 ! Ua /= Ub Hubbard U values - limiting case, NOTE SIGN CHANGE! 141 expGamma = & 142 & -0.5_dp*((tauA*tauB)/(tauA+tauB) + (tauA*tauB)**2/(tauA+tauB)**3) 143 end if 144 else if (abs(Ua - Ub) < MinHubDiff) then 145 ! R > 0 and same Hubbard U values 146 tauMean = 0.5_dp*(tauA + tauB) 147 expGamma = & 148 & exp(-tauMean*rab) * (1.0_dp/rab + 0.6875_dp*tauMean & 149 & + 0.1875_dp*rab*(tauMean**2) & 150 & + 0.02083333333333333333_dp*(rab**2)*(tauMean**3)) 151 else 152 ! using the sign convention in the review articles, not Joachim Elstner's thesis -- there's a 153 ! typo there 154 expGamma = gammaSubExprn(rab,tauA,tauB) + gammaSubExprn(rab,tauB,tauA) 155 end if 156 157 end function expGamma 158 159 160 !> Determines the value of the derivative of the short range contribution to gamma with the 161 !> exponential form 162 function expGammaPrime(rab,Ua,Ub) 163 164 !> separation of sites a and b 165 real(dp), intent(in) :: rab 166 167 !> Hubbard U for site a 168 real(dp), intent(in) :: Ua 169 170 !> Hubbard U for site b 171 real(dp), intent(in) :: Ub 172 173 !> returned contribution 174 real(dp) :: expGammaPrime 175 176 real(dp) :: tauA, tauB, tauMean 177 178 if (rab < 0.0_dp) then 17999060 format ('Failure in short-range gamma, r_ab negative :',f12.6) 180 write(error_string, 99060) rab 181 call error(error_string) 182 else if (Ua < MinHubTol) then 18399070 format ('Failure in short-range gamma, U too small :',f12.6) 184 write(error_string, 99070) Ua 185 call error(error_string) 186 else if (Ub < MinHubTol) then 18799080 format ('Failure in short-range gamma, U too small : ',f12.6) 188 write(error_string, 99080) Ub 189 call error(error_string) 190 end if 191 192 ! on-site case with R~0 193 if (rab < tolSameDist) then 194 expGammaPrime = 0.0_dp 195 else if (abs(Ua - Ub) < MinHubDiff) then 196 ! R > 0 and same Hubbard U values 197 ! 16/5 * U, see review papers 198 tauMean = 3.2_dp * 0.5_dp * (Ua + Ub) 199 expGammaPrime = & 200 & -tauMean * exp(-tauMean*rab) * & 201 & ( 1.0_dp/rab + 0.6875_dp*tauMean & 202 & + 0.1875_dp*rab*(tauMean**2) & 203 & + 0.02083333333333333333_dp*(rab**2)*(tauMean**3) ) & 204 & + exp(-tauMean*rab) * & 205 & ( -1.0_dp/rab**2 + 0.1875_dp*(tauMean**2) & 206 & + 2.0_dp*0.02083333333333333333_dp*rab*(tauMean**3) ) 207 else 208 ! 16/5 * U, see review papers 209 tauA = 3.2_dp*Ua 210 tauB = 3.2_dp*Ub 211 ! using the sign convention in the review articles, not Joachim Elstner's thesis -- there's a 212 ! typo there 213 expGammaPrime = gammaSubExprnPrime(rab,tauA,tauB) + & 214 & gammaSubExprnPrime(rab,tauB,tauA) 215 end if 216 end function expGammaPrime 217 218 219 !> Determines the value of the short range contribution to gamma with the exponential form with 220 !> damping. 221 !> See J. Phys. Chem. A, 111, 10865 (2007). 222 function expGammaDamped(rab, Ua, Ub, dampExp) 223 224 !> separation of sites a and b 225 real(dp), intent(in) :: rab 226 227 !> Hubbard U for site a 228 real(dp), intent(in) :: Ua 229 230 !> Hubbard U for site b 231 real(dp), intent(in) :: Ub 232 233 !> Damping exponent 234 real(dp), intent(in) :: dampExp 235 236 !> returned contribution 237 real(dp) :: expGammaDamped 238 239 real(dp) :: rTmp 240 241 rTmp = -1.0_dp * (0.5_dp * (Ua + Ub))**dampExp 242 expGammaDamped = expGamma(rab, Ua, Ub) * exp(rTmp * rab**2) 243 244 end function expGammaDamped 245 246 247 !> Determines the value of the derivative of the short range contribution to gamma with the 248 !> exponential form with damping 249 function expGammaDampedPrime(rab, Ua, Ub, dampExp) 250 251 !> separation of sites a and b 252 real(dp), intent(in) :: rab 253 254 !> Hubbard U for site a 255 real(dp), intent(in) :: Ua 256 257 !> Hubbard U for site b 258 real(dp), intent(in) :: Ub 259 260 !> Damping exponent 261 real(dp), intent(in) :: dampExp 262 263 !> returned contribution 264 real(dp) :: expGammaDampedPrime 265 266 real(dp) :: rTmp 267 268 rTmp = -1.0_dp * (0.5_dp *(Ua + Ub))**dampExp 269 expGammaDampedPrime = expGammaPrime(rab, Ua, Ub) * exp(rTmp * rab**2) & 270 &+ 2.0_dp * expGamma(rab, Ua, Ub) * exp(rTmp * rab**2) * rab * rTmp 271 272 end function expGammaDampedPrime 273 274 275 !> Determines the value of the short range contribution to gamma using the old Ohno/Klopman form 276 !> Caveat: This is too long ranged to use in a periodic calculation 277 function OhnoKlopman(rab,Ua,Ub) 278 279 !> separation of sites a and b 280 real(dp), intent(in) :: rab 281 282 !> Hubbard U for site a 283 real(dp), intent(in) :: Ua 284 285 !> Hubbard U for site b 286 real(dp), intent(in) :: Ub 287 288 !> contribution 289 real(dp) :: OhnoKlopman 290 291 if (Ua < MinHubTol) then 29299090 format ('Failure in short-range gamma, U too small : ',f12.6) 293 write(error_string, 99090) Ua 294 call error(error_string) 295 else if (Ub < MinHubTol) then 29699100 format ('Failure in short-range gamma, U too small : ',f12.6) 297 write(error_string, 99100) Ub 298 call error(error_string) 299 end if 300 301 OhnoKlopman = 1.0_dp/sqrt(rab**2 + 0.25_dp*(1.0_dp/Ua + 1.0_dp/Ub)**2) 302 303 end function OhnoKlopman 304 305 306 !> Determines the value of the derivative of the short range contribution to gamma using the old 307 !> Ohno/Klopman form. Caveat: This is too long ranged to use in a periodic calculation 308 function OhnoKlopmanPrime(rab,Ua,Ub) 309 310 !> separation of sites a and b 311 real(dp), intent(in) :: rab 312 313 !> Hubbard U for site a 314 real(dp), intent(in) :: Ua 315 316 !> Hubbard U for site b 317 real(dp), intent(in) :: Ub 318 319 !> contribution 320 real(dp) :: OhnoKlopmanPrime 321 322 if (Ua < MinHubTol) then 32399110 format ('Failure in short-range gamma, U too small : ',f12.6) 324 write(error_string, 99110) Ua 325 call error(error_string) 326 else if (Ub < MinHubTol) then 32799120 format ('Failure in short-range gamma, U too small : ',f12.6) 328 write(error_string, 99120) Ub 329 call error(error_string) 330 end if 331 332 OhnoKlopmanPrime = -rab / & 333 & (sqrt(rab**2 + 0.25_dp*(1.0_dp/Ua + 1.0_dp/Ub)**2)**3) 334 335 end function OhnoKlopmanPrime 336 337 338 !> Determines the value of a part of the short range contribution to the exponential gamma, when 339 !> Ua /= Ub and R > 0 340 function gammaSubExprn(rab,tau1,tau2) 341 342 !> separation of sites a and b 343 real(dp), intent(in) :: rab 344 345 !> Charge fluctuation for site a 346 real(dp), intent(in) :: tau1 347 348 !> Charge fluctuation U for site b 349 real(dp), intent(in) :: tau2 350 351 !> contribution 352 real(dp) :: gammaSubExprn 353 354 if (abs(tau1 - tau2) < 3.2_dp*MinHubDiff) then 35599130 format ('Failure in gammaSubExprn, both tau degenerate ',f12.6,f12.6) 356 write(error_string, 99130) tau1,tau2 357 call error(error_string) 358 else if (rab < tolSameDist) then 35999140 format ('Atoms on top of each other in gammaSubExprn') 360 write(error_string, 99140) 361 call error(error_string) 362 end if 363 364 gammaSubExprn = exp(-tau1 * rab) * & 365 & ( (0.5_dp*tau2**4*tau1/(tau1**2-tau2**2)**2) - & 366 & (tau2**6-3.0_dp*tau2**4*tau1**2)/(rab*(tau1**2-tau2**2)**3) ) 367 368 end function gammaSubExprn 369 370 371 !> Determines the derivative of the value of a part of the short range contribution to the 372 !> exponential gamma, when Ua /= Ub and R > 0 373 function gammaSubExprnPrime(rab,tau1,tau2) 374 375 !> separation of sites a and b 376 real(dp), intent(in) :: rab 377 378 !> Charge fluctuation for site a 379 real(dp), intent(in) :: tau1 380 381 !> Charge fluctuation U for site b 382 real(dp), intent(in) :: tau2 383 384 !> contribution 385 real(dp) :: gammaSubExprnPrime 386 387 if (abs(tau1 - tau2) < 3.2_dp*MinHubDiff) then 38899150 format ('Failure in gammaSubExprn, both tau degenerate ',f12.6,f12.6) 389 write(error_string, 99150) tau1,tau2 390 call error(error_string) 391 else if (rab < tolSameDist) then 39299160 format ('Atoms on top of each other in gammaSubExprn') 393 write(error_string, 99160) 394 call error(error_string) 395 end if 396 397 gammaSubExprnPrime = -tau1 * exp(- tau1 * rab) * & 398 &( (0.5_dp*tau2**4*tau1/(tau1**2-tau2**2)**2) - & 399 & (tau2**6-3.0_dp*tau2**4*tau1**2)/(rab*(tau1**2-tau2**2)**3) ) + & 400 & exp(- tau1 * rab) * (tau2**6-3.0_dp*tau2**4*tau1**2) & 401 & / (rab**2 *(tau1**2-tau2**2)**3) 402 403 end function gammaSubExprnPrime 404 405end module dftbp_shortgamma 406