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