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#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
12      subroutine xc_bnl(tol_rho, fac, lfac, nlfac, rho, Amat, nq,
13     &                    ipol, Ex, qwght, ldew, func)
14#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
15c     For locations of 2nd derivatives of functionals in array
16#include "dft2drv.fh"
17      subroutine xc_bnl_d2(tol_rho, fac, lfac, nlfac, rho, Amat,
18     &                       Amat2, nq, ipol, Ex, qwght, ldew, func)
19#else
20#include "dft3drv.fh"
21      subroutine xc_bnl_d3(tol_rho, fac, lfac, nlfac, rho, Amat,
22     1                Amat2, Amat3, nq, ipol, Ex, qwght, ldew, func)
23#endif
24c
25      implicit none
26c
27#include "errquit.fh"
28#include "stdio.fh"
29c
30ccase...start
31#include "case.fh"
32ccase...end
33c
34      integer nq, ipol, n
35      double precision fac, Ex, total
36      logical ldew, lfac, nlfac
37      double precision func(*) ! value of the functional [output]
38      double precision tol_rho
39      double precision rho(nq,(ipol*(ipol+1))/2) ! charge density
40      double precision qwght(nq)                 ! quadrature weights
41      double precision Amat(nq,ipol)             ! partial first derivatives
42      double precision F(nq),RA(nq),RB(nq)
43      double precision rhoA, rhoB, rhoTotal, rhoA1, rhoB1
44      double precision gamma
45c      double precision fA, fB, fpA, fpB, fppA, fppB
46      double precision fA, fB, fpA, fpB
47      double precision EpsX
48      double precision EpsXprime
49c      double precision EpsTwoXprime
50c
51#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
52c      double precision Amat2(nq,*)               ! partial second derivatives
53      double precision Amat2(nq,NCOL_AMAT2)       ! partial second derivatives
54      double precision fppA, fppB
55      double precision EpsTwoXprime
56#endif
57#ifdef THIRD_DERIV
58      double precision Amat3(nq,NCOL_AMAT3)
59      double precision fpppA, fpppB
60      double precision EpsThreeXprime
61#endif
62c
63c     -----------------------------------------------------------------------
64c     Preliminaries
65c     -----------------------------------------------------------------------
66c
67      gamma = cam_omega
68c
69      do n = 1,nq
70         if (ipol.eq.1) then   ! spin-restricted
71            rA(n) = rho(n,1)
72            rB(n) = 0.d0
73         else                  ! spin-unrestricted
74            rA(n) = rho(n,2)
75            rB(n) = rho(n,3)
76         end if
77      end do
78c
79c     -----------------------------------------------------------------------
80c     Calculate the first and second derivatives
81c     -----------------------------------------------------------------------
82c
83      total = 0.d0
84      do n = 1,nq
85         rhoA = rA(n)
86         rhoB = rB(n)
87         rhoTotal  = rhoA + rhoB   ! total density at point
88         if (rhoTotal.gt.tol_rho) then
89
90            if (ipol.eq.1) then    ! spin-restricted
91              rhoA1 = rhoA
92              rhoB1 = rhoB
93            else                   ! spin-unrestricted
94              rhoA1 = rhoA*2.0d0
95              rhoB1 = rhoB*2.0d0
96            end if
97
98            fA   = EpsX(rhoA1,gamma)
99            fB   = EpsX(rhoB1,gamma)
100            fpA  = EpsXprime(rhoA1,gamma)
101            fpB  = EpsXprime(rhoB1,gamma)
102
103            f(n) = fA * rhoA + fB * rhoB
104            Amat(n,1) = Amat(n,1) + (fpA*rhoA1+fA)*fac
105            if (ipol.gt.1) then
106              Amat(n,2) = Amat(n,2) + (fpB*rhoB1+fB)*fac
107            end if
108
109#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
110            fppA = EpsTwoXprime(rhoA1,gamma)
111            fppB = EpsTwoXprime(rhoB1,gamma)
112c
113            if (ipol.eq.1) then
114             Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) +
115     &         (fppA*rhoA+2.0d0*fpA)*fac*2.0d0
116            else
117             Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) +
118     &         (fppA*rhoA+fpA)*fac*4.0d0
119c           Guard against case of no beta electrons, e.g. H atom
120             if (rho(n,3).gt.tol_rho) then
121c              Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) +
122c     &          ((fppB*rhoB+fpB)*4)*fac
123              Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) +
124     &          (fppB*rhoB+fpB)*fac*4.0d0
125             end if
126            end if
127#endif
128#ifdef THIRD_DERIV
129            fpppA = EpsThreeXprime(rhoA1,gamma)
130            fpppB = EpsThreeXprime(rhoB1,gamma)
131c
132            if (ipol.eq.1) then
133             Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA)
134     1         + ( fpppA*rhoA + 3.0d0*fppA )*fac*4.0d0
135            else
136c             Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA)
137c     1         + ( fpppA*rhoA + 3.0d0*fppA )*fac*4.0d0
138             Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA)
139     1         + ( fpppA*rhoA + 1.5d0*fppA )*fac*8.0d0
140c           Guard against case of no beta electrons, e.g. H atom
141             if (rho(n,3).gt.tol_rho) then
142c              Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB)
143c     &          + ( fpppB*rhoB + 3.0d0*fppB )*fac*4.0d0
144              Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB)
145     &          + ( fpppB*rhoB + 1.5d0*fppB )*fac*8.0d0
146             end if
147            end if
148#endif
149            if (ldew) func(n) = func(n) + f(n)*fac
150            total = total + f(n)*qwght(n)
151         end if
152      end do
153
154      Ex = Ex + total*fac
155
156      return
157      end
158c
159c
160#ifndef SECOND_DERIV
161#define SECOND_DERIV
162c
163#include "xc_bnl.F"
164c
165c     ---------------------------------------------------------------------------------------
166c     Utility functions
167c     ---------------------------------------------------------------------------------------
168c
169c     ---------------------------------------------------------------------------------------
170c     Return the value of pi
171c     ---------------------------------------------------------------------------------------
172c
173      double precision function ValueOfPi()
174c
175      implicit none
176#include "xc_params.fh"
177c
178      ValueOfPi = pi
179
180      return
181      end
182c
183c     ---------------------------------------------------------------------------------------
184c     Evaluates the actual function
185c     ---------------------------------------------------------------------------------------
186c
187      double precision function HqBNL(q)
188
189      implicit none
190#include "xc_params.fh"
191
192      double precision q,TwoSqrtPi,OneOverQ,q2,DERF
193
194      if (q .lt. 1D-15) then
195         HqBNL=1.d0
196         return
197      end if
198
199      OneOverQ = 1.0d0/q
200      TwoSqrtPi = 2.0d0*dsqrt(pi)
201      q2 = q**2.0d0
202
203      if (q .lt. 0.1d0) then
204         HqBNL=1.0d0-q*2.0d0/3.0d0*(TwoSqrtPi-q+q*(q2-2.0d0))
205         return
206      end if
207
208      HqBNL=1.0d0-q*2.0d0/3.0d0*(TwoSqrtPi*DERF(OneOverQ)-q+
209     $     q*(q2-2.0d0)*(1.0d0-exp(-OneOverQ*OneOverQ)))
210
211      return
212      end
213c
214c     ---------------------------------------------------------------------------------------
215c     Calculate the local Fermi vector for the provided density
216c     ---------------------------------------------------------------------------------------
217c
218      double precision function FermiK(den)
219
220      implicit none
221#include "xc_params.fh"
222
223      double precision F13, den
224
225      F13 = 1.0D0 / 3.0D0
226      FermiK = (3.d0*pi*pi*den)**F13
227
228      return
229      end
230c
231c     ---------------------------------------------------------------------------------------
232c     Calculate the function EpsX at the given density value and gamma
233c     ---------------------------------------------------------------------------------------
234c
235      double precision function EpsX(Rho,gamma)
236
237      implicit none
238#include "xc_params.fh"
239
240      double precision  kF,RHO,gamma,Cs
241      double precision HqBNL
242      double precision FermiK
243
244      if (RHO.le.0D0) then
245         EpsX = 0.0D0
246         return
247      end if
248
249      kF = FermiK(Rho)
250      Cs = -3.0D0/(4.0d0*pi)
251      EpsX = Cs * kF * HqBNL(gamma/kF)
252
253      return
254      end
255c
256c     ---------------------------------------------------------------------------------------
257c     Calculate the first derivative of the function
258c     ---------------------------------------------------------------------------------------
259c
260      double precision function HqBNLPrime(q)
261
262      implicit none
263#include "xc_params.fh"
264
265      double precision q,OneOverQ,q2,q3,DERF
266
267      q2 = q**2.0d0
268      q3 = q**3.0d0
269      if (q .lt. 0.1d0) then
270        HqBNLPrime = -4.0d0/3.0d0*(dsqrt(Pi)+2.0d0*q3-3.0d0*q)
271        return
272      end if
273
274      OneOverQ = 1.0d0/q
275
276      HqBNLPrime = 4.0d0/3.0d0*(q*(exp(-OneOverQ*OneOverQ)*(2.0d0*q2
277     $     -1.0d0)+(3.0d0-2.0d0*q2))-dsqrt(Pi)*DERF(OneOverQ))
278
279      return
280      end
281c
282c     ---------------------------------------------------------------------------------------
283c     Calculate the first derivative of the local Fermi vector (it depends on the density)
284c     ---------------------------------------------------------------------------------------
285c
286      double precision function FermiKPrime(den)
287
288      implicit none
289#include "xc_params.fh"
290
291      double precision F23, den
292
293      F23 = 2.0D0 / 3.0D0
294      FermiKPrime = (Pi/(3.0d0*den))**F23
295
296      return
297      end
298c
299c     ---------------------------------------------------------------------------------------
300c     Calculate the first derivative of q (q=gamma/kf) (it implicitly depends on the density)
301c     ---------------------------------------------------------------------------------------
302c
303      double precision function QPrime(gamma,kF)
304
305      implicit none
306
307      double precision  kF, FermiK2, gamma
308
309      FermiK2 = kF**2.0d0
310      QPrime = -gamma/FermiK2
311
312      return
313      end
314c
315c     ---------------------------------------------------------------------------------------
316c     Calculate the first derivative of EpsX
317c     ---------------------------------------------------------------------------------------
318c
319      double precision function EpsXprime(Rho,gamma)
320
321      implicit none
322#include "xc_params.fh"
323
324      double precision Rho,gamma
325      double precision Cs,kF,CsPrime
326
327      double precision HqBNL
328      double precision HqBNLPrime
329      double precision QPrime
330      double precision FermiK
331      double precision FermiKPrime
332
333      kF = FermiK(Rho)
334      CsPrime = -3.0D0/(4.0d0*Pi)
335      Cs = CsPrime*kF
336
337      if (Rho.le.0d0) then
338         EpsXprime = 0.0d0
339         return
340      end if
341
342      EpsXprime = FermiKPrime(Rho)*(CsPrime*HqBNL(gamma/kF)+
343     $     QPrime(gamma,kF)*HqBNLPrime(gamma/kF)*Cs)
344
345      return
346      end
347c
348c     ---------------------------------------------------------------------------------------
349c     Calculate the second derivative of the main function that consititutes the functional
350c     ---------------------------------------------------------------------------------------
351c
352      double precision function HqBNLTwoPrime(q)
353
354      implicit none
355#include "xc_params.fh"
356
357      double precision q,OneOverQ,q2
358
359      q2 = q**2.0d0
360      if (q .lt. 0.1d0) then
361         HqBNLTwoPrime = 4.0d0-8.0d0*q2
362         return
363      end if
364
365      OneOverQ = 1.0d0/q
366
367      HqBNLTwoPrime = exp(-OneOverQ*OneOverQ)*(4.0d0+8.0d0*q2)
368     $     -8.0d0*q2+4.0d0
369
370      return
371      end
372c
373c     ---------------------------------------------------------------------------------------
374c     Calculate the second derivative of the local Fermi vector
375c     ---------------------------------------------------------------------------------------
376c
377      double precision function FermiKTwoPrime(den)
378
379      implicit none
380#include "xc_params.fh"
381
382      double precision F13, den
383
384      F13 = 1.0D0/3.0D0
385      FermiKTwoPrime =  -(8.0d0*Pi**2.0d0/(243.0d0*den**5.0d0))**F13
386
387      return
388      end
389c
390c     ---------------------------------------------------------------------------------------
391c     Calculate the second derivative of q
392c     ---------------------------------------------------------------------------------------
393c
394      double precision function QTwoPrime(gamma,kF)
395
396      implicit none
397
398      double precision gamma, kF, FermiK3
399
400      FermiK3 = kF**3.0d0
401      QTwoPrime = (2.0d0*gamma)/FermiK3
402
403      return
404      end
405c
406c     ---------------------------------------------------------------------------------------
407c     Calculate the second derivative of EpsX
408c     ---------------------------------------------------------------------------------------
409c
410      double precision function EpsTwoXprime(Rho,gamma)
411
412      implicit none
413#include "xc_params.fh"
414
415      double precision Rho,gamma
416      double precision kF,kFPrim,kFprim2,kF2prim
417      double precision q,qprim,qprim2,q2prim
418      double precision g,gprim,g2prim
419      double precision Cs,CsPrim
420
421      double precision FermiK
422      double precision FermiKPrime
423      double precision FermiKTwoPrime
424      double precision QPrime
425      double precision QTwoPrime
426      double precision HqBNL
427      double precision HqBNLPrime
428      double precision HqBNLTwoPrime
429
430      if (Rho.le.0d0) then
431         EpsTwoXprime = 0.0d0
432         return
433      end if
434
435      kF = FermiK(Rho)
436      kFPrim = FermiKPrime(Rho)
437      kFPrim2=kFPrim**2.0d0
438      kF2prim = FermiKTwoPrime(Rho)
439      CsPrim = -3.0d0/(4.0d0*Pi)
440      Cs = CsPrim * kF
441      q = gamma / kF
442      qprim = QPrime(gamma,kF)
443      Qprim2=qprim**2.0d0
444      q2prim = QTwoPrime(gamma,kF)
445      g = HqBNL(q)
446      gprim = HqBNLPrime(q)
447      g2prim = HqBNLTwoPrime(q)
448
449      EpsTwoXprime =
450     $     kFPrim2*(2.0d0*CsPrim*gprim*qprim
451     $     +Cs*(QPrim2*g2prim+gprim*Q2Prim))
452     $     +kF2Prim*(g*CsPrim+Cs*gprim*qprim)
453
454      return
455      end
456#endif
457#ifndef THIRD_DERIV
458#define THIRD_DERIV
459c
460#include "xc_bnl.F"
461c
462c
463c     ---------------------------------------------------------------------------------------
464c     Calculate the third derivative of the main function that consititutes the functional
465c     ---------------------------------------------------------------------------------------
466c
467      double precision function HqBNLThreePrime(q)
468
469      implicit none
470#include "xc_params.fh"
471
472      double precision q,OneOverQ
473      double precision q2, q3, q4
474
475      if (q .lt. 0.1d0) then
476         HqBNLThreePrime = -16.0d0*q
477         return
478      end if
479
480      OneOverQ = 1.0d0/q
481      q2 = q*q
482      q3 = q2*q
483      q4 = q3*q
484
485      HqBNLThreePrime = 8.0d0*( exp(-OneOverQ*OneOverQ)
486     1                        + 2.0d0*q2*exp(-OneOverQ*OneOverQ)
487     2                        - 2.0d0*q4
488     3                        + 2.0d0*q4*exp(-OneOverQ*OneOverQ) ) / q3
489
490      return
491      end
492c
493c     ---------------------------------------------------------------------------------------
494c     Calculate the third derivative of the local Fermi vector
495c     ---------------------------------------------------------------------------------------
496c
497      double precision function FermiKThreePrime(den)
498
499      implicit none
500#include "xc_params.fh"
501
502      double precision F13, den
503
504      F13 = 1.0D0/3.0D0
505      FermiKThreePrime =  (10.0d0/9.0d0)*
506     1                  (Pi**2.0d0/(9.0d0*den**8.0d0))**F13
507
508      return
509      end
510c
511c     ---------------------------------------------------------------------------------------
512c     Calculate the third derivative of q
513c     ---------------------------------------------------------------------------------------
514c
515      double precision function QThreePrime(gamma,kF)
516
517      implicit none
518
519      double precision gamma, kF, FermiK4
520
521      FermiK4 = kF**4.0d0
522      QThreePrime = -(6.0d0*gamma)/FermiK4
523
524      return
525      end
526c
527c     ---------------------------------------------------------------------------------------
528c     Calculate the third derivative of EpsX
529c     ---------------------------------------------------------------------------------------
530c
531      double precision function EpsThreeXprime(Rho,gamma)
532
533      implicit none
534#include "xc_params.fh"
535
536      double precision Rho, gamma
537c Fermi wavevector stuff
538      double precision kF, kFPrim, kF2prim, kF3prim
539      double precision kFPrim2, kFPrim3
540c q stuff
541      double precision q, qprim, q2prim, q3prim
542      double precision qprim2, qprim3
543c H(q) stuff
544      double precision g, gprim, g2prim, g3prim
545      double precision Cs
546
547      double precision FermiK
548      double precision FermiKPrime
549      double precision FermiKTwoPrime
550      double precision FermiKThreePrime
551      double precision QPrime
552      double precision QTwoPrime
553      double precision QThreePrime
554      double precision HqBNL
555      double precision HqBNLPrime
556      double precision HqBNLTwoPrime
557      double precision HqBNLThreePrime
558
559      if (Rho.le.0d0) then
560         EpsThreeXprime = 0.0d0
561         return
562      end if
563
564      kF = FermiK(Rho)
565      kFPrim = FermiKPrime(Rho)
566      kF2Prim = FermiKTwoPrime(Rho)
567      kF3Prim = FermiKThreePrime(Rho)
568c
569      kFPrim2 = kFPrim**2.0d0
570      kFPrim3 = kFPrim**3.0d0
571c
572      Cs = -3.0d0/(4.0d0*Pi)
573c
574      q = gamma / kF
575      qprim = QPrime(gamma,kF)
576      q2prim = QTwoPrime(gamma,kF)
577      q3prim = QThreePrime(gamma,kF)
578c
579      qprim2 = qprim**2.0d0
580      qprim3 = qprim**3.0d0
581c
582      g = HqBNL(q)
583      gprim = HqBNLPrime(q)
584      g2prim = HqBNLTwoPrime(q)
585      g3prim = HqBNLThreePrime(q)
586
587      EpsThreeXprime = Cs*kFPrim3*( 3.0d0*qprim2*g2prim
588     1                            + 3.0d0*kF*qprim*g2prim*q2prim
589     2                            + kF*qprim3*g3prim
590     3                            + gprim*( 3.0d0*q2prim
591     4                                    + kF*q3prim ) )
592     5               + 3.0d0*Cs*kFPrim*kF2Prim*( kF*qprim2*g2prim
593     6                                         + gprim*( 2.0d0*qprim
594     7                                                 + kF*q2prim ) )
595     8               + Cs*kF3Prim*( g + kF*gprim*qprim )
596
597      return
598      end
599#endif
600c
601c $Id$
602