1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
2#if !defined(NWAD_PRINT)
3!>
4!> \ingroup nwxc
5!> @{
6!>
7!> \file nwxcP_xc_b97.c
8!> The code for the B97 family of functionals
9!>
10!> \brief Evaluate functionals from the B97 family of functionals
11!>
12!> These functions evaluate the exchange and correlation terms of the B97
13!> family of functionals [1]. The code was provided by Tobias Risthaus [2]
14!> and based on [3].
15!> The code was initially modified to integrate it into the NWDFT
16!> framework. Later it was modified to integrate it into the NWXC
17!> framework. These latter changes included separating out the correlation
18!> and exchange terms as well as the usual changes to the top-level
19!> interface.
20!>
21!> With the introduction of the NWAD framework for automatic differentiation
22!> the routines were translated into Fortran.
23!>
24!> Obviously we maintained the structure that allows 3 and 5 term functionals.
25!>
26!> ### References ###
27!>
28!> [1] A.D. Becke, "Density-functional thermochemistry. V. Systematic
29!>     optimization of exchange-correlation functionals", J. Chem. Phys.
30!>     107 (1997) 8554-8560, DOI:
31!>     <a href="https://doi.org/10.1063/1.475007">
32!>     10.1063/1.475007</a>.
33!>
34!> [2] Tobias Risthaus (tobias.risthaus@thch.uni-bonn.de),
35!>     Mulliken Center for Theoretical Chemistry,
36!>     Institut für Physikalische und Theoretische Chemie,
37!>     Universität Bonn, Beringstr. 4, D-53115 Bonn, Germany.
38!>
39!> [3] S. Grimme, "Semiempirical GGA-type density functional constructed
40!>     with a long-range dispersion correction", J. Comput. Chem. 27 (2006)
41!>     1787-1799, DOI:
42!>     <a href="https://doi.org/10.1002/jcc.20495">
43!>     10.1002/jcc.20495</a>.
44!> @}
45!>
46#endif
47#endif
48!>
49!> \ingroup nwxc_priv
50!> \@{
51#include "nwxcP_xc_b97.fh"
52
53#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
54#ifdef NWAD_PRINT
55      subroutine nwxcP_xc_gss3tar_p(r, g, alpha, a, f)
56#else
57      subroutine nwxcP_xc_gss3tar(r, g, alpha, a, f)
58#endif
59#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
60      subroutine nwxcP_xc_gss3tar_d2(r, g, alpha, a, f)
61#else
62      subroutine nwxcP_xc_gss3tar_d3(r, g, alpha, a, f)
63#endif
64#include "nwad.fh"
65      implicit none
66      type(nwad_dble):: r, g, f
67      double precision alpha
68#ifdef NWAD_PRINT
69#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
70      type(nwad_dble)::a(3)
71#else
72      double precision a(3)
73#endif
74#else
75      double precision a(3)
76#endif
77!     B97 like
78!     S. Grimme, J. Comput. Chem., 2006, 27, 1787-1799
79!     doi: 10.1002/jcc.20495
80!     for now, this is only valid for 3 terms.
81
82      type(nwad_dble):: r83,r113,r143,s2,s2p1,u
83
84      r83=r**(8.0d0/3.0d0)
85!     r113=r83*r
86!     r143=r113*r
87      s2=g/r83
88      s2p1=1.0d0+alpha*s2
89
90      u=alpha*s2/s2p1
91
92      f=a(1) + a(2)*u + a(3)*u**2
93
94      end
95
96#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
97#ifdef NWAD_PRINT
98      subroutine nwxcP_xc_gss5tar_p(r, g, alpha, a, f)
99#else
100      subroutine nwxcP_xc_gss5tar(r, g, alpha, a, f)
101#endif
102#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
103      subroutine nwxcP_xc_gss5tar_d2(r, g, alpha, a, f)
104#else
105      subroutine nwxcP_xc_gss5tar_d3(r, g, alpha, a, f)
106#endif
107#include "nwad.fh"
108      implicit none
109      type(nwad_dble):: r, g, f
110      double precision alpha
111#ifdef NWAD_PRINT
112#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
113      type(nwad_dble)::a(5)
114#else
115      double precision a(5)
116#endif
117#else
118      double precision a(5)
119#endif
120
121!     B97 like
122!     S. Grimme, J. Comput. Chem., 2006, 27, 1787-1799
123!     doi: 10.1002/jcc.20495
124!     for now, this is only valid for 5 terms.
125
126      type(nwad_dble):: r83,r113,r143,s2,s2p1,u
127
128      r83=r**(8.0d0/3.0d0)
129!     r113=r83*r
130!     r143=r113*r
131      s2=g/r83
132      s2p1=1.0d0+alpha*s2
133
134      u=alpha*s2/s2p1
135
136      f=a(1) + a(2)*u + a(3)*(u**2) + a(4)*(u**3) + a(5)*(u**4)
137
138      end
139
140#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
141#ifdef NWAD_PRINT
142      subroutine nwxcP_xc_gab3tar_p(ra, rb, ga, gb, alpha, a, f)
143#else
144      subroutine nwxcP_xc_gab3tar(ra, rb, ga, gb, alpha, a, f)
145#endif
146#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
147      subroutine nwxcP_xc_gab3tar_d2(ra, rb, ga, gb, alpha, a, f)
148#else
149      subroutine nwxcP_xc_gab3tar_d3(ra, rb, ga, gb, alpha, a, f)
150#endif
151#include "nwad.fh"
152      implicit none
153      double precision alpha
154#ifdef NWAD_PRINT
155#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
156      type(nwad_dble)::a(3)
157#else
158      double precision a(3)
159#endif
160#else
161      double precision a(3)
162#endif
163      type(nwad_dble):: ra, rb, ga, gb, f
164
165!     B97 like
166!     S. Grimme, J. Comput. Chem., 2006, 27, 1787-1799
167!     doi: 10.1002/jcc.20495
168!     for now, this is only valid for 3 terms.
169
170      type(nwad_dble):: ra83,rb83,ra113,rb113,ra143,rb143,sa2,sb2,
171     +                  sv2,sv2p1,u
172
173      ra83=ra**(8.0d0/3.0d0)
174!     ra113=ra83*ra
175!     ra143=ra113*ra
176      rb83=rb**(8.0d0/3.0d0)
177!     rb113=rb83*rb
178!     rb143=rb113*rb
179      sa2=ga/ra83
180      sb2=gb/rb83
181      sv2=0.5d0*(sa2 + sb2)                  ! called s_av in some papers
182      sv2p1=1.0d0+alpha*sv2
183
184      u=alpha*sv2/sv2p1
185      f=a(1) + a(2)*u + a(3)*u**2
186
187      end
188
189#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
190#ifdef NWAD_PRINT
191      subroutine nwxcP_xc_gab5tar_p(ra, rb, ga, gb, alpha, a, f)
192#else
193      subroutine nwxcP_xc_gab5tar(ra, rb, ga, gb, alpha, a, f)
194#endif
195#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
196      subroutine nwxcP_xc_gab5tar_d2(ra, rb, ga, gb, alpha, a, f)
197#else
198      subroutine nwxcP_xc_gab5tar_d3(ra, rb, ga, gb, alpha, a, f)
199#endif
200#include "nwad.fh"
201      implicit none
202      double precision alpha
203#ifdef NWAD_PRINT
204#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
205      type(nwad_dble)::a(5)
206#else
207      double precision a(5)
208#endif
209#else
210      double precision a(5)
211#endif
212      type(nwad_dble):: ra, rb, ga, gb, f
213
214!     B97 like
215!     S. Grimme, J. Comput. Chem., 2006, 27, 1787-1799
216!     doi: 10.1002/jcc.20495
217!     for now, this is only valid for 5 terms.
218
219      type(nwad_dble):: ra83,rb83,ra113,rb113,ra143,rb143,sa2,sb2,
220     +                  sv2,sv2p1,u
221
222      ra83=ra**(8.0d0/3.0d0)
223!     ra113=ra83*ra
224!     ra143=ra113*ra
225      rb83=rb**(8.0d0/3.0d0)
226!     rb113=rb83*rb
227!     rb143=rb113*rb
228      sa2=ga/ra83
229      sb2=gb/rb83
230      sv2=0.5d0*(sa2 + sb2)                  ! called s_av in some papers
231      sv2p1=1.0d0+alpha*sv2
232
233      u=alpha*sv2/sv2p1
234      f=a(1) + a(2)*u + a(3)*u**2 + a(4)*u**3 + a(5)*u**4
235
236      end
237
238
239!/ =============================================================================
240!/
241!/                      B97 exchange
242!/
243!/ =============================================================================
244
245
246#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
247#ifdef NWAD_PRINT
248      subroutine nwxcp_x_b97_p(rho_a,rho_b,ScalGGAX,tol_rho,FX,sol)
249#else
250      subroutine nwxcp_x_b97(rho_a,rho_b,ScalGGAX,tol_rho,FX,sol)
251#endif
252#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
253      subroutine nwxcp_x_b97_d2(rho_a,rho_b,ScalGGAX,tol_rho,FX,sol)
254#else
255      subroutine nwxcp_x_b97_d3(rho_a,rho_b,ScalGGAX,tol_rho,FX,sol)
256#endif
257#include "nwad.fh"
258      implicit none
259      type(nwad_dble):: rho_a(2), rho_b(2), FX(0:0)
260      double precision ScalGGAX, tol_rho
261#ifdef NWAD_PRINT
262#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
263      type(nwad_dble)::sol(6), xss(5)
264#else
265      double precision sol(6), xss(5)
266#endif
267#else
268      double precision sol(6), xss(5)
269#endif
270      integer max_pow_u, i
271      double precision gamma_xss
272      parameter(gamma_xss = 0.004d0)
273
274      type(nwad_dble):: ra, rb, gaa, gbb
275      type(nwad_dble):: ra43, rb43
276      ! enhancement factors and derivatives thereof
277      type(nwad_dble):: exaa, exbb
278      double precision CX,CX43,CX49
279
280      ra   = rho_a(1)
281      rb   = rho_b(1)
282      gaa  = rho_a(2) ! gamma_aa = G_AA
283      gbb  = rho_b(2) ! gamma_bb = G_BB
284#ifdef NWAD_PRINT
285#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
286      max_pow_u = value(sol(1))
287#else
288      max_pow_u = sol(1)
289#endif
290#else
291      max_pow_u = sol(1)
292#endif
293
294
295      do i = 1, max_pow_u+1
296        xss(i)=sol(i+1)
297      enddo
298
299      !  CX    = 0.930525736349
300      CX = 0.9305257363490993d0
301      ! CX43  = (4.0d0/3.0d0)*CX
302      ! CX49  = (4.0d0/9.0d0)*CX
303
304      if(ra > tol_rho ) then
305        if(max_pow_u == 2 ) then
306#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
307#ifdef NWAD_PRINT
308          call nwxcP_xc_gss3tar_p(ra,gaa,gamma_xss,xss,exaa)
309#else
310          call nwxcP_xc_gss3tar(ra,gaa,gamma_xss,xss,exaa)
311#endif
312#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
313          call nwxcP_xc_gss3tar_d2(ra,gaa,gamma_xss,xss,exaa)
314#else
315          call nwxcP_xc_gss3tar_d3(ra,gaa,gamma_xss,xss,exaa)
316#endif
317        else
318#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
319#ifdef NWAD_PRINT
320          call nwxcP_xc_gss5tar_p(ra,gaa,gamma_xss,xss,exaa)
321#else
322          call nwxcP_xc_gss5tar(ra,gaa,gamma_xss,xss,exaa)
323#endif
324#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
325          call nwxcP_xc_gss5tar_d2(ra,gaa,gamma_xss,xss,exaa)
326#else
327          call nwxcP_xc_gss5tar_d3(ra,gaa,gamma_xss,xss,exaa)
328#endif
329        endif
330        ! Slater-X and derivatives
331        ra43 = -CX  *ra**(4.0d0/3.0d0)
332        ! ra13 = -CX43*ra**(1.0d0/3.0d0)
333        ! raM23= -CX49*ra**(-2.0d0/3.0d0)
334      else
335        ra43 = 0.0d0
336        exaa = 0.0d0
337        ! ra13 = 0.0
338        ! raM23= 0.0
339      endif
340
341      if(rb > tol_rho ) then
342        if(max_pow_u == 2 )then
343#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
344#ifdef NWAD_PRINT
345          call nwxcP_xc_gss3tar_p(rb,gbb,gamma_xss,xss,exbb)
346#else
347          call nwxcP_xc_gss3tar(rb,gbb,gamma_xss,xss,exbb)
348#endif
349#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
350          call nwxcP_xc_gss3tar_d2(rb,gbb,gamma_xss,xss,exbb)
351#else
352          call nwxcP_xc_gss3tar_d3(rb,gbb,gamma_xss,xss,exbb)
353#endif
354        else
355#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
356#ifdef NWAD_PRINT
357          call nwxcP_xc_gss5tar_p(rb,gbb,gamma_xss,xss,exbb)
358#else
359          call nwxcP_xc_gss5tar(rb,gbb,gamma_xss,xss,exbb)
360#endif
361#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
362          call nwxcP_xc_gss5tar_d2(rb,gbb,gamma_xss,xss,exbb)
363#else
364          call nwxcP_xc_gss5tar_d3(rb,gbb,gamma_xss,xss,exbb)
365#endif
366        endif
367        ! Slater-X, more or less
368        rb43 = -CX  *rb**(4.0d0/3.0d0)
369        ! rb13 = -CX43*rb**(1.0d0/3.0d0)
370        ! rbM23= -CX49*rb**(-2.0d0/3.0)
371      else
372        rb43 = 0.0d0
373        exbb = 0.0d0
374        ! rb13 = 0.0
375        ! rbM23 = 0.0
376      endif
377
378      ! ScalGGAX is the ACM parameter for the hybrid
379      ! recall that B97 does not obey the UEG limit
380
381      FX( _FXC_E ) = ScalGGAX*(ra43*exaa + rb43*exbb)
382
383      end
384
385! =============================================================================
386!
387!                      B97 correlation
388!
389! =============================================================================
390
391#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
392#ifdef NWAD_PRINT
393      subroutine nwxcp_c_b97_p(rho_a,rho_b,ScalGGAC,tol_rho,FC,sol)
394#else
395      subroutine nwxcp_c_b97(rho_a,rho_b,ScalGGAC,tol_rho,FC,sol)
396#endif
397#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
398      subroutine nwxcp_c_b97_d2(rho_a,rho_b,ScalGGAC,tol_rho,FC,sol)
399#else
400      subroutine nwxcp_c_b97_d3(rho_a,rho_b,ScalGGAC,tol_rho,FC,sol)
401#endif
402#include "nwad.fh"
403      implicit none
404#include "intf_nwxcP_c_pwlda.fh"
405#include "intf_nwxcP_xc_gab3tar.fh"
406#include "intf_nwxcP_xc_gss3tar.fh"
407#include "intf_nwxcP_xc_gab5tar.fh"
408#include "intf_nwxcP_xc_gss5tar.fh"
409      type(nwad_dble):: rho_a(2), rho_b(2), FC(0:0)
410      double precision ScalGGAC, tol_rho
411#ifdef NWAD_PRINT
412#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
413      type(nwad_dble)::sol(11), css(5), cab(5)
414#else
415      double precision sol(11), css(5), cab(5)
416#endif
417#else
418      double precision sol(11), css(5), cab(5)
419#endif
420
421      type(nwad_dble):: ra,rb,gaa,gbb, ec_a0, ec_b0, ec_ab
422      type(nwad_dble):: gcaa, gcbb, gcab
423      integer max_pow_u, i
424      type(nwad_dble)::  FCLDA(0:0), t
425
426      double precision gamma_css, gamma_cab
427      parameter(gamma_css = 0.2d0, gamma_cab = 0.006d0)
428
429      ra   = rho_a(1)
430      rb   = rho_b(1)
431      gaa  = rho_a(2)
432      gbb  = rho_b(2)
433#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
434#ifdef NWAD_PRINT
435      max_pow_u = value(sol(1))
436#else
437      max_pow_u = sol(1)
438#endif
439#else
440      max_pow_u = sol(1)
441#endif
442
443
444      do i = 1, max_pow_u+1
445        css(i)=sol((i-1)*2+2)
446        cab(i)=sol((i-1)*2+3)
447      enddo
448
449      ! normal
450      if((ra > tol_rho) .or. (rb > tol_rho)) then
451#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
452#ifdef NWAD_PRINT
453        call nwxcP_c_pwlda_p(ra, rb, FCLDA)
454#else
455        call nwxcP_c_pwlda(ra, rb, FCLDA)
456#endif
457#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
458        call nwxcP_c_pwlda_d2(ra, rb, FCLDA)
459#else
460        call nwxcP_c_pwlda_d3(ra, rb, FCLDA)
461#endif
462        ! LDA derivatives
463        ec_ab = FCLDA(_FXC_E)
464      else
465        ec_ab = 0.0d0
466      endif
467
468      ! alpha density only
469      if(ra > tol_rho) then
470        t = 0.0d0
471#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
472#ifdef NWAD_PRINT
473        call nwxcP_c_pwlda_p(ra, t, FCLDA)
474#else
475        call nwxcP_c_pwlda(ra, t, FCLDA)
476#endif
477#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
478        call nwxcP_c_pwlda_d2(ra, t, FCLDA)
479#else
480        call nwxcP_c_pwlda_d3(ra, t, FCLDA)
481#endif
482        ec_a0 = FCLDA(_FXC_E)
483      else
484        ec_a0 = 0.0d0
485      endif
486
487      ! beta density only
488      if(rb > tol_rho) then
489        t = 0.0d0
490#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
491#ifdef NWAD_PRINT
492        call nwxcP_c_pwlda_p(t, rb, FCLDA)
493#else
494        call nwxcP_c_pwlda(t, rb, FCLDA)
495#endif
496#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
497        call nwxcP_c_pwlda_d2(t, rb, FCLDA)
498#else
499        call nwxcP_c_pwlda_d3(t, rb, FCLDA)
500#endif
501        ec_b0 = FCLDA(_FXC_E)
502      else
503        ec_b0 = 0.0d0
504      endif
505
506      ec_ab = ec_ab - ec_a0 - ec_b0
507
508      if(ra > tol_rho) then
509        if(max_pow_u == 2 ) then
510#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
511#ifdef NWAD_PRINT
512          call nwxcP_xc_gss3tar_p(ra,gaa,gamma_css,css,gcaa)
513#else
514          call nwxcP_xc_gss3tar(ra,gaa,gamma_css,css,gcaa)
515#endif
516#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
517          call nwxcP_xc_gss3tar_d2(ra,gaa,gamma_css,css,gcaa)
518#else
519          call nwxcP_xc_gss3tar_d3(ra,gaa,gamma_css,css,gcaa)
520#endif
521        else
522#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
523#ifdef NWAD_PRINT
524          call nwxcP_xc_gss5tar_p(ra,gaa,gamma_css,css,gcaa)
525#else
526          call nwxcP_xc_gss5tar(ra,gaa,gamma_css,css,gcaa)
527#endif
528#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
529          call nwxcP_xc_gss5tar_d2(ra,gaa,gamma_css,css,gcaa)
530#else
531          call nwxcP_xc_gss5tar_d3(ra,gaa,gamma_css,css,gcaa)
532#endif
533        endif
534      else
535        ! gcaa = css(1) + css(2) + css(3) !  u -> 1 if s -> infinity
536        gcaa = 0.0d0 !  Elda(r)*gcaa -> 0 if r -> 0 hence OK to zero gcaa
537      endif
538
539      if(rb > tol_rho) then
540        if(max_pow_u == 2 ) then
541#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
542#ifdef NWAD_PRINT
543          call nwxcP_xc_gss3tar_p(rb,gbb,gamma_css,css,gcbb)
544#else
545          call nwxcP_xc_gss3tar(rb,gbb,gamma_css,css,gcbb)
546#endif
547#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
548          call nwxcP_xc_gss3tar_d2(rb,gbb,gamma_css,css,gcbb)
549#else
550          call nwxcP_xc_gss3tar_d3(rb,gbb,gamma_css,css,gcbb)
551#endif
552        else
553#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
554#ifdef NWAD_PRINT
555          call nwxcP_xc_gss5tar_p(rb,gbb,gamma_css,css,gcbb)
556#else
557          call nwxcP_xc_gss5tar(rb,gbb,gamma_css,css,gcbb)
558#endif
559#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
560          call nwxcP_xc_gss5tar_d2(rb,gbb,gamma_css,css,gcbb)
561#else
562          call nwxcP_xc_gss5tar_d3(rb,gbb,gamma_css,css,gcbb)
563#endif
564        endif
565      else
566        ! gcbb = css(1) + css(2) + css(3) !  u -> 1 if s -> infinity
567        gcbb = 0.0d0 !  Elda(r)*gcbb -> 0 if r -> 0 hence OK to zero gcbb
568      endif
569
570      if(ra > tol_rho .and. rb > tol_rho) then
571        if(max_pow_u == 2 ) then
572#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
573#ifdef NWAD_PRINT
574          call nwxcP_xc_gab3tar_p(ra,rb,gaa,gbb,gamma_cab,cab,gcab)
575#else
576          call nwxcP_xc_gab3tar(ra,rb,gaa,gbb,gamma_cab,cab,gcab)
577#endif
578#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
579          call nwxcP_xc_gab3tar_d2(ra,rb,gaa,gbb,gamma_cab,cab,gcab)
580#else
581          call nwxcP_xc_gab3tar_d3(ra,rb,gaa,gbb,gamma_cab,cab,gcab)
582#endif
583        else
584#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
585#ifdef NWAD_PRINT
586          call nwxcP_xc_gab5tar_p(ra,rb,gaa,gbb,gamma_cab,cab,gcab)
587#else
588          call nwxcP_xc_gab5tar(ra,rb,gaa,gbb,gamma_cab,cab,gcab)
589#endif
590#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
591          call nwxcP_xc_gab5tar_d2(ra,rb,gaa,gbb,gamma_cab,cab,gcab)
592#else
593          call nwxcP_xc_gab5tar_d3(ra,rb,gaa,gbb,gamma_cab,cab,gcab)
594#endif
595        endif
596      else
597        ! This term is not necessary for exchange as there is no cross
598        ! term between the spin-channels in that case.
599        ! gcab = cab(1) + cab(2) + cab(3) !  u == 1 if s->inf.
600        gcab = 0.0d0 !  Elda(r)*gcab -> 0 if r -> 0 hence OK to zero gcab
601      endif
602
603      ! ScalGGAC is the ACM D parameter as used in double hybrids
604
605      FC( _FXC_E ) = ScalGGAC*(ec_a0*gcaa + ec_b0*gcbb + ec_ab*gcab)
606
607      end
608#ifndef NWAD_PRINT
609#define NWAD_PRINT
610c
611c     Compile source again for Maxima
612c
613#include "nwxcP_xc_b97.F"
614#endif
615#ifndef SECOND_DERIV
616#define SECOND_DERIV
617c
618c     Compile source again for the 2nd derivative case
619c
620#include "nwxcP_xc_b97.F"
621#endif
622#ifndef THIRD_DERIV
623#define THIRD_DERIV
624c
625c     Compile source again for the 3rd derivative case
626c
627#include "nwxcP_xc_b97.F"
628#endif
629#undef NWAD_PRINT
630!>
631!> @}
632!>
633!> $Id$
634