1#ifndef SECOND_DERIV
2      Subroutine xc_spbe96(tol_rho, cfac, lcfac, nlcfac, rho, delrho,
3     &                    Amat, Cmat, nq, ipol, Ec, qwght,
4     &                    ldew,ffunc)
5#else
6      Subroutine xc_spbe96_d2(tol_rho, cfac, lcfac, nlcfac,
7     &                    rho, delrho, Amat, Amat2, Cmat, Cmat2,
8     &                    nq, ipol, Ec, qwght, ldew, ffunc)
9#endif
10c
11c$Id$
12c
13      Implicit none
14#include "dft2drv.fh"
15c
16c
17c     Input and other parameters
18c
19      integer lnq
20      integer ipol, nq
21      double precision dummy(1)
22      double precision cfac(*)
23      logical lcfac(*), nlcfac(*), ldew
24      logical lfac, nlfac, lfacl, nlfacl
25      double precision ffunc(*)
26      double precision fac, facl
27      double precision lqwght
28      double precision tol_rho
29c
30c     Constants in PBE functional
31c
32      double precision GAMMA, BETA, PI
33      parameter (GAMMA = 0.03109069086965489503494086371273d0)
34      parameter (BETA = 0.06672455060314922d0)
35      parameter (PI = 3.1415926535897932385d0)
36c
37c     Threshold parameters
38c
39      double precision TOLL, EXPTOL
40      double precision EPS
41      parameter (TOLL = 1.0D-40, EXPTOL = 40.0d0)
42      parameter (EPS = 1.0e-8)
43c
44c     Correlation energy
45c
46      double precision Ec
47c
48c     Charge Density
49c
50      double precision rho(nq,ipol*(ipol+1)/2)
51      double precision rho_t(3)
52c
53c     Charge Density Gradient
54c
55      double precision delrho(nq,3,ipol)
56      double precision dsqgamma
57c
58c     Quadrature Weights
59c
60      double precision qwght(nq)
61c
62c     Sampling Matrices for the XC Potential
63c
64      double precision Amat(nq,ipol), Cmat(nq,*)
65#ifdef SECOND_DERIV
66c
67c     Sampling Matrices for the XC Kernel
68c
69      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
70#endif
71c
72c     Intermediate derivative results, etc.
73c
74      integer n
75      double precision rhoval, gammaval
76      double precision nepsc, dnepscdn(2)
77      double precision epsc, depscdna, depscdnb
78      double precision H0, dH0dna, dH0dnb, dH0dg
79      double precision phi, dphidna, dphidnb, dphidzeta
80      double precision zeta, dzetadna, dzetadnb
81      double precision arglog, darglogdna, darglogdnb, darglogdg
82      double precision fAt, dfAtdt, dfAtdA
83      double precision fAtnum
84      double precision fAtden, dfAtdendt, dfAtdendA
85      double precision dfAtdna, dfAtdnb, dfAtdg
86      double precision A, dAdna, dAdnb
87      double precision t, dtdna, dtdnb, dtdg
88      double precision ks, dksdna, dksdnb
89      double precision argexp, dargexpdna, dargexpdnb
90      double precision expinA
91#ifdef SECOND_DERIV
92      double precision d2nepscdn2(NCOL_AMAT2)
93      double precision d2epscdna2, d2epscdnadnb, d2epscdnb2
94      double precision d2H0dna2, d2H0dnadnb, d2H0dnb2
95      double precision d2H0dnadg, d2H0dnbdg, d2H0dg2
96      double precision d2phidzeta2, d2phidna2, d2phidnadnb, d2phidnb2
97      double precision d2zetadna2, d2zetadnadnb, d2zetadnb2
98      double precision d2arglogdna2, d2arglogdnb2, d2arglogdnadnb
99      double precision d2arglogdnadg, d2arglogdnbdg, d2arglogdg2
100      double precision d2fAtdt2, d2fAtdA2, d2fAtdtdA, d2fAtdg2
101      double precision d2fAtdendt2, d2fAtdendtdA, d2fAtdendA2
102      double precision d2fAtdna2, d2fAtdnb2, d2fAtdnadnb
103      double precision d2fAtdnadg, d2fAtdnbdg
104      double precision d2Adna2, d2Adnadnb, d2Adnb2
105      double precision d2tdna2, d2tdnb2, d2tdnadnb
106      double precision d2tdg2, d2tdnadg, d2tdnbdg
107      double precision d2ksdna2, d2ksdnb2, d2ksdnadnb
108      double precision d2argexpdna2, d2argexpdnb2, d2argexpdnadnb
109#endif
110c
111c References:
112c [a] J. P. Perdew, K. Burke, and M. Ernzerhof,
113c     {\it Generalized gradient approximation made simple},
114c     Phys.\ Rev.\ Lett. {\bf 77,} 3865 (1996).
115c [b] J. P. Perdew, K. Burke, and Y. Wang, {\it Real-space cutoff
116c     construction of a generalized gradient approximation: The PW91
117c     density functional}, submitted to Phys.\ Rev.\ B, Feb. 1996.
118c [c] J. P. Perdew and Y. Wang, Phys.\ Rev.\ B {\bf 45}, 13244 (1992).
119c
120c  E_c(PBE) = Int n (epsilon_c + H0) dxdydz
121c
122c  n*epsilon_c                <=== supplied by another subroutine
123c  d(n*epsilon_c)/d(na)       <=== supplied by another subroutine
124c  d2(n*epsilon_c)/d(na)d(na) <=== supplied by another subroutine
125c  d2(n*epsilon_c)/d(na)d(nb) <=== supplied by another subroutine
126c  d2(n*epsilon_c)/d(nb)d(nb) <=== supplied by another subroutine
127c
128c  H0 = GAMMA * phi**3 * log{ 1 + BETA/GAMMA * t**2 * [ ... ]}
129c
130c  phi = (1/2)[(1+zeta)**(2/3)+(1-zeta)**(2/3)]
131c
132c  zeta = (na - nb)/n
133c
134c  [ ... ] = (1 + A * t**2)/(1 + A * t**2 + A**2 * t**4)
135c
136c  A = BETA/GAMMA [exp{-epsilon_c/(GAMMA*phi**3)}-1]**(-1)
137c
138c  t = |Nabla n|/(2*phi*ks*n)
139c
140c  ks = 2 * (3 * PI**2 * n)**(1/6) / sqrt(PI)
141c
142c  |Nabla n| = sqrt(g_aa + g_bb + 2*g_ab)
143c
144c  Names of variables
145c
146c  E_c(PBE)                  : Ec
147c  n (alpha+beta density)    : rhoval
148c  na, nb                    : rho(*,2), rho(*,3)
149c  epsilon_c                 : epsc
150c  H0                        : H0
151c  n*epsilon_c               : nepsc
152c  phi                       : phi
153c  zeta                      : zeta
154c  { ... }                   : arglog
155c  [ ... ]                   : fAt
156c  (1 + A * t**2)            : fAtnum
157c  (1 + A * t**2 + A**2 * t**4) : fAtden
158c  A                         : A
159c  t                         : t
160c  |Nabla n|                 : gammaval
161c  ks                        : ks
162c  {-epsilon_c ... }         : argexp
163c  g_aa, g_bb, g_ab          : g
164c
165c  Derivatives of these are named like d...dna, d2...dnadnb,
166c  d2...dna2, etc.
167c
168      fac = cfac(46)
169      lfac = lcfac(46)
170      nlfac = nlcfac(46)
171c
172c     ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <======
173c
174      do 20 n = 1, nq
175c
176c        n and zeta = (na - nb)/n
177c
178         rhoval = rho(n,1)
179         rho_t(1) = rho(n,1)
180         if (ipol.eq.2) then
181            rho_t(2) = rho(n,2)
182            rho_t(3) = rho(n,3)
183         endif
184         if (rhoval.le.tol_rho) goto 20
185         if (ipol.eq.1) then
186            gammaval = delrho(n,1,1)*delrho(n,1,1) +
187     &                 delrho(n,2,1)*delrho(n,2,1) +
188     &                 delrho(n,3,1)*delrho(n,3,1)
189         else
190            gammaval = delrho(n,1,1)*delrho(n,1,1) +
191     &                 delrho(n,1,2)*delrho(n,1,2) +
192     &                 delrho(n,2,1)*delrho(n,2,1) +
193     &                 delrho(n,2,2)*delrho(n,2,2) +
194     &                 delrho(n,3,1)*delrho(n,3,1) +
195     &                 delrho(n,3,2)*delrho(n,3,2) +
196     &           2.d0*(delrho(n,1,1)*delrho(n,1,2) +
197     &                 delrho(n,2,1)*delrho(n,2,2) +
198     &                 delrho(n,3,1)*delrho(n,3,2))
199         endif
200         dsqgamma = max(dsqrt(gammaval),tol_rho)
201         nepsc = 0.0d0
202         dnepscdn(1) = 0.0d0
203         if (ipol.eq.2) dnepscdn(2) = 0.0d0
204#ifdef SECOND_DERIV
205         d2nepscdn2(D2_RA_RA)=0.0d0
206         d2nepscdn2(D2_RA_RB)=0.0d0
207         if (ipol.eq.2) d2nepscdn2(D2_RB_RB)=0.0d0
208#endif
209c
210c        call for LDA bit
211c
212         lnq = 1
213         lqwght = 1.0d0
214c
215         if (abs(cfac(1)).gt.EPS) then
216            facl = cfac(1)
217#ifndef SECOND_DERIV
218            call xc_vwn_5(tol_rho,facl,rho_t,
219     &         dnepscdn,lnq,ipol,nepsc,lqwght,.false.,dummy)
220#else
221            call xc_vwn_5_d2(tol_rho,facl,rho_t,
222     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
223     &         .false.,dummy)
224#endif
225         endif
226c
227         if (abs(cfac(3)).gt.EPS) then
228            facl = cfac(3)
229            lfacl = lcfac(3)
230            nlfacl = nlcfac(3)
231#ifndef SECOND_DERIV
232            call xc_p81(tol_rho,facl,lfacl,nlfacl,rho_t,dnepscdn,
233     &         lnq,ipol,nepsc,lqwght,.false.,dummy)
234#else
235            call xc_p81_d2(tol_rho,facl,lfacl,nlfacl,rho_t,
236     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
237     &         .false.,dummy)
238#endif
239         endif
240c
241         if (abs(cfac(6)).gt.EPS.or.lfac) then
242            facl = cfac(6)
243            if (abs(facl).lt.EPS) facl = 1.0d0
244            lfacl = lcfac(6)
245            nlfacl = nlcfac(6)
246#ifndef SECOND_DERIV
247            call xc_pw91lda(tol_rho,facl,lfacl,nlfacl,rho_t,
248     &         dnepscdn,lnq,ipol,nepsc,lqwght,
249     &         .false.,dummy)
250#else
251            call xc_pw91lda_d2(tol_rho,facl,lfacl,nlfacl,rho_t,
252     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
253     &         .false.,dummy)
254#endif
255         endif
256c
257         if (abs(cfac(7)).gt.EPS) then
258            facl = cfac(7)
259#ifndef SECOND_DERIV
260            call xc_vwn_1_rpa(tol_rho,facl,rho_t,
261     &         dnepscdn,lnq,ipol,nepsc,lqwght,
262     &         .false.,dummy)
263#else
264            call xc_vwn_1_rpa_d2(tol_rho,facl,rho_t,
265     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
266     &         .false.,dummy)
267#endif
268         endif
269c
270         if (abs(cfac(8)).gt.EPS) then
271            facl = cfac(8)
272#ifndef SECOND_DERIV
273            call xc_vwn_1(tol_rho,facl,rho_t,
274     &         dnepscdn,lnq,ipol,nepsc,lqwght,
275     &         .false.,dummy)
276#else
277            call xc_vwn_1_d2(tol_rho,facl,rho_t,
278     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
279     &         .false.,dummy)
280#endif
281         endif
282c
283         if (abs(cfac(9)).gt.EPS) then
284            facl = cfac(9)
285#ifndef SECOND_DERIV
286            call xc_vwn_2(tol_rho,facl,rho_t,
287     &         dnepscdn,lnq,ipol,nepsc,lqwght,
288     &         .false.,dummy)
289#else
290            call xc_vwn_2_d2(tol_rho,facl,rho_t,
291     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
292     &         .false.,dummy)
293#endif
294         endif
295c
296         if (abs(cfac(10)).gt.EPS) then
297            facl = cfac(10)
298#ifndef SECOND_DERIV
299            call xc_vwn_3(tol_rho,facl,rho_t,
300     &         dnepscdn,lnq,ipol,nepsc,lqwght,
301     &         .false.,dummy)
302#else
303            call xc_vwn_3_d2(tol_rho,facl,rho_t,
304     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
305     &         .false.,dummy)
306#endif
307         endif
308c
309         if (abs(cfac(11)).gt.EPS) then
310            facl = cfac(11)
311#ifndef SECOND_DERIV
312            call xc_vwn_4(tol_rho,facl,rho_t,
313     &         dnepscdn,lnq,ipol,nepsc,lqwght,
314     &         .false.,dummy)
315#else
316            call xc_vwn_4_d2(tol_rho,facl,rho_t,
317     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
318     &         .false.,dummy)
319#endif
320         endif
321c        ==================
322c        PBE non-local part
323c        ==================
324         if(abs(nepsc).lt.tol_rho*tol_rho) goto 20
325c
326c        epsilon_c = n*epsilon_c / n
327c
328         epsc = nepsc/rhoval
329         if (ipol.eq.1) then
330            depscdna = dnepscdn(1)/rhoval-nepsc/(rhoval**2)
331            depscdnb = depscdna
332         else
333            depscdna = dnepscdn(1)/rhoval-nepsc/(rhoval**2)
334            depscdnb = dnepscdn(2)/rhoval-nepsc/(rhoval**2)
335         endif
336#ifdef SECOND_DERIV
337         if (ipol.eq.1) then
338            d2epscdna2   = d2nepscdn2(D2_RA_RA)/rhoval
339     &                     -dnepscdn(1)/(rhoval**2)
340     &                     -dnepscdn(1)/(rhoval**2)
341     &                     +2.0d0*nepsc/(rhoval**3)
342            d2epscdnadnb = d2nepscdn2(D2_RA_RB)/rhoval
343     &                     -dnepscdn(1)/(rhoval**2)
344     &                     -dnepscdn(1)/(rhoval**2)
345     &                     +2.0d0*nepsc/(rhoval**3)
346            d2epscdnb2   = d2epscdna2
347         else
348            d2epscdna2   = d2nepscdn2(D2_RA_RA)/rhoval
349     &                     -dnepscdn(1)/(rhoval**2)
350     &                     -dnepscdn(1)/(rhoval**2)
351     &                     +2.0d0*nepsc/(rhoval**3)
352            d2epscdnadnb = d2nepscdn2(D2_RA_RB)/rhoval
353     &                     -dnepscdn(1)/(rhoval**2)
354     &                     -dnepscdn(2)/(rhoval**2)
355     &                     +2.0d0*nepsc/(rhoval**3)
356            d2epscdnb2   = d2nepscdn2(D2_RB_RB)/rhoval
357     &                     -dnepscdn(2)/(rhoval**2)
358     &                     -dnepscdn(2)/(rhoval**2)
359     &                     +2.0d0*nepsc/(rhoval**3)
360         endif
361#endif
362c
363c        ks = 2*(3*PI**2*n)**(1/6)/sqrt(PI) and its derivs
364c
365         ks = 2.0d0*(3.0d0*PI*PI*rhoval)**(1.0d0/6.0d0)/dsqrt(PI)
366         dksdna = (1.0d0/6.0d0)*ks/rhoval
367         dksdnb = dksdna
368#ifdef SECOND_DERIV
369         d2ksdna2   = (1.0d0/6.0d0)*dksdna/rhoval
370     &              - (1.0d0/6.0d0)*ks/(rhoval**2)
371         d2ksdnadnb = d2ksdna2
372         d2ksdnb2   = d2ksdna2
373#endif
374c
375c        zeta = (na-nb)/n and its derivs
376c
377         if (ipol.eq.1) then
378            zeta = 0.0d0
379         else
380            zeta = (rho(n,2)-rho(n,3))/rhoval
381         endif
382         if(zeta.lt.-1.0d0) zeta=-1.0d0
383         if(zeta.gt. 1.0d0) zeta= 1.0d0
384         if (ipol.eq.1) then
385            dzetadna = 1.0d0/rhoval
386            dzetadnb = -1.0d0/rhoval
387#ifdef SECOND_DERIV
388            d2zetadna2   = -2.0d0/(rhoval**2)
389            d2zetadnadnb = 0.0d0
390            d2zetadnb2   = -2.0d0/(rhoval**2)
391#endif
392         else
393            dzetadna =  2.0d0*rho(n,3)/(rhoval**2)
394            dzetadnb = -2.0d0*rho(n,2)/(rhoval**2)
395#ifdef SECOND_DERIV
396            d2zetadna2   = -4.0d0*rho(n,3)/(rhoval**3)
397            d2zetadnadnb = 2.0d0*(rho(n,2)-rho(n,3))/(rhoval**3)
398            d2zetadnb2   = -4.0d0*rho(n,2)/(rhoval**3)
399#endif
400         endif
401c
402c        phi = (1/2)[(1+zeta)**(2/3)+(1-zeta)**(2/3)] and its derivs
403c
404         phi = 0.5d0*((1.0d0+zeta)**(2.0d0/3.0d0)
405     &               +(1.0d0-zeta)**(2.0d0/3.0d0))
406         if ((1.0d0-zeta).lt.tol_rho) then
407            dphidzeta = 0.5d0*(2.0d0/3.0d0)*(
408     &             (1.0d0+zeta)**(2.0d0/3.0d0)/(1.0d0+zeta))
409         else if ((1.0d0+zeta).lt.tol_rho) then
410            dphidzeta = 0.5d0*(2.0d0/3.0d0)*(
411     &            -(1.0d0-zeta)**(2.0d0/3.0d0)/(1.0d0-zeta))
412         else
413            dphidzeta = 0.5d0*(2.0d0/3.0d0)*(
414     &         (1.0d0+zeta)**(2.0d0/3.0d0)/(1.0d0+zeta)
415     &        -(1.0d0-zeta)**(2.0d0/3.0d0)/(1.0d0-zeta))
416         endif
417         dphidna = dphidzeta*dzetadna
418         dphidnb = dphidzeta*dzetadnb
419#ifdef SECOND_DERIV
420         if ((1.0d0-zeta).lt.tol_rho) then
421            d2phidzeta2 = -(1.0d0/9.0d0)*(
422     &         (1.0d0+zeta)**(2.0d0/3.0d0)/((1.0d0+zeta)**2))
423         else if ((1.0d0+zeta).lt.tol_rho) then
424            d2phidzeta2 = -(1.0d0/9.0d0)*(
425     &         (1.0d0-zeta)**(2.0d0/3.0d0)/((1.0d0-zeta)**2))
426         else
427            d2phidzeta2 = -(1.0d0/9.0d0)*(
428     &         (1.0d0+zeta)**(2.0d0/3.0d0)/((1.0d0+zeta)**2)
429     &        +(1.0d0-zeta)**(2.0d0/3.0d0)/((1.0d0-zeta)**2))
430         endif
431         d2phidna2   = d2phidzeta2*dzetadna*dzetadna
432     &               + dphidzeta*d2zetadna2
433         d2phidnadnb = d2phidzeta2*dzetadna*dzetadnb
434     &               + dphidzeta*d2zetadnadnb
435         d2phidnb2   = d2phidzeta2*dzetadnb*dzetadnb
436     &               + dphidzeta*d2zetadnb2
437#endif
438c
439c        t = |Nabla n|/(2*phi*ks*n) and its derivs
440c
441         t = dsqgamma/(2.0d0*phi*ks*rhoval)
442         dtdna = -t/rhoval-t/phi*dphidna-t/ks*dksdna
443         dtdnb = -t/rhoval-t/phi*dphidnb-t/ks*dksdnb
444#ifdef SECOND_DERIV
445         d2tdna2 = - dtdna/rhoval
446     &           + t/(rhoval**2)
447     &           - dtdna/phi*dphidna
448     &           + t/(phi**2)*(dphidna**2)
449     &           - t/phi*d2phidna2
450     &           - dtdna/ks*dksdna
451     &           + t/(ks**2)*(dksdna**2)
452     &           - t/ks*d2ksdna2
453         d2tdnadnb = - dtdnb/rhoval
454     &           + t/(rhoval**2)
455     &           - dtdnb/phi*dphidna
456     &           + t/(phi**2)*(dphidna*dphidnb)
457     &           - t/phi*d2phidnadnb
458     &           - dtdnb/ks*dksdna
459     &           + t/(ks**2)*(dksdna*dksdnb)
460     &           - t/ks*d2ksdnadnb
461         d2tdnb2 = - dtdna/rhoval
462     &           + t/(rhoval**2)
463     &           - dtdnb/phi*dphidnb
464     &           + t/(phi**2)*(dphidnb**2)
465     &           - t/phi*d2phidnb2
466     &           - dtdnb/ks*dksdnb
467     &           + t/(ks**2)*(dksdnb**2)
468     &           - t/ks*d2ksdnb2
469#endif
470c
471c        { ... } in A (see below) and its derivs
472c
473         argexp = -epsc/GAMMA/(phi**3)
474         dargexpdna = -depscdna/GAMMA/(phi**3)
475     &                +3.0d0*epsc/GAMMA/(phi**4)*dphidna
476         dargexpdnb = -depscdnb/GAMMA/(phi**3)
477     &                +3.0d0*epsc/GAMMA/(phi**4)*dphidnb
478#ifdef SECOND_DERIV
479         d2argexpdna2 = -d2epscdna2/GAMMA/(phi**3)
480     &        +3.0d0*depscdna/GAMMA/(phi**4)*dphidna
481     &        +3.0d0*depscdna/GAMMA/(phi**4)*dphidna
482     &        -12.0d0*epsc/GAMMA/(phi**5)*dphidna**2
483     &        +3.0d0*epsc/GAMMA/(phi**4)*d2phidna2
484         d2argexpdnadnb = -d2epscdnadnb/GAMMA/(phi**3)
485     &        +3.0d0*depscdna/GAMMA/(phi**4)*dphidnb
486     &        +3.0d0*depscdnb/GAMMA/(phi**4)*dphidna
487     &        -12.0d0*epsc/GAMMA/(phi**5)*dphidna*dphidnb
488     &        +3.0d0*epsc/GAMMA/(phi**4)*d2phidnadnb
489         d2argexpdnb2 = -d2epscdnb2/GAMMA/(phi**3)
490     &        +3.0d0*depscdnb/GAMMA/(phi**4)*dphidnb
491     &        +3.0d0*depscdnb/GAMMA/(phi**4)*dphidnb
492     &        -12.0d0*epsc/GAMMA/(phi**5)*dphidnb**2
493     &        +3.0d0*epsc/GAMMA/(phi**4)*d2phidnb2
494#endif
495c
496c        A = BETA/GAMMA [exp{-epsilon_c/(GAMMA*phi**3)}-1]**(-1)
497c
498         if (dabs(argexp).lt.EXPTOL) then
499            expinA=dexp(argexp)
500         else
501            expinA=0.0d0
502         endif
503         A = BETA/GAMMA/(expinA-1.0d0)
504         dAdna = -BETA/GAMMA*dargexpdna*expinA/(expinA-1.0d0)**2
505         dAdnb = -BETA/GAMMA*dargexpdnb*expinA/(expinA-1.0d0)**2
506#ifdef SECOND_DERIV
507         d2Adna2   = -BETA/GAMMA*d2argexpdna2
508     &               *expinA/(expinA-1.0d0)**2
509     &             - BETA/GAMMA*dargexpdna
510     &               *dargexpdna*expinA/(expinA-1.0d0)**2
511     &             + 2.0d0*BETA/GAMMA*dargexpdna*dargexpdna
512     &               *expinA*expinA/(expinA-1.0d0)**3
513         d2Adnadnb  = -BETA/GAMMA*d2argexpdnadnb
514     &               *expinA/(expinA-1.0d0)**2
515     &             - BETA/GAMMA*dargexpdna
516     &               *dargexpdnb*expinA/(expinA-1.0d0)**2
517     &             + 2.0d0*BETA/GAMMA*dargexpdna*dargexpdnb
518     &               *expinA*expinA/(expinA-1.0d0)**3
519         d2Adnb2   = -BETA/GAMMA*d2argexpdnb2
520     &               *expinA/(expinA-1.0d0)**2
521     &             - BETA/GAMMA*dargexpdnb
522     &               *dargexpdnb*expinA/(expinA-1.0d0)**2
523     &             + 2.0d0*BETA/GAMMA*dargexpdnb*dargexpdnb
524     &               *expinA*expinA/(expinA-1.0d0)**3
525#endif
526c
527c        fAt = (1 + A * t**2)/(1 + A * t**2 + A**2 * t**4) and its derivs
528c
529         fAtnum = 1.0d0
530         fAtden = 1.0d0+A*t**2
531         fAt = fAtnum/fAtden
532         dfAtdendt = 2.0d0*A*t
533         dfAtdendA = t**2
534         dfAtdt = (-fAtnum*dfAtdendt)/(fAtden**2)
535         dfAtdA = (-fAtnum*dfAtdendA)/(fAtden**2)
536         dfAtdna = dfAtdt * dtdna + dfAtdA * dAdna
537         dfAtdnb = dfAtdt * dtdnb + dfAtdA * dAdnb
538#ifdef SECOND_DERIV
539         d2fAtdendt2 = 2.0d0*A
540         d2fAtdendtdA = 2.0d0*t
541         d2fAtdendA2 = 0.0d0
542         d2fAtdt2  = (-fAtnum*d2fAtdendt2)
543     &               /(fAtden**2)
544     &               -2.0d0*(-fAtnum*dfAtdendt)
545     &               /(fAtden**3)*dfAtdendt
546         d2fAtdtdA = (-fAtnum*d2fAtdendtdA)
547     &               /(fAtden**2)
548     &               -2.0d0*(-fAtnum*dfAtdendt)
549     &               /(fAtden**3)*dfAtdendA
550         d2fAtdA2  = (-fAtnum*d2fAtdendA2)
551     &               /(fAtden**2)
552     &               -2.0d0*(-fAtnum*dfAtdendA)
553     &               /(fAtden**3)*dfAtdendA
554         d2fAtdna2 = d2fAtdt2*dtdna*dtdna + d2fAtdtdA*dtdna*dAdna
555     &             + dfAtdt*d2tdna2 + d2fAtdtdA*dAdna*dtdna
556     &             + d2fAtdA2*dAdna*dAdna + dfAtdA*d2Adna2
557         d2fAtdnb2 = d2fAtdt2*dtdnb*dtdnb + d2fAtdtdA*dtdnb*dAdnb
558     &             + dfAtdt*d2tdnb2 + d2fAtdtdA*dAdnb*dtdnb
559     &             + d2fAtdA2*dAdnb*dAdnb + dfAtdA*d2Adnb2
560         d2fAtdnadnb = d2fAtdt2*dtdna*dtdnb + d2fAtdtdA*dtdna*dAdnb
561     &             + dfAtdt*d2tdnadnb + d2fAtdtdA*dAdna*dtdnb
562     &             + d2fAtdA2*dAdna*dAdnb + dfAtdA*d2Adnadnb
563#endif
564c
565c        arglog = 1 + BETA/GAMMA * t**2 * fAt and its derivs
566c
567         arglog = 1.0d0 + BETA/GAMMA*t**2*fAt
568         darglogdna = BETA/GAMMA*(2.0d0*t*dtdna*fAt
569     &                            +t*t*dfAtdna)
570         darglogdnb = BETA/GAMMA*(2.0d0*t*dtdnb*fAt
571     &                            +t*t*dfAtdnb)
572#ifdef SECOND_DERIV
573         d2arglogdna2 = BETA/GAMMA*(2.0d0*dtdna*dtdna*fAt
574     &                             +2.0d0*t*d2tdna2*fAt
575     &                             +2.0d0*t*dtdna*dfAtdna
576     &                             +2.0d0*t*dtdna*dfAtdna
577     &                             +t*t*d2fAtdna2)
578         d2arglogdnb2 = BETA/GAMMA*(2.0d0*dtdnb*dtdnb*fAt
579     &                             +2.0d0*t*d2tdnb2*fAt
580     &                             +2.0d0*t*dtdnb*dfAtdnb
581     &                             +2.0d0*t*dtdnb*dfAtdnb
582     &                             +t*t*d2fAtdnb2)
583         d2arglogdnadnb = BETA/GAMMA*(2.0d0*dtdna*dtdnb*fAt
584     &                             +2.0d0*t*d2tdnadnb*fAt
585     &                             +2.0d0*t*dtdna*dfAtdnb
586     &                             +2.0d0*t*dtdnb*dfAtdna
587     &                             +t*t*d2fAtdnadnb)
588#endif
589c
590c        H0 = GAMMA * phi**3 * log{arglog} and its derivs
591c
592         H0 = GAMMA*(phi**3)*dlog(arglog)
593         dH0dna = GAMMA*(3.0d0*(phi**2)*dphidna*dlog(arglog)
594     &                  +(phi**3)*darglogdna/arglog)
595         dH0dnb = GAMMA*(3.0d0*(phi**2)*dphidnb*dlog(arglog)
596     &                  +(phi**3)*darglogdnb/arglog)
597#ifdef SECOND_DERIV
598         d2H0dna2 = GAMMA*(6.0d0*phi*dphidna*dphidna*dlog(arglog)
599     &                +3.0d0*(phi**2)*d2phidna2*dlog(arglog)
600     &                +6.0d0*(phi**2)*dphidna*darglogdna/arglog
601     &                +(phi**3)*d2arglogdna2/arglog
602     &                -(phi**3)*darglogdna*darglogdna/arglog/arglog)
603         d2H0dnadnb = GAMMA*(6.0d0*phi*dphidna*dphidnb*dlog(arglog)
604     &                +3.0d0*(phi**2)*d2phidnadnb*dlog(arglog)
605     &                +3.0d0*(phi**2)*dphidna*darglogdnb/arglog
606     &                +3.0d0*(phi**2)*dphidnb*darglogdna/arglog
607     &                +(phi**3)*d2arglogdnadnb/arglog
608     &                -(phi**3)*darglogdna*darglogdnb/arglog/arglog)
609         d2H0dnb2 = GAMMA*(6.0d0*phi*dphidnb*dphidnb*dlog(arglog)
610     &                +3.0d0*(phi**2)*d2phidnb2*dlog(arglog)
611     &                +6.0d0*(phi**2)*dphidnb*darglogdnb/arglog
612     &                +(phi**3)*d2arglogdnb2/arglog
613     &                -(phi**3)*darglogdnb*darglogdnb/arglog/arglog)
614#endif
615c
616c        Now we update Ec, Amat, and Amat2
617c
618         if (lfac) then
619            Ec = Ec+nepsc*qwght(n)*fac
620            if (ldew) ffunc(n)=ffunc(n)+nepsc*fac
621         endif
622         if (nlfac) then
623            Ec = Ec+(H0*rhoval)*qwght(n)*fac
624            if (ldew) ffunc(n)=ffunc(n)+(H0*rhoval)*fac
625         endif
626         if (lfac) then
627            Amat(n,1) = Amat(n,1) + dnepscdn(1)*fac
628            if (ipol.eq.2) Amat(n,2) = Amat(n,2) + dnepscdn(2)*fac
629#ifdef SECOND_DERIV
630            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
631     &                        + d2nepscdn2(D2_RA_RA)*fac
632            Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB)
633     &                        + d2nepscdn2(D2_RA_RB)*fac
634            if (ipol.eq.2)
635     &      Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB)
636     &                        + d2nepscdn2(D2_RB_RB)*fac
637#endif
638         endif
639         if (nlfac)then
640            Amat(n,1) = Amat(n,1) +  (H0 +
641     &                  rhoval*dH0dna)*fac
642            if (ipol.eq.2) Amat(n,2) = Amat(n,2) +  (H0 +
643     &                  rhoval*dH0dnb)*fac
644#ifdef SECOND_DERIV
645            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
646     &                + (2.d0*dH0dna + rhoval*d2H0dna2)*fac
647            Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB)
648     &                + (dH0dna + dH0dnb + rhoval*d2H0dnadnb)*fac
649            if (ipol.eq.2)
650     &      Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB)
651     &                + (2.d0*dH0dnb + rhoval*d2H0dnb2)*fac
652#endif
653         endif
654c
655c        Now we go into gradient-correction parts
656c        Note that the functional depends on |Nabla n| through "t" only
657c
658         if (dsqgamma.gt.TOLL)then
659            dtdg = 0.25d0/(phi*ks*rhoval)/dsqgamma
660            dfAtdg = dfAtdt*dtdg
661            darglogdg = BETA/GAMMA*(2.0d0*t*dtdg*fAt+t*t*dfAtdg)
662            dH0dg = GAMMA*(phi**3)*darglogdg/arglog
663            if (ipol.eq.1) then
664               Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dH0dg*rhoval*fac
665               Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + dH0dg*rhoval*fac*2.0d0
666            else
667               Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dH0dg*rhoval*fac
668               Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + dH0dg*rhoval*fac*2.0d0
669               Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + dH0dg*rhoval*fac
670            endif
671#ifdef SECOND_DERIV
672            d2tdg2 = -0.125d0/(phi*ks*rhoval)/(dsqgamma**3)
673            d2tdnadg = -dtdg/rhoval-dtdg/phi*dphidna
674     &                 -dtdg/ks*dksdna
675            d2tdnbdg = -dtdg/rhoval-dtdg/phi*dphidnb
676     &                 -dtdg/ks*dksdnb
677            d2fAtdg2 = d2fAtdt2*(dtdg**2)+dfAtdt*d2tdg2
678            d2fAtdnadg = d2fAtdt2*dtdg*dtdna
679     &                  +d2fAtdtdA*dtdg*dAdna
680     &                  +dfAtdt*d2tdnadg
681            d2fAtdnbdg = d2fAtdt2*dtdg*dtdnb
682     &                  +d2fAtdtdA*dtdg*dAdnb
683     &                  +dfAtdt*d2tdnbdg
684            d2arglogdnadg = BETA/GAMMA*(2.0d0*dtdna*dtdg*fAt
685     &                                 +2.0d0*t*d2tdnadg*fAt
686     &                                 +2.0d0*t*dtdg*dfAtdna
687     &                                 +2.0d0*t*dtdna*dfAtdg
688     &                                 +t*t*d2fAtdnadg)
689            d2arglogdnbdg = BETA/GAMMA*(2.0d0*dtdnb*dtdg*fAt
690     &                                 +2.0d0*t*d2tdnbdg*fAt
691     &                                 +2.0d0*t*dtdg*dfAtdnb
692     &                                 +2.0d0*t*dtdnb*dfAtdg
693     &                                 +t*t*d2fAtdnbdg)
694            d2arglogdg2 = BETA/GAMMA*(2.0d0*dtdg*dtdg*fAt
695     &                               +2.0d0*t*d2tdg2*fAt
696     &                               +2.0d0*t*dtdg*dfAtdg
697     &                               +2.0d0*t*dtdg*dfAtdg
698     &                               +t*t*d2fAtdg2)
699            d2H0dnadg = GAMMA*3.0d0*phi**2*dphidna*darglogdg/arglog
700     &                + GAMMA*phi**3*d2arglogdnadg/arglog
701     &                - GAMMA*phi**3*darglogdg*darglogdna/arglog**2
702            d2H0dnbdg = GAMMA*3.0d0*phi**2*dphidnb*darglogdg/arglog
703     &                + GAMMA*phi**3*d2arglogdnbdg/arglog
704     &                - GAMMA*phi**3*darglogdg*darglogdnb/arglog**2
705            d2H0dg2 = GAMMA*phi**3*d2arglogdg2/arglog
706     &              - GAMMA*phi**3*darglogdg*darglogdg/arglog**2
707            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA)
708     &             + (dH0dg + d2H0dnadg*rhoval)*fac
709            Cmat2(n,D2_RA_GAB) = Cmat2(n,D2_RA_GAB)
710     &             + 2.0d0*(dH0dg + d2H0dnadg*rhoval)*fac
711            Cmat2(n,D2_RA_GBB) = Cmat2(n,D2_RA_GBB)
712     &             + (dH0dg + d2H0dnadg*rhoval)*fac
713            Cmat2(n,D2_RB_GAA) = Cmat2(n,D2_RB_GAA)
714     &             + (dH0dg + d2H0dnbdg*rhoval)*fac
715            Cmat2(n,D2_RB_GAB) = Cmat2(n,D2_RB_GAB)
716     &             + 2.0d0*(dH0dg + d2H0dnbdg*rhoval)*fac
717            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB)
718     &             + (dH0dg + d2H0dnbdg*rhoval)*fac
719            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA)
720     &             + d2H0dg2*rhoval*fac
721            Cmat2(n,D2_GAA_GAB) = Cmat2(n,D2_GAA_GAB)
722     &             + 2.0d0*d2H0dg2*rhoval*fac
723            Cmat2(n,D2_GAA_GBB) = Cmat2(n,D2_GAA_GBB)
724     &             + d2H0dg2*rhoval*fac
725            Cmat2(n,D2_GAB_GAB) = Cmat2(n,D2_GAB_GAB)
726     &             + 4.0d0*d2H0dg2*rhoval*fac
727            Cmat2(n,D2_GAB_GBB) = Cmat2(n,D2_GAB_GBB)
728     &             + 2.0d0*d2H0dg2*rhoval*fac
729            Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB)
730     &             + d2H0dg2*rhoval*fac
731#endif
732         endif
733   20 continue
734c
735      return
736      end
737c
738#ifndef SECOND_DERIV
739#define SECOND_DERIV
740c
741c     Compile source again for the 2nd derivative case
742c
743#include "xc_spbe96.F"
744#endif
745