1c     PBE exchange functional
2c
3c     References:
4c     [a] J.P. Perdew, K. Burke, and M. Ernzerhof, PRL 77, 3865 (1996).
5c     [b] J.P. Perdew and Y. Wang, Phys. Rev. B 33, 8800 (1986).;
6c                                               40, 3399 (1989) (E).
7c     Hammer, Hansen and Norskov, PRB 59, 7413 (1999) [RPBE]
8c     Zhang and Yang, PRL 80, 890 (1998) [RevPBE]
9c
10#if !defined SECOND_DERIV && !defined THIRD_DERIV
11      Subroutine xc_xpbe96(whichf,
12     W     tol_rho, fac, lfac, nlfac, rho, delrho,
13     &                     Amat, Cmat, nq, ipol, Ex, qwght,ldew,func)
14#elif defined(SECOND_DERIV) && !defined THIRD_DERIV
15      Subroutine xc_xpbe96_d2(whichf,
16     W     tol_rho, fac, lfac, nlfac, rho, delrho,
17     &                        Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
18     &                        qwght,ldew,func)
19#else
20      Subroutine xc_xpbe96_d3(whichf,
21     W     tol_rho, fac, lfac, nlfac, rho, delrho,
22     &                        Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3,
23     &                        nq, ipol, Ex, qwght,ldew,func)
24#endif
25c
26      implicit none
27c
28#include "dft2drv.fh"
29#include "dft3drv.fh"
30c
31      character*4 whichf
32      double precision fac, Ex
33      integer nq, ipol
34      logical lfac, nlfac,ldew
35      double precision func(*)  ! value of the functional [output]
36c
37c     Charge Density & Its Cube Root
38c
39      double precision rho(nq,ipol*(ipol+1)/2)
40c
41c     Charge Density Gradient
42c
43      double precision delrho(nq,3,ipol)
44c
45c     Quadrature Weights
46c
47      double precision qwght(nq)
48c
49c     Sampling Matrices for the XC Potential & Energy
50c
51      double precision amat(nq,ipol), cmat(nq,*)
52c
53#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
54      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
55#endif
56#ifdef THIRD_DERIV
57      double precision Amat3(nq,NCOL_AMAT3), Cmat3(nq,NCOL_CMAT3)
58#endif
59c
60      double precision tol_rho, pi, um, uk, umk,ukrev,umkrev
61      double precision C, Cs
62      double precision F43, F13, F23
63c
64#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
65      double precision F73
66#endif
67#ifdef THIRD_DERIV
68      double precision F10d3
69#endif
70      parameter(um=0.2195149727645171d0, uk=0.8040d0, umk=um/uk)
71      parameter(ukrev=1.245d0, umkrev=um/ukrev)
72      parameter (F43=4.d0/3.d0, F13=1.d0/3.d0, F23=2.0d0/3.0d0)
73#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
74      parameter (F73=7.d0/3.d0)
75#endif
76#ifdef THIRD_DERIV
77      parameter (F10d3=10.0d0/3.0d0)
78#endif
79c
80      integer n
81      double precision rrho, rho43, rho13, gamma, gam12, s, d1s(2),
82     &      d, g, gp, d1g(2)
83#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
84      double precision rhom23, d2s(3), gpp, d2g(3)
85#endif
86#ifdef THIRD_DERIV
87      double precision d3s(4), d3g(4), rhom53, gppp
88#endif
89      double precision gpbe0, gpbe1, gpbe2, gpbe3
90      double precision grpbe0, grpbe1, grpbe2, grpbe3
91      double precision grevpbe0, grevpbe1, grevpbe2, grevpbe3
92c Original PBE
93      gpbe0(s)= uk*(1d0 - 1d0/(1d0+umk*s*s))
94      gpbe1(s)= 2d0*um*s/(1d0+umk*s*s)**2
95      gpbe2(s)= 2d0*um*(1d0-4d0*umk*s*s/(1d0+umk*s*s))/(1d0+umk*s*s)**2
96      gpbe3(s)= 24.0d0*umk*um*s*
97     1  (2.0d0*umk*s*s/(1.0d0+umk*s*s)-1.0d0)/(1.0d0+umk*s*s)**3
98c revPBE by Zhang et al.
99      grevpbe0(s)= ukrev*(1d0 - 1d0/(1d0+umkrev*s*s))
100      grevpbe1(s)= 2d0*um*s/(1d0+umkrev*s*s)**2
101      grevpbe2(s)= 2d0*um*(1d0-4d0*umkrev*s*s/(1d0+umkrev*s*s))/
102     /     (1d0+umkrev*s*s)**2
103      grevpbe3(s)= 24.0d0*umkrev*um*s*
104     1  (2.0d0*umkrev*s*s/(1.0d0+umkrev*s*s)-1.0d0)/
105     2  (1.0d0+umkrev*s*s)**3
106c RPBE by Hammer et al.
107      grpbe0(s)= uk*(1d0 - exp(-umk*s*s))
108      grpbe1(s)= 2d0*um*s*exp(-umk*s*s)
109      grpbe2(s)= 2d0*um*exp(-umk*s*s)*(1d0-2d0*umk*s*s)
110      grpbe3(s)= -4.0d0*umk*um*s*exp(-umk*s*s)*(3d0-2d0*umk*s*s)
111c
112      pi = acos(-1.d0)
113      C = -3d0/(4d0*pi)*(3d0*pi*pi)**F13
114      Cs = 0.5d0/(3d0*pi*pi)**F13
115      Cs = Cs * C               ! account for including C in rho43
116c
117      if (ipol.eq.1 )then
118c
119c        ======> SPIN-RESTRICTED <======
120c
121c Daniel (9-28-12): There are somewhat mysterious coefficients involved
122c in the evaluation of the functional and its derivatives.  We must
123c recall that the exchange energy is always written based on the
124c spin-scaling relationship for exchange:
125c
126c Ex[rho] = Ex[rho_a,rho_b] = 0.5*( Ex[2*rho_a] + Ex[2*rho_b] )
127c
128c Thus, the electron density is always written:
129c rho -> 2*rho_s
130c gamma -> 4*gamma_ss
131c
132c Rationalization for the coefficients is mathematically justified below:
133c
134c ----------------------------
135c Amat       -> 0.5*2 = 1
136c Cmat       -> 0.5*4 = 2
137c ----------------------------
138c Amat2      -> 0.5*2*2 = 2
139c Cmat2(rg)  -> 0.5*2*4 = 4
140c Cmat2(gg)  -> 0.5*4*4 = 8
141c ----------------------------
142c Amat3      -> 0.5*2*2*2 = 4
143c Cmat3(rrg) -> 0.5*2*2*4 = 8
144c Cmat3(rgg) -> 0.5*2*4*4 = 16
145c Cmat3(ggg) -> 0.5*4*4*4 = 32
146c ----------------------------
147c
148c If, instead, the author of this code had decided to divide the total
149c density (rho(n,1)) by 2 in constructing the density and gamma, those
150c coefficients would be unnecessary.
151c
152#ifdef IFCV81
153CDEC$ NOSWP
154#endif
155         do 10 n = 1, nq
156            if (rho(n,1).lt.tol_rho) goto 10
157c#ifdef THIRD_DERIV
158c            write(6,*) 'rho', rho(n,1)
159c#endif
160            rho43 = C*rho(n,1)**F43
161            rrho = 1d0/rho(n,1)
162            rho13 = F43*rho43*rrho
163#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
164            rhom23 = F13*rho13*rrho
165#endif
166#ifdef THIRD_DERIV
167            rhom53 = F23*rhom23*rrho
168#endif
169            if (lfac) then
170               Ex = Ex + rho43*qwght(n)*fac
171               if(ldew)func(n) = func(n) + rho43*fac
172               Amat(n,1) = Amat(n,1) + rho13*fac
173#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
174               Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 2d0*rhom23*fac
175#endif
176#ifdef THIRD_DERIV
177               Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA)
178     1                              - 4.0d0*rhom53*fac
179#endif
180            endif
181c
182            gamma = delrho(n,1,1)*delrho(n,1,1) +
183     &              delrho(n,2,1)*delrho(n,2,1) +
184     &              delrho(n,3,1)*delrho(n,3,1)
185            gam12 = dsqrt(gamma)
186            if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 10
187c
188
189            s = Cs*gam12/rho43
190            d1s(1) = -F43*s*rrho
191            d1s(2) = 0.5d0*s/gamma
192c
193c     Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1
194c
195            if(whichf.eq.'revp') then
196               g=grevpbe0(s)
197               gp=grevpbe1(s)
198            elseif(whichf.eq.'rpbe') then
199               g=grpbe0(s)
200               gp=grpbe1(s)
201            else
202               g=gpbe0(s)
203               gp=gpbe1(s)
204            endif
205c
206c Daniel (7-27-12): gp is the derivative of the rational function,
207c or whatever the function in the revision is.
208c First derivatives of the enhancement factor
209            d1g(1) = gp*d1s(1)
210            d1g(2) = gp*d1s(2)
211            Ex = Ex + rho43*g*qwght(n)*fac
212            if(ldew)func(n) = func(n) + rho43*g*fac
213            Amat(n,1) = Amat(n,1) + (rho13*g+rho43*d1g(1))*fac
214            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 2d0*rho43*d1g(2)*fac
215#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
216            d2s(1) = -F73*d1s(1)*rrho
217            d2s(2) = -F43*d1s(2)*rrho
218            d2s(3) = -0.5d0*d1s(2)/gamma
219            if(whichf.eq.'revp') then
220               gpp=grevpbe2(s)
221            elseif(whichf.eq.'rpbe') then
222               gpp=grpbe2(s)
223            else
224               gpp=gpbe2(s)
225            endif
226c Second derivatives of the enhancement factor
227            d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1)
228            d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2)
229            d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2)
230            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
231     &           +(rhom23*g
232     &           + 2.d0*rho13*d1g(1)
233     &           + rho43*d2g(1))*fac*2d0
234            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA)
235     &           +(rho13*d1g(2)
236     &           + rho43*d2g(2))*fac*4d0
237            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA)
238     &           + rho43*d2g(3)*fac*8d0
239#endif
240#ifdef THIRD_DERIV
241c 1 = drdrdr, 2 = drdrdg, 3 = drdgdg, 4 = dgdgdg
242            d3s(1) = -F10d3*d2s(1)*rrho
243            d3s(2) = 0.5d0*d2s(1)/gamma
244            d3s(3) = -F43*d2s(3)*rrho
245            d3s(4) = -1.5d0*d2s(3)/gamma
246            if(whichf.eq.'revp') then
247               gppp = grevpbe3(s)
248            elseif(whichf.eq.'rpbe') then
249               gppp = grpbe3(s)
250            else
251               gppp = gpbe3(s)
252            endif
253c Third derivatives of the enhancement factor
254            d3g(1) = gp*d3s(1) + 3.0d0*gpp*d1s(1)*d2s(1)
255     1             + gppp*d1s(1)*d1s(1)*d1s(1)
256            d3g(2) = gp*d3s(2)
257     1             + gpp*d1s(2)*d2s(1)
258     2             + 2.0d0*gpp*d1s(1)*d2s(2)
259     3             + gppp*d1s(1)*d1s(1)*d1s(2)
260            d3g(3) = gp*d3s(3)
261     1             + gpp*d1s(1)*d2s(3)
262     2             + 2.0d0*gpp*d1s(2)*d2s(2)
263     3             + gppp*d1s(1)*d1s(2)*d1s(2)
264            d3g(4) = gp*d3s(4) + 3.0d0*gpp*d1s(2)*d2s(3)
265     1             + gppp*d1s(2)*d1s(2)*d1s(2)
266c
267            Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA)
268     1           + (-rhom53*g
269     2           +  3.0d0*rhom23*d1g(1)
270     3           +  3.0d0*rho13*d2g(1)
271     4           +  rho43*d3g(1))*fac*4.0d0
272            Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA)
273     1           + (rhom23*d1g(2)
274     2           +  2.0d0*rho13*d2g(2)
275     3           +  rho43*d3g(2))*fac*8.0d0
276            Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA)
277     1           + (rho13*d2g(3)
278     2           +  rho43*d3g(3))*fac*16.0d0
279            Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA)
280     1           + (rho43*d3g(4))*fac*32.0d0
281#endif
282 10      continue
283c
284      else
285c
286c        ======> SPIN-UNRESTRICTED <======
287c
288c Daniel (9-28-12): There are somewhat mysterious coefficients involved
289c in the evaluation of the functional and its derivatives.  We must
290c recall that the exchange energy is always written based on the
291c spin-scaling relationship for exchange:
292c
293c Ex[rho] = Ex[rho_a,rho_b] = 0.5*( Ex[2*rho_a] + Ex[2*rho_b] )
294c
295c Thus, the electron density is always written:
296c rho -> 2*rho_s
297c gamma -> 4*gamma_ss
298c
299c It seems like the derivatives should be correctly balanced by the
300c following coefficients:
301c
302c -----------------------------
303c Amat       -> 0.5*2 = 1
304c Cmat       -> 0.5*1 = 0.5
305c -----------------------------
306c Amat2      -> 0.5*2*2 = 2
307c Cmat2(rg)  -> 0.5*2*1 = 1
308c Cmat2(gg)  -> 0.5*1*1 = 0.5
309c -----------------------------
310c Amat3      -> 0.5*2*2*2 = 4
311c Cmat3(rrg) -> 0.5*2*2*1 = 2
312c Cmat3(rgg) -> 0.5*2*1*1 = 1
313c Cmat3(ggg) -> 0.5*1*1*1 = 0.5
314c -----------------------------
315c
316#ifdef IFCV81
317CDEC$ NOSWP
318#endif
319         do 20 n = 1, nq
320            if (rho(n,1).lt.tol_rho) goto 20
321c
322c     Alpha
323c
324            if (rho(n,2).lt.tol_rho) goto 25
325            rho43 = C*(2d0*rho(n,2))**F43
326            rrho = 0.5d0/rho(n,2)
327            rho13 = F43*rho43*rrho
328#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
329            rhom23 = F13*rho13*rrho
330#endif
331#ifdef THIRD_DERIV
332            rhom53 = F23*rhom23*rrho
333#endif
334            if (lfac) then
335               Ex = Ex + rho43*qwght(n)*fac*0.5d0
336               if(ldew)func(n) = func(n) + rho43*fac*0.5d0
337               Amat(n,1) = Amat(n,1) + rho13*fac
338#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
339               Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 2d0*rhom23*fac
340#endif
341#ifdef THIRD_DERIV
342               Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA)
343     1                              - 4.0d0*rhom53*fac
344#endif
345            endif
346c
347            gamma = delrho(n,1,1)*delrho(n,1,1) +
348     &              delrho(n,2,1)*delrho(n,2,1) +
349     &              delrho(n,3,1)*delrho(n,3,1)
350            gam12 = 2d0*dsqrt(gamma)
351            if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 25
352c
353            s = Cs*gam12/rho43
354            d1s(1) = -F43*s*rrho
355            d1s(2) = 0.5d0*s/gamma
356c
357c     Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1
358c
359
360            if(whichf.eq.'revp') then
361               g=grevpbe0(s)
362               gp=grevpbe1(s)
363            elseif(whichf.eq.'rpbe') then
364               g=grpbe0(s)
365               gp=grpbe1(s)
366            else
367               g=gpbe0(s)
368               gp=gpbe1(s)
369            endif
370c Daniel (9-28-12): Factors of 2 are inconsistent with the restricted
371c calculations because a gam12 is doubled above and
372            d1g(1) = gp*d1s(1)
373            d1g(2) = gp*d1s(2)
374            Ex = Ex + rho43*g*qwght(n)*fac*0.5d0
375            if(ldew)func(n) = func(n) + rho43*g*fac*0.5d0
376            Amat(n,1) = Amat(n,1) + (rho13*g+rho43*d1g(1))*fac
377            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 0.5d0*rho43*d1g(2)*fac
378#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
379            d2s(1) = -F73*d1s(1)*rrho
380            d2s(2) = -F43*d1s(2)*rrho
381            d2s(3) = -0.5d0*d1s(2)/gamma
382            if(whichf.eq.'revp') then
383               gpp=grevpbe2(s)
384            elseif(whichf.eq.'rpbe') then
385               gpp=grpbe2(s)
386            else
387               gpp=gpbe2(s)
388            endif
389            d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1)
390            d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2)
391            d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2)
392c
393            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
394     &           +(rhom23*g
395     &           + 2.d0*rho13*d1g(1)
396     &           + rho43*d2g(1))*fac*2d0
397            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA)
398     &           +(rho13*d1g(2)
399     &           + rho43*d2g(2))*fac
400            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA)
401     &           + rho43*d2g(3)*fac*0.5d0
402#endif
403#ifdef THIRD_DERIV
404c 1 = drdrdr, 2 = drdrdg, 3 = drdgdg, 4 = dgdgdg
405            d3s(1) = -F10d3*d2s(1)*rrho
406            d3s(2) = 0.5d0*d2s(1)/gamma
407            d3s(3) = -F43*d2s(3)*rrho
408            d3s(4) = -1.5d0*d2s(3)/gamma
409            if(whichf.eq.'revp') then
410               gppp = grevpbe3(s)
411            elseif(whichf.eq.'rpbe') then
412               gppp = grpbe3(s)
413            else
414               gppp = gpbe3(s)
415            endif
416c Third derivatives of the enhancement factor
417            d3g(1) = gp*d3s(1) + 3.0d0*gpp*d1s(1)*d2s(1)
418     1             + gppp*d1s(1)*d1s(1)*d1s(1)
419            d3g(2) = gp*d3s(2)
420     1             + gpp*d1s(2)*d2s(1)
421     2             + 2.0d0*gpp*d1s(1)*d2s(2)
422     3             + gppp*d1s(1)*d1s(1)*d1s(2)
423            d3g(3) = gp*d3s(3)
424     1             + gpp*d1s(1)*d2s(3)
425     2             + 2.0d0*gpp*d1s(2)*d2s(2)
426     3             + gppp*d1s(1)*d1s(2)*d1s(2)
427            d3g(4) = gp*d3s(4) + 3.0d0*gpp*d1s(2)*d2s(3)
428     1             + gppp*d1s(2)*d1s(2)*d1s(2)
429c
430            Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA)
431     1           + (-rhom53*g
432     2           +  3.0d0*rhom23*d1g(1)
433     3           +  3.0d0*rho13*d2g(1)
434     4           +  rho43*d3g(1))*fac*4.0d0
435            Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA)
436     1           + (rhom23*d1g(2)
437     2           +  2.0d0*rho13*d2g(2)
438     3           +  rho43*d3g(2))*fac*2.0d0
439            Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA)
440     1           + (rho13*d2g(3)
441     2           +  rho43*d3g(3))*fac
442            Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA)
443     1           + (rho43*d3g(4))*fac*0.5d0
444#endif
445c
446c     Beta
447c
448 25         continue
449            if (rho(n,3).lt.tol_rho) goto 20
450            rho43 = C*(2d0*rho(n,3))**F43
451            rrho = 0.5d0/rho(n,3)
452            rho13 = F43*rho43*rrho
453#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
454            rhom23 = F13*rho13*rrho
455#endif
456#ifdef THIRD_DERIV
457            rhom53 = F23*rhom23*rrho
458#endif
459            if (lfac) then
460               Ex = Ex + rho43*qwght(n)*fac*0.5d0
461               if(ldew)func(n) = func(n) + rho43*fac*0.5d0
462               Amat(n,2) = Amat(n,2) + rho13*fac
463#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
464               Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + 2d0*rhom23*fac
465#endif
466#ifdef THIRD_DERIV
467               Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB)
468     1                              - 4.0d0*rhom53*fac
469#endif
470            endif
471c
472            gamma = delrho(n,1,2)*delrho(n,1,2) +
473     &              delrho(n,2,2)*delrho(n,2,2) +
474     &              delrho(n,3,2)*delrho(n,3,2)
475            gam12 = 2d0*dsqrt(gamma)
476            if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 20
477c
478            s = Cs*gam12/rho43
479            d1s(1) = -F43*s*rrho
480            d1s(2) = 0.5d0*s/gamma
481c
482c     Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1
483c
484            if(whichf.eq.'revp') then
485               g=grevpbe0(s)
486               gp=grevpbe1(s)
487            elseif(whichf.eq.'rpbe') then
488               g=grpbe0(s)
489               gp=grpbe1(s)
490            else
491               g=gpbe0(s)
492               gp=gpbe1(s)
493            endif
494c
495            d1g(1) = gp*d1s(1)
496            d1g(2) = gp*d1s(2)
497            Ex = Ex + rho43*g*qwght(n)*fac*0.5d0
498            if(ldew)func(n) = func(n) + rho43*g*fac*0.5d0
499            Amat(n,2) = Amat(n,2) + (rho13*g+rho43*d1g(1))*fac
500            Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + 0.5d0*rho43*d1g(2)*fac
501#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
502            d2s(1) = -F73*d1s(1)*rrho
503            d2s(2) = -F43*d1s(2)*rrho
504            d2s(3) = -0.5d0*d1s(2)/gamma
505            if(whichf.eq.'revp') then
506               gpp=grevpbe2(s)
507            elseif(whichf.eq.'rpbe') then
508               gpp=grpbe2(s)
509            else
510               gpp=gpbe2(s)
511            endif
512            d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1)
513            d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2)
514            d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2)
515            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB)
516     &           +(rhom23*g
517     &           + 2.d0*rho13*d1g(1)
518     &           + rho43*d2g(1))*fac*2d0
519            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB)
520     &           +(rho13*d1g(2)
521     &           + rho43*d2g(2))*fac
522            Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB)
523     &           + rho43*d2g(3)*fac*0.5d0
524#endif
525#ifdef THIRD_DERIV
526c 1 = drdrdr, 2 = drdrdg, 3 = drdgdg, 4 = dgdgdg
527            d3s(1) = -F10d3*d2s(1)*rrho
528            d3s(2) = 0.5d0*d2s(1)/gamma
529            d3s(3) = -F43*d2s(3)*rrho
530            d3s(4) = -1.5d0*d2s(3)/gamma
531            if(whichf.eq.'revp') then
532               gppp = grevpbe3(s)
533            elseif(whichf.eq.'rpbe') then
534               gppp = grpbe3(s)
535            else
536               gppp = gpbe3(s)
537            endif
538c Third derivatives of the enhancement factor
539            d3g(1) = gp*d3s(1) + 3.0d0*gpp*d1s(1)*d2s(1)
540     1             + gppp*d1s(1)*d1s(1)*d1s(1)
541            d3g(2) = gp*d3s(2)
542     1             + gpp*d1s(2)*d2s(1)
543     2             + 2.0d0*gpp*d1s(1)*d2s(2)
544     3             + gppp*d1s(1)*d1s(1)*d1s(2)
545            d3g(3) = gp*d3s(3)
546     1             + gpp*d1s(1)*d2s(3)
547     2             + 2.0d0*gpp*d1s(2)*d2s(2)
548     3             + gppp*d1s(1)*d1s(2)*d1s(2)
549            d3g(4) = gp*d3s(4) + 3.0d0*gpp*d1s(2)*d2s(3)
550     1             + gppp*d1s(2)*d1s(2)*d1s(2)
551c
552            Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB)
553     1           + (-rhom53*g
554     2           +  3.0d0*rhom23*d1g(1)
555     3           +  3.0d0*rho13*d2g(1)
556     4           +  rho43*d3g(1))*fac*4.0d0
557            Cmat3(n,D3_RB_RB_GBB) = Cmat3(n,D3_RB_RB_GBB)
558     1           + (rhom23*d1g(2)
559     2           +  2.0d0*rho13*d2g(2)
560     3           +  rho43*d3g(2))*fac*2.0d0
561            Cmat3(n,D3_RB_GBB_GBB) = Cmat3(n,D3_RB_GBB_GBB)
562     1           + (rho13*d2g(3)
563     2           +  rho43*d3g(3))*fac
564            Cmat3(n,D3_GBB_GBB_GBB) = Cmat3(n,D3_GBB_GBB_GBB)
565     1           + (rho43*d3g(4))*fac*0.5d0
566#endif
567c
568 20      continue
569      endif
570c
571      return
572      end
573#ifndef SECOND_DERIV
574#define SECOND_DERIV
575c
576c     Compile source again for the 2nd derivative case
577c
578#include "xc_xpbe96.F"
579#endif
580#ifndef THIRD_DERIV
581#define THIRD_DERIV
582c
583c     Compile source again for the 3rd derivative case
584c
585#include "xc_xpbe96.F"
586#endif
587c $Id$
588