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