1#ifndef SECOND_DERIV
2      Subroutine xc_perdew91(tol_rho, cfac, lcfac, nlcfac, rho, delrho,
3     &                    Amat, Cmat, nq, ipol, Ec, qwght,
4     &                    ldew,ffunc)
5#else
6      Subroutine xc_perdew91_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     Input and other parameters
17c
18      integer lnq
19      integer ipol, nq
20      double precision dummy(1)
21      double precision cfac(*)
22      logical lcfac(*), nlcfac(*), ldew
23      logical lfac, nlfac, lfacl, nlfacl
24      double precision ffunc(*)
25      double precision fac, facl
26      double precision lqwght
27      double precision tol_rho
28c
29c     Constants in PBE functional
30c
31      double precision ALPHA, BETA, NU, CCZERO, CX, PI
32      double precision cnoname, ca, cb, cc, cd
33      parameter (ALPHA = 0.09d0)
34      parameter (CCZERO = 0.004235d0)
35      parameter (CX = -0.001667d0)
36      parameter (PI = 3.1415926535897932385d0)
37      parameter (cnoname = 2.568d0)
38      parameter (ca = 23.266d0)
39      parameter (cb = 7.389d-3)
40      parameter (cc = 8.723d0)
41      parameter (cd = 0.472d0)
42c
43c     Threshold parameters
44c
45      double precision TOLL, EXPTOL
46      double precision EPS
47      parameter (TOLL = 1.0D-40, EXPTOL = 40.0d0)
48      parameter (EPS = 1.0e-8)
49c
50c     Correlation energy
51c
52      double precision Ec
53c
54c     Charge Density
55c
56      double precision rho(nq,ipol*(ipol+1)/2)
57      double precision rho_t(3)
58c
59c     Charge Density Gradient
60c
61      double precision delrho(nq,3,ipol)
62      double precision dsqgamma
63c
64c     Quadrature Weights
65c
66      double precision qwght(nq)
67c
68c     Sampling Matrices for the XC Potential
69c
70      double precision Amat(nq,ipol), Cmat(nq,*)
71#ifdef SECOND_DERIV
72c
73c     Sampling Matrices for the XC Kernel
74c
75      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
76#endif
77c
78c     Intermediate derivative results, etc.
79c
80      integer n
81      double precision rhoval, gammaval
82      double precision nepsc, dnepscdn(2)
83      double precision epsc, depscdna, depscdnb
84      double precision H0, dH0dna, dH0dnb, dH0dg
85      double precision phi, dphidna, dphidnb, dphidzeta
86      double precision zeta, dzetadna, dzetadnb
87      double precision arglog, darglogdna, darglogdnb, darglogdg
88      double precision fAt, dfAtdt, dfAtdA
89      double precision fAtnum, dfAtnumdt, dfAtnumdA
90      double precision fAtden, dfAtdendt, dfAtdendA
91      double precision dfAtdna, dfAtdnb, dfAtdg
92      double precision A, dAdna, dAdnb
93      double precision t, dtdna, dtdnb, dtdg
94      double precision ks, dksdna, dksdnb
95      double precision argexp, dargexpdna, dargexpdnb
96      double precision expinA
97      double precision rs, drsdna, drsdnb
98      double precision ccrs, dccrsdna, dccrsdnb
99      double precision cnum, dcnumdna, dcnumdnb
100      double precision cden, dcdendna, dcdendnb
101      double precision kf, dkfdna, dkfdnb
102      double precision H1argexp, dH1argexpdna, dH1argexpdnb
103      double precision dH1argexpdg
104      double precision H1prefac, dH1prefacdna, dH1prefacdnb
105      double precision dH1prefacdg
106      double precision expinH1
107      double precision H1, dH1dna, dH1dnb, dH1dg
108#ifdef SECOND_DERIV
109      double precision d2nepscdn2(NCOL_AMAT2)
110      double precision d2epscdna2, d2epscdnadnb, d2epscdnb2
111      double precision d2H0dna2, d2H0dnadnb, d2H0dnb2
112      double precision d2H0dnadg, d2H0dnbdg, d2H0dg2
113      double precision d2phidzeta2, d2phidna2, d2phidnadnb, d2phidnb2
114      double precision d2zetadna2, d2zetadnadnb, d2zetadnb2
115      double precision d2arglogdna2, d2arglogdnb2, d2arglogdnadnb
116      double precision d2arglogdnadg, d2arglogdnbdg, d2arglogdg2
117      double precision d2fAtdt2, d2fAtdA2, d2fAtdtdA, d2fAtdg2
118      double precision d2fAtnumdt2, d2fAtnumdtdA, d2fAtnumdA2
119      double precision d2fAtdendt2, d2fAtdendtdA, d2fAtdendA2
120      double precision d2fAtdna2, d2fAtdnb2, d2fAtdnadnb
121      double precision d2fAtdnadg, d2fAtdnbdg
122      double precision d2Adna2, d2Adnadnb, d2Adnb2
123      double precision d2tdna2, d2tdnb2, d2tdnadnb
124      double precision d2tdg2, d2tdnadg, d2tdnbdg
125      double precision d2ksdna2, d2ksdnb2, d2ksdnadnb
126      double precision d2argexpdna2, d2argexpdnb2, d2argexpdnadnb
127      double precision d2rsdna2, d2rsdnb2, d2rsdnadnb
128      double precision d2ccrsdna2, d2ccrsdnb2, d2ccrsdnadnb
129      double precision d2cnumdna2, d2cnumdnb2, d2cnumdnadnb
130      double precision d2cdendna2, d2cdendnb2, d2cdendnadnb
131      double precision d2kfdna2, d2kfdnb2, d2kfdnadnb
132      double precision d2H1argexpdna2, d2H1argexpdnadnb, d2H1argexpdnb2
133      double precision d2H1argexpdnadg, d2H1argexpdnbdg, d2H1argexpdg2
134      double precision d2H1prefacdna2, d2H1prefacdnadnb, d2H1prefacdnb2
135      double precision d2H1prefacdnadg, d2H1prefacdnbdg, d2H1prefacdg2
136      double precision d2H1dna2, d2H1dnb2, d2H1dnadnb
137      double precision d2H1dnadg, d2H1dnbdg, d2H1dg2
138#endif
139c
140c References:
141c [a] J. P. Perdew et al., Phys. Rev. B 46, 6671 (1992).
142c [b] M. Rasolt and D. J. W. Geldart, Phys. Rev. B 34, 1325 (1986).
143c
144c  E_c(PW91) = Int n (epsilon_c + H0 + H1) dxdydz
145c
146c  n*epsilon_c                <=== supplied by another subroutine
147c  d(n*epsilon_c)/d(na)       <=== supplied by another subroutine
148c  d2(n*epsilon_c)/d(na)d(na) <=== supplied by another subroutine
149c  d2(n*epsilon_c)/d(na)d(nb) <=== supplied by another subroutine
150c  d2(n*epsilon_c)/d(nb)d(nb) <=== supplied by another subroutine
151c
152c  H0 = BETA**2/2/ALPHA * phi**3 * log{ 1 + 2*ALPHA/BETA * t**2 * [ ... ]}
153c
154c  BETA = NU * CCZERO
155c
156c  NU = (16/PI)(3*PI**2)**(1/3)
157c
158c  phi = (1/2)[(1+zeta)**(2/3)+(1-zeta)**(2/3)]
159c
160c  zeta = (na - nb)/n
161c
162c  [ ... ] = (1 + A * t**2)/(1 + A * t**2 + A**2 * t**4)
163c
164c  A = 2*ALPHA/BETA [exp{-2*ALPHA*epsilon_c/(BETA**2*phi**3)}-1]**(-1)
165c
166c  t = |Nabla n|/(2*phi*ks*n)
167c
168c  ks = 2 * (3 * PI**2 * n)**(1/6) / sqrt(PI)
169c
170c  |Nabla n| = sqrt(g_aa + g_bb + 2*g_ab)
171c
172c  H1 = NU * [CC(rs) - CCZERO - 3/7*CX] * Phi**3 * t**2
173c       * exp[-100 * phi**4 * (ks**2/kf**2) * t**2]
174c
175c  CC(rs) = 0.001 * (cnoname + ca*rs + cb*rs**2)/(1+cc*rs+cd*rs**2+10*cb*rs**3)
176c         - CX
177c
178c  rs = (3/4 / PI / n)**(1/3)
179c
180c  Names of variables
181c
182c  E_c(PW91)                 : Ec
183c  n (alpha+beta density)    : rhoval
184c  na, nb                    : rho(*,2), rho(*,3)
185c  epsilon_c                 : epsc
186c  H0                        : H0
187c  n*epsilon_c               : nepsc
188c  phi (also called "g")     : phi
189c  zeta                      : zeta
190c  { ... }                   : arglog
191c  [ ... ]                   : fAt
192c  (1 + A * t**2)            : fAtnum
193c  (1 + A * t**2 + A**2 * t**4) : fAtden
194c  A                         : A
195c  t                         : t
196c  |Nabla n|                 : gammaval
197c  ks                        : ks
198c  {-epsilon_c ... }         : argexp
199c  g_aa, g_bb, g_ab          : g
200c  H1                        : H1
201c  CC(rs)                    : ccrs
202c  CCZERO                    : CCZERO
203c  CX                        : CX
204c  kf                        : kf
205c  rs                        : rs
206c
207c  Derivatives of these are named like d...dna, d2...dnadnb,
208c  d2...dna2, etc.
209c
210      fac = cfac(5)
211      lfac = lcfac(5)
212      nlfac = nlcfac(5)
213c
214      NU = (16.0d0/PI)*(3.0d0*PI**2)**(1.0d0/3.0d0)
215      BETA = NU * CCZERO
216c
217c     ======> BOTH SPIN-RESTRICETED AND UNRESTRICTED <======
218c
219      do 20 n = 1, nq
220c
221c        n and zeta = (na - nb)/n
222c
223         rhoval = rho(n,1)
224         rho_t(1) = rho(n,1)
225         if (ipol.eq.2) then
226            rho_t(2) = rho(n,2)
227            rho_t(3) = rho(n,3)
228         endif
229         if (rhoval.le.tol_rho) goto 20
230         if (ipol.eq.1) then
231            gammaval = delrho(n,1,1)*delrho(n,1,1) +
232     &                 delrho(n,2,1)*delrho(n,2,1) +
233     &                 delrho(n,3,1)*delrho(n,3,1)
234         else
235            gammaval = delrho(n,1,1)*delrho(n,1,1) +
236     &                 delrho(n,1,2)*delrho(n,1,2) +
237     &                 delrho(n,2,1)*delrho(n,2,1) +
238     &                 delrho(n,2,2)*delrho(n,2,2) +
239     &                 delrho(n,3,1)*delrho(n,3,1) +
240     &                 delrho(n,3,2)*delrho(n,3,2) +
241     &           2.d0*(delrho(n,1,1)*delrho(n,1,2) +
242     &                 delrho(n,2,1)*delrho(n,2,2) +
243     &                 delrho(n,3,1)*delrho(n,3,2))
244         endif
245         dsqgamma = max(dsqrt(gammaval),tol_rho)
246         nepsc = 0.0d0
247         dnepscdn(1) = 0.0d0
248         if (ipol.eq.2) dnepscdn(2) = 0.0d0
249#ifdef SECOND_DERIV
250         d2nepscdn2(D2_RA_RA)=0.0d0
251         d2nepscdn2(D2_RA_RB)=0.0d0
252         if (ipol.eq.2) d2nepscdn2(D2_RB_RB)=0.0d0
253#endif
254c
255c        call for LDA bit
256c
257         lnq = 1
258         lqwght = 1.0d0
259c
260         if (abs(cfac(1)).gt.EPS) then
261            facl = cfac(1)
262#ifndef SECOND_DERIV
263            call xc_vwn_5(tol_rho,facl,rho_t,
264     &         dnepscdn,lnq,ipol,nepsc,lqwght,.false.,dummy)
265#else
266            call xc_vwn_5_d2(tol_rho,facl,rho_t,
267     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
268     &         .false.,dummy)
269#endif
270         endif
271c
272         if (abs(cfac(3)).gt.EPS) then
273            facl = cfac(3)
274            lfacl = lcfac(3)
275            nlfacl = nlcfac(3)
276#ifndef SECOND_DERIV
277            call xc_p81(tol_rho,facl,lfacl,nlfacl,rho_t,dnepscdn,
278     &         lnq,ipol,nepsc,lqwght,.false.,dummy)
279#else
280            call xc_p81_d2(tol_rho,facl,lfacl,nlfacl,rho_t,
281     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
282     &         .false.,dummy)
283#endif
284         endif
285c
286         if (abs(cfac(6)).gt.EPS.or.lfac) then
287            facl = cfac(6)
288            if (abs(facl).lt.EPS) facl = 1.0d0
289            lfacl = lcfac(6)
290            nlfacl = nlcfac(6)
291#ifndef SECOND_DERIV
292            call xc_pw91lda(tol_rho,facl,lfacl,nlfacl,rho_t,
293     &         dnepscdn,lnq,ipol,nepsc,lqwght,
294     &         .false.,dummy)
295#else
296            call xc_pw91lda_d2(tol_rho,facl,lfacl,nlfacl,rho_t,
297     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
298     &         .false.,dummy)
299#endif
300         endif
301c
302         if (abs(cfac(7)).gt.EPS) then
303            facl = cfac(7)
304#ifndef SECOND_DERIV
305            call xc_vwn_1_rpa(tol_rho,facl,rho_t,
306     &         dnepscdn,lnq,ipol,nepsc,lqwght,
307     &         .false.,dummy)
308#else
309            call xc_vwn_1_rpa_d2(tol_rho,facl,rho_t,
310     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
311     &         .false.,dummy)
312#endif
313         endif
314c
315         if (abs(cfac(8)).gt.EPS) then
316            facl = cfac(8)
317#ifndef SECOND_DERIV
318            call xc_vwn_1(tol_rho,facl,rho_t,
319     &         dnepscdn,lnq,ipol,nepsc,lqwght,
320     &         .false.,dummy)
321#else
322            call xc_vwn_1_d2(tol_rho,facl,rho_t,
323     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
324     &         .false.,dummy)
325#endif
326         endif
327c
328         if (abs(cfac(9)).gt.EPS) then
329            facl = cfac(9)
330#ifndef SECOND_DERIV
331            call xc_vwn_2(tol_rho,facl,rho_t,
332     &         dnepscdn,lnq,ipol,nepsc,lqwght,
333     &         .false.,dummy)
334#else
335            call xc_vwn_2_d2(tol_rho,facl,rho_t,
336     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
337     &         .false.,dummy)
338#endif
339         endif
340c
341         if (abs(cfac(10)).gt.EPS) then
342            facl = cfac(10)
343#ifndef SECOND_DERIV
344            call xc_vwn_3(tol_rho,facl,rho_t,
345     &         dnepscdn,lnq,ipol,nepsc,lqwght,
346     &         .false.,dummy)
347#else
348            call xc_vwn_3_d2(tol_rho,facl,rho_t,
349     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
350     &         .false.,dummy)
351#endif
352         endif
353c
354         if (abs(cfac(11)).gt.EPS) then
355            facl = cfac(11)
356#ifndef SECOND_DERIV
357            call xc_vwn_4(tol_rho,facl,rho_t,
358     &         dnepscdn,lnq,ipol,nepsc,lqwght,
359     &         .false.,dummy)
360#else
361            call xc_vwn_4_d2(tol_rho,facl,rho_t,
362     &         dnepscdn,d2nepscdn2,lnq,ipol,nepsc,lqwght,
363     &         .false.,dummy)
364#endif
365         endif
366c        ============
367c        PW91 H0 part
368c        ============
369         if(abs(nepsc).lt.tol_rho*tol_rho) goto 20
370c
371c        epsilon_c = n*epsilon_c / n
372c
373         epsc = nepsc/rhoval
374         if (ipol.eq.1) then
375            depscdna = dnepscdn(1)/rhoval-nepsc/(rhoval**2)
376            depscdnb = depscdna
377         else
378            depscdna = dnepscdn(1)/rhoval-nepsc/(rhoval**2)
379            depscdnb = dnepscdn(2)/rhoval-nepsc/(rhoval**2)
380         endif
381#ifdef SECOND_DERIV
382         if (ipol.eq.1) then
383            d2epscdna2   = d2nepscdn2(D2_RA_RA)/rhoval
384     &                     -dnepscdn(1)/(rhoval**2)
385     &                     -dnepscdn(1)/(rhoval**2)
386     &                     +2.0d0*nepsc/(rhoval**3)
387            d2epscdnadnb = d2nepscdn2(D2_RA_RB)/rhoval
388     &                     -dnepscdn(1)/(rhoval**2)
389     &                     -dnepscdn(1)/(rhoval**2)
390     &                     +2.0d0*nepsc/(rhoval**3)
391            d2epscdnb2   = d2epscdna2
392         else
393            d2epscdna2   = d2nepscdn2(D2_RA_RA)/rhoval
394     &                     -dnepscdn(1)/(rhoval**2)
395     &                     -dnepscdn(1)/(rhoval**2)
396     &                     +2.0d0*nepsc/(rhoval**3)
397            d2epscdnadnb = d2nepscdn2(D2_RA_RB)/rhoval
398     &                     -dnepscdn(1)/(rhoval**2)
399     &                     -dnepscdn(2)/(rhoval**2)
400     &                     +2.0d0*nepsc/(rhoval**3)
401            d2epscdnb2   = d2nepscdn2(D2_RB_RB)/rhoval
402     &                     -dnepscdn(2)/(rhoval**2)
403     &                     -dnepscdn(2)/(rhoval**2)
404     &                     +2.0d0*nepsc/(rhoval**3)
405         endif
406#endif
407c
408c        ks = 2*(3*PI**2*n)**(1/6)/sqrt(PI) and its derivs
409c
410         ks = 2.0d0*(3.0d0*PI*PI*rhoval)**(1.0d0/6.0d0)/dsqrt(PI)
411         dksdna = (1.0d0/6.0d0)*ks/rhoval
412         dksdnb = dksdna
413#ifdef SECOND_DERIV
414         d2ksdna2   = (1.0d0/6.0d0)*dksdna/rhoval
415     &              - (1.0d0/6.0d0)*ks/(rhoval**2)
416         d2ksdnadnb = d2ksdna2
417         d2ksdnb2   = d2ksdna2
418#endif
419c
420c        zeta = (na-nb)/n and its derivs
421c
422         if (ipol.eq.1) then
423            zeta = 0.0d0
424         else
425            zeta = (rho(n,2)-rho(n,3))/rhoval
426         endif
427         if(zeta.lt.-1.0d0) zeta=-1.0d0
428         if(zeta.gt. 1.0d0) zeta= 1.0d0
429         if (ipol.eq.1) then
430            dzetadna = 1.0d0/rhoval
431            dzetadnb = -1.0d0/rhoval
432#ifdef SECOND_DERIV
433            d2zetadna2   = -2.0d0/(rhoval**2)
434            d2zetadnadnb = 0.0d0
435            d2zetadnb2   = -2.0d0/(rhoval**2)
436#endif
437         else
438            dzetadna =  2.0d0*rho(n,3)/(rhoval**2)
439            dzetadnb = -2.0d0*rho(n,2)/(rhoval**2)
440#ifdef SECOND_DERIV
441            d2zetadna2   = -4.0d0*rho(n,3)/(rhoval**3)
442            d2zetadnadnb = 2.0d0*(rho(n,2)-rho(n,3))/(rhoval**3)
443            d2zetadnb2   = -4.0d0*rho(n,2)/(rhoval**3)
444#endif
445         endif
446c
447c        phi = (1/2)[(1+zeta)**(2/3)+(1-zeta)**(2/3)] and its derivs
448c
449         phi = 0.5d0*((1.0d0+zeta)**(2.0d0/3.0d0)
450     &               +(1.0d0-zeta)**(2.0d0/3.0d0))
451         if ((1.0d0-zeta).lt.tol_rho) then
452            dphidzeta = 0.5d0*(2.0d0/3.0d0)*(
453     &             (1.0d0+zeta)**(2.0d0/3.0d0)/(1.0d0+zeta))
454         else if ((1.0d0+zeta).lt.tol_rho) then
455            dphidzeta = 0.5d0*(2.0d0/3.0d0)*(
456     &            -(1.0d0-zeta)**(2.0d0/3.0d0)/(1.0d0-zeta))
457         else
458            dphidzeta = 0.5d0*(2.0d0/3.0d0)*(
459     &         (1.0d0+zeta)**(2.0d0/3.0d0)/(1.0d0+zeta)
460     &        -(1.0d0-zeta)**(2.0d0/3.0d0)/(1.0d0-zeta))
461         endif
462         dphidna = dphidzeta*dzetadna
463         dphidnb = dphidzeta*dzetadnb
464#ifdef SECOND_DERIV
465         if ((1.0d0-zeta).lt.tol_rho) then
466            d2phidzeta2 = -(1.0d0/9.0d0)*(
467     &         (1.0d0+zeta)**(2.0d0/3.0d0)/((1.0d0+zeta)**2))
468         else if ((1.0d0+zeta).lt.tol_rho) then
469            d2phidzeta2 = -(1.0d0/9.0d0)*(
470     &         (1.0d0-zeta)**(2.0d0/3.0d0)/((1.0d0-zeta)**2))
471         else
472            d2phidzeta2 = -(1.0d0/9.0d0)*(
473     &         (1.0d0+zeta)**(2.0d0/3.0d0)/((1.0d0+zeta)**2)
474     &        +(1.0d0-zeta)**(2.0d0/3.0d0)/((1.0d0-zeta)**2))
475         endif
476         d2phidna2   = d2phidzeta2*dzetadna*dzetadna
477     &               + dphidzeta*d2zetadna2
478         d2phidnadnb = d2phidzeta2*dzetadna*dzetadnb
479     &               + dphidzeta*d2zetadnadnb
480         d2phidnb2   = d2phidzeta2*dzetadnb*dzetadnb
481     &               + dphidzeta*d2zetadnb2
482#endif
483c
484c        t = |Nabla n|/(2*phi*ks*n) and its derivs
485c
486         t = dsqgamma/(2.0d0*phi*ks*rhoval)
487         dtdna = -t/rhoval-t/phi*dphidna-t/ks*dksdna
488         dtdnb = -t/rhoval-t/phi*dphidnb-t/ks*dksdnb
489#ifdef SECOND_DERIV
490         d2tdna2 = - dtdna/rhoval
491     &           + t/(rhoval**2)
492     &           - dtdna/phi*dphidna
493     &           + t/(phi**2)*(dphidna**2)
494     &           - t/phi*d2phidna2
495     &           - dtdna/ks*dksdna
496     &           + t/(ks**2)*(dksdna**2)
497     &           - t/ks*d2ksdna2
498         d2tdnadnb = - dtdnb/rhoval
499     &           + t/(rhoval**2)
500     &           - dtdnb/phi*dphidna
501     &           + t/(phi**2)*(dphidna*dphidnb)
502     &           - t/phi*d2phidnadnb
503     &           - dtdnb/ks*dksdna
504     &           + t/(ks**2)*(dksdna*dksdnb)
505     &           - t/ks*d2ksdnadnb
506         d2tdnb2 = - dtdna/rhoval
507     &           + t/(rhoval**2)
508     &           - dtdnb/phi*dphidnb
509     &           + t/(phi**2)*(dphidnb**2)
510     &           - t/phi*d2phidnb2
511     &           - dtdnb/ks*dksdnb
512     &           + t/(ks**2)*(dksdnb**2)
513     &           - t/ks*d2ksdnb2
514#endif
515c
516c        { ... } in A (see below) and its derivs
517c
518         argexp = -2.0d0*ALPHA*epsc/BETA**2/(phi**3)
519         dargexpdna = -2.0d0*ALPHA*depscdna/BETA**2/(phi**3)
520     &                +6.0d0*ALPHA*epsc/BETA**2/(phi**4)*dphidna
521         dargexpdnb = -2.0d0*ALPHA*depscdnb/BETA**2/(phi**3)
522     &                +6.0d0*ALPHA*epsc/BETA**2/(phi**4)*dphidnb
523#ifdef SECOND_DERIV
524         d2argexpdna2 = -2.0d0*ALPHA*d2epscdna2/BETA**2/(phi**3)
525     &        +6.0d0*ALPHA*depscdna/BETA**2/(phi**4)*dphidna
526     &        +6.0d0*ALPHA*depscdna/BETA**2/(phi**4)*dphidna
527     &        -24.0d0*ALPHA*epsc/BETA**2/(phi**5)*dphidna**2
528     &        +6.0d0*ALPHA*epsc/BETA**2/(phi**4)*d2phidna2
529         d2argexpdnadnb = -2.0d0*ALPHA*d2epscdnadnb/BETA**2/(phi**3)
530     &        +6.0d0*ALPHA*depscdna/BETA**2/(phi**4)*dphidnb
531     &        +6.0d0*ALPHA*depscdnb/BETA**2/(phi**4)*dphidna
532     &        -24.0d0*ALPHA*epsc/BETA**2/(phi**5)*dphidna*dphidnb
533     &        +6.0d0*ALPHA*epsc/BETA**2/(phi**4)*d2phidnadnb
534         d2argexpdnb2 = -2.0d0*ALPHA*d2epscdnb2/BETA**2/(phi**3)
535     &        +6.0d0*ALPHA*depscdnb/BETA**2/(phi**4)*dphidnb
536     &        +6.0d0*ALPHA*depscdnb/BETA**2/(phi**4)*dphidnb
537     &        -24.0d0*ALPHA*epsc/BETA**2/(phi**5)*dphidnb**2
538     &        +6.0d0*ALPHA*epsc/BETA**2/(phi**4)*d2phidnb2
539#endif
540c
541c        A = 2*ALPHA/BETA [exp{-2*ALPHA*epsilon_c/(BETA**2*phi**3)}-1]**(-1)
542c
543         if (dabs(argexp).lt.EXPTOL) then
544            expinA=dexp(argexp)
545         else
546            expinA=0.0d0
547         endif
548         A = 2.0d0*ALPHA/BETA/(expinA-1.0d0)
549         dAdna = -2.0d0*ALPHA/BETA*dargexpdna*expinA/(expinA-1.0d0)**2
550         dAdnb = -2.0d0*ALPHA/BETA*dargexpdnb*expinA/(expinA-1.0d0)**2
551#ifdef SECOND_DERIV
552         d2Adna2   = -2.0d0*ALPHA/BETA*d2argexpdna2
553     &               *expinA/(expinA-1.0d0)**2
554     &             - 2.0d0*ALPHA/BETA*dargexpdna
555     &               *dargexpdna*expinA/(expinA-1.0d0)**2
556     &             + 4.0d0*ALPHA/BETA*dargexpdna*dargexpdna
557     &               *expinA*expinA/(expinA-1.0d0)**3
558         d2Adnadnb  = -2.0d0*ALPHA/BETA*d2argexpdnadnb
559     &               *expinA/(expinA-1.0d0)**2
560     &             - 2.0d0*ALPHA/BETA*dargexpdna
561     &               *dargexpdnb*expinA/(expinA-1.0d0)**2
562     &             + 4.0d0*ALPHA/BETA*dargexpdna*dargexpdnb
563     &               *expinA*expinA/(expinA-1.0d0)**3
564         d2Adnb2   = -2.0d0*ALPHA/BETA*d2argexpdnb2
565     &               *expinA/(expinA-1.0d0)**2
566     &             - 2.0d0*ALPHA/BETA*dargexpdnb
567     &               *dargexpdnb*expinA/(expinA-1.0d0)**2
568     &             + 4.0d0*ALPHA/BETA*dargexpdnb*dargexpdnb
569     &               *expinA*expinA/(expinA-1.0d0)**3
570#endif
571c
572c        fAt = (1 + A * t**2)/(1 + A * t**2 + A**2 * t**4) and its derivs
573c
574         fAtnum = 1.0d0+A*t**2
575         fAtden = 1.0d0+A*t**2+A**2*t**4
576         fAt = fAtnum/fAtden
577         dfAtnumdt = 2.0d0*A*t
578         dfAtnumdA = t**2
579         dfAtdendt = 2.0d0*A*t+4.0d0*A**2*t**3
580         dfAtdendA = t**2+2.0d0*A*t**4
581         dfAtdt = (dfAtnumdt*fAtden-fAtnum*dfAtdendt)/(fAtden**2)
582         dfAtdA = (dfAtnumdA*fAtden-fAtnum*dfAtdendA)/(fAtden**2)
583         dfAtdna = dfAtdt * dtdna + dfAtdA * dAdna
584         dfAtdnb = dfAtdt * dtdnb + dfAtdA * dAdnb
585#ifdef SECOND_DERIV
586         d2fAtnumdt2 = 2.0d0*A
587         d2fAtdendt2 = 2.0d0*A+12.0d0*A**2*t**2
588         d2fAtnumdtdA = 2.0d0*t
589         d2fAtdendtdA = 2.0d0*t+8.0d0*A*t**3
590         d2fAtnumdA2 = 0.0d0
591         d2fAtdendA2 = 2.0d0*t**4
592         d2fAtdt2  = (d2fAtnumdt2*fAtden-fAtnum*d2fAtdendt2)
593     &               /(fAtden**2)
594     &               -2.0d0*(dfAtnumdt*fAtden-fAtnum*dfAtdendt)
595     &               /(fAtden**3)*dfAtdendt
596         d2fAtdtdA = (d2fAtnumdtdA*fAtden+dfAtnumdt*dfAtdendA
597     &               -dfAtnumdA*dfAtdendt-fAtnum*d2fAtdendtdA)
598     &               /(fAtden**2)
599     &               -2.0d0*(dfAtnumdt*fAtden-fAtnum*dfAtdendt)
600     &               /(fAtden**3)*dfAtdendA
601         d2fAtdA2  = (d2fAtnumdA2*fAtden-fAtnum*d2fAtdendA2)
602     &               /(fAtden**2)
603     &               -2.0d0*(dfAtnumdA*fAtden-fAtnum*dfAtdendA)
604     &               /(fAtden**3)*dfAtdendA
605         d2fAtdna2 = d2fAtdt2*dtdna*dtdna + d2fAtdtdA*dtdna*dAdna
606     &             + dfAtdt*d2tdna2 + d2fAtdtdA*dAdna*dtdna
607     &             + d2fAtdA2*dAdna*dAdna + dfAtdA*d2Adna2
608         d2fAtdnb2 = d2fAtdt2*dtdnb*dtdnb + d2fAtdtdA*dtdnb*dAdnb
609     &             + dfAtdt*d2tdnb2 + d2fAtdtdA*dAdnb*dtdnb
610     &             + d2fAtdA2*dAdnb*dAdnb + dfAtdA*d2Adnb2
611         d2fAtdnadnb = d2fAtdt2*dtdna*dtdnb + d2fAtdtdA*dtdna*dAdnb
612     &             + dfAtdt*d2tdnadnb + d2fAtdtdA*dAdna*dtdnb
613     &             + d2fAtdA2*dAdna*dAdnb + dfAtdA*d2Adnadnb
614#endif
615c
616c        arglog = 1 + 2*ALPHA/BETA * t**2 * fAt and its derivs
617c
618         arglog = 1.0d0 + 2.0d0*ALPHA/BETA*t**2*fAt
619         darglogdna = 2.0d0*ALPHA/BETA*(2.0d0*t*dtdna*fAt
620     &                            +t*t*dfAtdna)
621         darglogdnb = 2.0d0*ALPHA/BETA*(2.0d0*t*dtdnb*fAt
622     &                            +t*t*dfAtdnb)
623#ifdef SECOND_DERIV
624         d2arglogdna2 = 2.0d0*ALPHA/BETA*(2.0d0*dtdna*dtdna*fAt
625     &                             +2.0d0*t*d2tdna2*fAt
626     &                             +2.0d0*t*dtdna*dfAtdna
627     &                             +2.0d0*t*dtdna*dfAtdna
628     &                             +t*t*d2fAtdna2)
629         d2arglogdnb2 = 2.0d0*ALPHA/BETA*(2.0d0*dtdnb*dtdnb*fAt
630     &                             +2.0d0*t*d2tdnb2*fAt
631     &                             +2.0d0*t*dtdnb*dfAtdnb
632     &                             +2.0d0*t*dtdnb*dfAtdnb
633     &                             +t*t*d2fAtdnb2)
634         d2arglogdnadnb = 2.0d0*ALPHA/BETA*(2.0d0*dtdna*dtdnb*fAt
635     &                             +2.0d0*t*d2tdnadnb*fAt
636     &                             +2.0d0*t*dtdna*dfAtdnb
637     &                             +2.0d0*t*dtdnb*dfAtdna
638     &                             +t*t*d2fAtdnadnb)
639#endif
640c
641c        H0 = BETA**2/2/ALPHA * phi**3 * log{arglog} and its derivs
642c
643         H0 = 0.5d0*BETA**2/ALPHA*(phi**3)*dlog(arglog)
644         dH0dna = 0.5d0*BETA**2/ALPHA*(3.0d0*(phi**2)
645     &                  *dphidna*dlog(arglog)
646     &                  +(phi**3)*darglogdna/arglog)
647         dH0dnb = 0.5d0*BETA**2/ALPHA*(3.0d0*(phi**2)
648     &                  *dphidnb*dlog(arglog)
649     &                  +(phi**3)*darglogdnb/arglog)
650#ifdef SECOND_DERIV
651         d2H0dna2 = 0.5d0*BETA**2/ALPHA
652     &                *(6.0d0*phi*dphidna*dphidna*dlog(arglog)
653     &                +3.0d0*(phi**2)*d2phidna2*dlog(arglog)
654     &                +6.0d0*(phi**2)*dphidna*darglogdna/arglog
655     &                +(phi**3)*d2arglogdna2/arglog
656     &                -(phi**3)*darglogdna*darglogdna/arglog/arglog)
657         d2H0dnadnb = 0.5d0*BETA**2/ALPHA
658     &                *(6.0d0*phi*dphidna*dphidnb*dlog(arglog)
659     &                +3.0d0*(phi**2)*d2phidnadnb*dlog(arglog)
660     &                +3.0d0*(phi**2)*dphidna*darglogdnb/arglog
661     &                +3.0d0*(phi**2)*dphidnb*darglogdna/arglog
662     &                +(phi**3)*d2arglogdnadnb/arglog
663     &                -(phi**3)*darglogdna*darglogdnb/arglog/arglog)
664         d2H0dnb2 = 0.5d0*BETA**2/ALPHA
665     &                *(6.0d0*phi*dphidnb*dphidnb*dlog(arglog)
666     &                +3.0d0*(phi**2)*d2phidnb2*dlog(arglog)
667     &                +6.0d0*(phi**2)*dphidnb*darglogdnb/arglog
668     &                +(phi**3)*d2arglogdnb2/arglog
669     &                -(phi**3)*darglogdnb*darglogdnb/arglog/arglog)
670#endif
671c        ============
672c        PW91 H1 part
673c        ============
674c
675c        rs = (3/4 / PI / n)**(1/3) and its derivs
676c
677         rs = (0.75d0 / PI / dabs(rhoval))**(1.0d0/3.0d0)
678         drsdna = (-1.0d0/3.0d0)*rs/rhoval
679         drsdnb = drsdna
680#ifdef SECOND_DERIV
681         d2rsdna2 = (-4.0d0/3.0d0)*drsdna/rhoval
682         d2rsdnb2 = d2rsdna2
683         d2rsdnadnb = d2rsdna2
684#endif
685c
686c        CC(rs) and its derivs
687c
688         cnum = cnoname + ca*rs + cb*rs**2
689         cden = 1.0d0 + cc*rs + cd*rs**2 + 10.0d0*cb*rs**3
690         ccrs = 1.0d-3*cnum/cden - CX
691         dcnumdna = (ca + 2.0d0*cb*rs)*drsdna
692         dcnumdnb = (ca + 2.0d0*cb*rs)*drsdnb
693         dcdendna = (cc + 2.0d0*cd*rs + 30.0d0*cb*rs**2)*drsdna
694         dcdendnb = (cc + 2.0d0*cd*rs + 30.0d0*cb*rs**2)*drsdnb
695         dccrsdna = 1.0d-3*(dcnumdna*cden-cnum*dcdendna)/cden**2
696         dccrsdnb = 1.0d-3*(dcnumdnb*cden-cnum*dcdendnb)/cden**2
697#ifdef SECOND_DERIV
698         d2cnumdna2 = (ca + 2.0d0*cb*rs)*d2rsdna2
699     &              + (2.0d0*cb)*drsdna*drsdna
700         d2cnumdnb2 = (ca + 2.0d0*cb*rs)*d2rsdnb2
701     &              + (2.0d0*cb)*drsdnb*drsdnb
702         d2cnumdnadnb = (ca + 2.0d0*cb*rs)*d2rsdnadnb
703     &              + (2.0d0*cb)*drsdna*drsdnb
704         d2cdendna2 = (cc + 2.0d0*cd*rs + 30.0d0*cb*rs**2)*d2rsdna2
705     &              + (2.0d0*cd+60.0d0*cb*rs)*drsdna*drsdna
706         d2cdendnb2 = (cc + 2.0d0*cd*rs + 30.0d0*cb*rs**2)*d2rsdnb2
707     &              + (2.0d0*cd+60.0d0*cb*rs)*drsdnb*drsdnb
708         d2cdendnadnb = (cc + 2.0d0*cd*rs + 30.0d0*cb*rs**2)*d2rsdnadnb
709     &              + (2.0d0*cd+60.0d0*cb*rs)*drsdna*drsdnb
710         d2ccrsdna2 = 1.0d-3*(d2cnumdna2*cden-cnum*d2cdendna2)/cden**2
711     &      - 1.0d-3*2.0d0*(dcnumdna*cden-cnum*dcdendna)
712     &        /cden**3*dcdendna
713         d2ccrsdnb2 = 1.0d-3*(d2cnumdnb2*cden-cnum*d2cdendnb2)/cden**2
714     &      - 1.0d-3*2.0d0*(dcnumdnb*cden-cnum*dcdendnb)
715     &        /cden**3*dcdendnb
716         d2ccrsdnadnb
717     &      = 1.0d-3*(d2cnumdnadnb*cden-cnum*d2cdendnadnb)/cden**2
718     &      + 1.0d-3*(dcnumdna*dcdendnb-dcnumdnb*dcdendna)/cden**2
719     &      - 1.0d-3*2.0d0*(dcnumdna*cden-cnum*dcdendna)
720     &        /cden**3*dcdendnb
721#endif
722c
723c        kf = (3 * PI**2 * n)**(1/3)
724c
725         kf = (3.0d0*PI**2*rhoval)**(1.0d0/3.0d0)
726         dkfdna = (1.0d0/3.0d0)*kf/rhoval
727         dkfdnb = dkfdna
728#ifdef SECOND_DERIV
729         d2kfdna2 = (-2.0d0/3.0d0)*dkfdna/rhoval
730         d2kfdnadnb = (-2.0d0/3.0d0)*dkfdna/rhoval
731         d2kfdnb2 = (-2.0d0/3.0d0)*dkfdnb/rhoval
732#endif
733c
734c        H1argexp = -100 * phi**4 * (ks**2/kf**2) * t**2 and its derivs
735c
736         H1argexp = -100.0d0*phi**4*(ks**2/kf**2)*t**2
737         dH1argexpdna = 4.0d0*H1argexp/phi*dphidna
738     &                 +2.0d0*H1argexp/ks*dksdna
739     &                 -2.0d0*H1argexp/kf*dkfdna
740     &                 +2.0d0*H1argexp/t*dtdna
741         dH1argexpdnb = 4.0d0*H1argexp/phi*dphidnb
742     &                 +2.0d0*H1argexp/ks*dksdnb
743     &                 -2.0d0*H1argexp/kf*dkfdnb
744     &                 +2.0d0*H1argexp/t*dtdnb
745#ifdef SECOND_DERIV
746         d2H1argexpdna2 = 4.0d0*dH1argexpdna/phi*dphidna
747     &                   -4.0d0*H1argexp/phi**2*dphidna**2
748     &                   +4.0d0*H1argexp/phi*d2phidna2
749     &                   +2.0d0*dH1argexpdna/ks*dksdna
750     &                   -2.0d0*H1argexp/ks**2*dksdna**2
751     &                   +2.0d0*H1argexp/ks*d2ksdna2
752     &                   -2.0d0*dH1argexpdna/kf*dkfdna
753     &                   +2.0d0*H1argexp/kf**2*dkfdna**2
754     &                   -2.0d0*H1argexp/kf*d2kfdna2
755     &                   +2.0d0*dH1argexpdna/t*dtdna
756     &                   -2.0d0*H1argexp/t**2*dtdna**2
757     &                   +2.0d0*H1argexp/t*d2tdna2
758         d2H1argexpdnb2 = 4.0d0*dH1argexpdnb/phi*dphidnb
759     &                   -4.0d0*H1argexp/phi**2*dphidnb**2
760     &                   +4.0d0*H1argexp/phi*d2phidnb2
761     &                   +2.0d0*dH1argexpdnb/ks*dksdnb
762     &                   -2.0d0*H1argexp/ks**2*dksdnb**2
763     &                   +2.0d0*H1argexp/ks*d2ksdnb2
764     &                   -2.0d0*dH1argexpdnb/kf*dkfdnb
765     &                   +2.0d0*H1argexp/kf**2*dkfdnb**2
766     &                   -2.0d0*H1argexp/kf*d2kfdnb2
767     &                   +2.0d0*dH1argexpdnb/t*dtdnb
768     &                   -2.0d0*H1argexp/t**2*dtdnb**2
769     &                   +2.0d0*H1argexp/t*d2tdnb2
770         d2H1argexpdnadnb = 4.0d0*dH1argexpdna/phi*dphidnb
771     &                   -4.0d0*H1argexp/phi**2*dphidna*dphidnb
772     &                   +4.0d0*H1argexp/phi*d2phidnadnb
773     &                   +2.0d0*dH1argexpdna/ks*dksdnb
774     &                   -2.0d0*H1argexp/ks**2*dksdna*dksdnb
775     &                   +2.0d0*H1argexp/ks*d2ksdnadnb
776     &                   -2.0d0*dH1argexpdna/kf*dkfdnb
777     &                   +2.0d0*H1argexp/kf**2*dkfdna*dkfdnb
778     &                   -2.0d0*H1argexp/kf*d2kfdnadnb
779     &                   +2.0d0*dH1argexpdna/t*dtdnb
780     &                   -2.0d0*H1argexp/t**2*dtdna*dtdnb
781     &                   +2.0d0*H1argexp/t*d2tdnadnb
782#endif
783c
784c        H1prefac = NU*[CC(rs) - CC(0)-3/7*CX] * g**3 * t**2
785c
786         H1prefac = NU*(ccrs - CCZERO - (3.0d0/7.0d0)*CX)
787     &            * phi**3 * t**2
788         dH1prefacdna = NU*dccrsdna * phi**3 * t**2
789     &                + 3.0d0*H1prefac/phi*dphidna
790     &                + 2.0d0*H1prefac/t*dtdna
791         dH1prefacdnb = NU*dccrsdnb * phi**3 * t**2
792     &                + 3.0d0*H1prefac/phi*dphidnb
793     &                + 2.0d0*H1prefac/t*dtdnb
794#ifdef SECOND_DERIV
795         d2H1prefacdna2 = NU*d2ccrsdna2 * phi**3 * t**2
796     &                  + 3.0d0* NU*dccrsdna * phi**2*dphidna * t**2
797     &                  + 2.0d0* NU*dccrsdna * phi**3 * t*dtdna
798     &                  + 3.0d0*dH1prefacdna/phi*dphidna
799     &                  - 3.0d0*H1prefac/phi**2*dphidna**2
800     &                  + 3.0d0*H1prefac/phi*d2phidna2
801     &                  + 2.0d0*dH1prefacdna/t*dtdna
802     &                  - 2.0d0*H1prefac/t**2*dtdna**2
803     &                  + 2.0d0*H1prefac/t*d2tdna2
804         d2H1prefacdnb2 = NU*d2ccrsdnb2 * phi**3 * t**2
805     &                  + 3.0d0* NU*dccrsdnb * phi**2*dphidnb * t**2
806     &                  + 2.0d0* NU*dccrsdnb * phi**3 * t*dtdnb
807     &                  + 3.0d0*dH1prefacdnb/phi*dphidnb
808     &                  - 3.0d0*H1prefac/phi**2*dphidnb**2
809     &                  + 3.0d0*H1prefac/phi*d2phidnb2
810     &                  + 2.0d0*dH1prefacdnb/t*dtdnb
811     &                  - 2.0d0*H1prefac/t**2*dtdnb**2
812     &                  + 2.0d0*H1prefac/t*d2tdnb2
813         d2H1prefacdnadnb = NU*d2ccrsdnadnb * phi**3 * t**2
814     &                  + 3.0d0* NU*dccrsdna * phi**2*dphidnb * t**2
815     &                  + 2.0d0* NU*dccrsdna * phi**3 * t*dtdnb
816     &                  + 3.0d0*dH1prefacdna/phi*dphidnb
817     &                  - 3.0d0*H1prefac/phi**2*dphidna*dphidnb
818     &                  + 3.0d0*H1prefac/phi*d2phidnadnb
819     &                  + 2.0d0*dH1prefacdna/t*dtdnb
820     &                  - 2.0d0*H1prefac/t**2*dtdna*dtdnb
821     &                  + 2.0d0*H1prefac/t*d2tdnadnb
822#endif
823c
824c        H1 = H1prefac * exp(H1argexp)
825c
826         if (dabs(H1argexp).lt.EXPTOL) then
827            expinH1 = dexp(H1argexp)
828         else
829            expinH1 = 0.0d0
830         endif
831         H1 = H1prefac * expinH1
832         dH1dna = dH1prefacdna * expinH1
833     &          + H1prefac * dH1argexpdna * expinH1
834         dH1dnb = dH1prefacdnb * expinH1
835     &          + H1prefac * dH1argexpdnb * expinH1
836#ifdef SECOND_DERIV
837         d2H1dna2 = d2H1prefacdna2 * expinH1
838     &            + dH1prefacdna * dH1argexpdna * expinH1
839     &            + dH1prefacdna * dH1argexpdna * expinH1
840     &            + H1prefac * d2H1argexpdna2 * expinH1
841     &            + H1prefac * dH1argexpdna * dH1argexpdna * expinH1
842         d2H1dnb2 = d2H1prefacdnb2 * expinH1
843     &            + dH1prefacdnb * dH1argexpdnb * expinH1
844     &            + dH1prefacdnb * dH1argexpdnb * expinH1
845     &            + H1prefac * d2H1argexpdnb2 * expinH1
846     &            + H1prefac * dH1argexpdnb * dH1argexpdnb * expinH1
847         d2H1dnadnb = d2H1prefacdnadnb * expinH1
848     &            + dH1prefacdna * dH1argexpdnb * expinH1
849     &            + dH1prefacdnb * dH1argexpdna * expinH1
850     &            + H1prefac * d2H1argexpdnadnb * expinH1
851     &            + H1prefac * dH1argexpdna * dH1argexpdnb * expinH1
852#endif
853c
854c        Now we update Ec, Amat, and Amat2
855c
856         if (lfac) then
857            Ec = Ec+nepsc*qwght(n)*fac
858            if (ldew) ffunc(n)=ffunc(n)+nepsc*fac
859         endif
860         if (nlfac) then
861            Ec = Ec+(H0*rhoval+H1*rhoval)*qwght(n)*fac
862            if (ldew) ffunc(n)=ffunc(n)+(H0*rhoval+H1*rhoval)*fac
863         endif
864         if (lfac) then
865            Amat(n,1) = Amat(n,1) + dnepscdn(1)*fac
866            if (ipol.eq.2) Amat(n,2) = Amat(n,2) + dnepscdn(2)*fac
867#ifdef SECOND_DERIV
868            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
869     &                        + d2nepscdn2(D2_RA_RA)*fac
870            Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB)
871     &                        + d2nepscdn2(D2_RA_RB)*fac
872            if (ipol.eq.2)
873     &      Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB)
874     &                        + d2nepscdn2(D2_RB_RB)*fac
875#endif
876         endif
877         if (nlfac)then
878            Amat(n,1) = Amat(n,1) + (H0 + H1 +
879     &         rhoval*dH0dna + rhoval*dH1dna)*fac
880            if (ipol.eq.2) Amat(n,2) = Amat(n,2) + (H0 + H1 +
881     &         rhoval*dH0dnb + rhoval*dH1dnb)*fac
882#ifdef SECOND_DERIV
883            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA)
884     &         + (2.0d0*dH0dna + rhoval*d2H0dna2
885     &         +  2.0d0*dH1dna + rhoval*d2H1dna2)*fac
886            Amat2(n,D2_RA_RB) = Amat2(n,D2_RA_RB)
887     &         + (dH0dna + dH0dnb + rhoval*d2H0dnadnb
888     &         +  dH1dna + dH1dnb + rhoval*d2H1dnadnb)*fac
889            if (ipol.eq.2)
890     &      Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB)
891     &         + (2.0d0*dH0dnb + rhoval*d2H0dnb2
892     &         +  2.0d0*dH1dnb + rhoval*d2H1dnb2)*fac
893#endif
894         endif
895c
896c        Now we go into gradient-correction parts
897c        Note that the functional depends on |Nabla n| through "t" only
898c
899cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
900c        goto 20
901cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
902         if (dsqgamma.gt.TOLL)then
903c
904c           H0 part
905c
906            dtdg = 0.25d0/(phi*ks*rhoval)/dsqgamma
907            dfAtdg = dfAtdt*dtdg
908            darglogdg = 2.0d0*ALPHA/BETA*(2.0d0*t*dtdg*fAt+t*t*dfAtdg)
909            dH0dg = 0.5d0*BETA**2/ALPHA*(phi**3)*darglogdg/arglog
910c
911c           H1 part
912c
913            dH1argexpdg = 2.0d0*H1argexp/t*dtdg
914            dH1prefacdg = 2.0d0*H1prefac/t*dtdg
915            dH1dg = dH1prefacdg * expinH1
916     &            + H1prefac * dH1argexpdg * expinH1
917c
918c           Now form Cmat
919c
920            if (ipol.eq.1) then
921               Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dH0dg*rhoval*fac
922     &                                         + dH1dg*rhoval*fac
923               Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + dH0dg*rhoval*fac*2.0d0
924     &                                         + dH1dg*rhoval*fac*2.0d0
925            else
926               Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + dH0dg*rhoval*fac
927     &                                         + dH1dg*rhoval*fac
928               Cmat(n,D1_GAB) = Cmat(n,D1_GAB) + dH0dg*rhoval*fac*2.0d0
929     &                                         + dH1dg*rhoval*fac*2.0d0
930               Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + dH0dg*rhoval*fac
931     &                                         + dH1dg*rhoval*fac
932            endif
933#ifdef SECOND_DERIV
934c
935c           H0 part
936c
937            d2tdg2 = -0.125d0/(phi*ks*rhoval)/(dsqgamma**3)
938            d2tdnadg = -dtdg/rhoval-dtdg/phi*dphidna
939     &                 -dtdg/ks*dksdna
940            d2tdnbdg = -dtdg/rhoval-dtdg/phi*dphidnb
941     &                 -dtdg/ks*dksdnb
942            d2fAtdg2 = d2fAtdt2*(dtdg**2)+dfAtdt*d2tdg2
943            d2fAtdnadg = d2fAtdt2*dtdg*dtdna
944     &                  +d2fAtdtdA*dtdg*dAdna
945     &                  +dfAtdt*d2tdnadg
946            d2fAtdnbdg = d2fAtdt2*dtdg*dtdnb
947     &                  +d2fAtdtdA*dtdg*dAdnb
948     &                  +dfAtdt*d2tdnbdg
949            d2arglogdnadg = 2.0d0*ALPHA/BETA*(2.0d0*dtdna*dtdg*fAt
950     &                                 +2.0d0*t*d2tdnadg*fAt
951     &                                 +2.0d0*t*dtdg*dfAtdna
952     &                                 +2.0d0*t*dtdna*dfAtdg
953     &                                 +t*t*d2fAtdnadg)
954            d2arglogdnbdg = 2.0d0*ALPHA/BETA*(2.0d0*dtdnb*dtdg*fAt
955     &                                 +2.0d0*t*d2tdnbdg*fAt
956     &                                 +2.0d0*t*dtdg*dfAtdnb
957     &                                 +2.0d0*t*dtdnb*dfAtdg
958     &                                 +t*t*d2fAtdnbdg)
959            d2arglogdg2 = 2.0d0*ALPHA/BETA*(2.0d0*dtdg*dtdg*fAt
960     &                               +2.0d0*t*d2tdg2*fAt
961     &                               +2.0d0*t*dtdg*dfAtdg
962     &                               +2.0d0*t*dtdg*dfAtdg
963     &                               +t*t*d2fAtdg2)
964            d2H0dnadg = 0.5d0*BETA**2/ALPHA*3.0d0*phi**2
965     &                  *dphidna*darglogdg/arglog
966     &                + 0.5d0*BETA**2/ALPHA*phi**3
967     &                  *d2arglogdnadg/arglog
968     &                - 0.5d0*BETA**2/ALPHA*phi**3
969     &                  *darglogdg*darglogdna/arglog**2
970            d2H0dnbdg = 0.5d0*BETA**2/ALPHA*3.0d0
971     &                  *phi**2*dphidnb*darglogdg/arglog
972     &                + 0.5d0*BETA**2/ALPHA*phi**3
973     &                  *d2arglogdnbdg/arglog
974     &                - 0.5d0*BETA**2/ALPHA*phi**3
975     &                  *darglogdg*darglogdnb/arglog**2
976            d2H0dg2 = 0.5d0*BETA**2/ALPHA*phi**3
977     &                *d2arglogdg2/arglog
978     &              - 0.5d0*BETA**2/ALPHA*phi**3
979     &                *darglogdg*darglogdg/arglog**2
980c
981c           H1 part
982c
983            d2H1argexpdnadg = 2.0d0*dH1argexpdna/t*dtdg
984     &                      - 2.0d0*H1argexp/t**2*dtdg*dtdna
985     &                      + 2.0d0*H1argexp/t*d2tdnadg
986            d2H1argexpdnbdg = 2.0d0*dH1argexpdnb/t*dtdg
987     &                      - 2.0d0*H1argexp/t**2*dtdg*dtdnb
988     &                      + 2.0d0*H1argexp/t*d2tdnbdg
989            d2H1argexpdg2 = 2.0d0*dH1argexpdg/t*dtdg
990     &                    - 2.0d0*H1argexp/t**2*dtdg*dtdg
991     &                    + 2.0d0*H1argexp/t*d2tdg2
992            d2H1prefacdnadg = 2.0d0*dH1prefacdna/t*dtdg
993     &                  - 2.0d0*H1prefac/t**2*dtdna*dtdg
994     &                  + 2.0d0*H1prefac/t*d2tdnadg
995            d2H1prefacdnbdg = 2.0d0*dH1prefacdnb/t*dtdg
996     &                  - 2.0d0*H1prefac/t**2*dtdnb*dtdg
997     &                  + 2.0d0*H1prefac/t*d2tdnbdg
998            d2H1prefacdg2 = 2.0d0*dH1prefacdg/t*dtdg
999     &                  - 2.0d0*H1prefac/t**2*dtdg*dtdg
1000     &                  + 2.0d0*H1prefac/t*d2tdg2
1001            d2H1dnadg = d2H1prefacdnadg * expinH1
1002     &            + dH1prefacdna * dH1argexpdg * expinH1
1003     &            + dH1prefacdg * dH1argexpdna * expinH1
1004     &            + H1prefac * d2H1argexpdnadg * expinH1
1005     &            + H1prefac * dH1argexpdna * dH1argexpdg * expinH1
1006            d2H1dnbdg = d2H1prefacdnbdg * expinH1
1007     &            + dH1prefacdnb * dH1argexpdg * expinH1
1008     &            + dH1prefacdg * dH1argexpdnb * expinH1
1009     &            + H1prefac * d2H1argexpdnbdg * expinH1
1010     &            + H1prefac * dH1argexpdnb * dH1argexpdg * expinH1
1011            d2H1dg2 = d2H1prefacdg2 * expinH1
1012     &            + dH1prefacdg * dH1argexpdg * expinH1
1013     &            + dH1prefacdg * dH1argexpdg * expinH1
1014     &            + H1prefac * d2H1argexpdg2 * expinH1
1015     &            + H1prefac * dH1argexpdg * dH1argexpdg * expinH1
1016c
1017c           Now form Cmat2
1018c
1019            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA)
1020     &             + (dH0dg + d2H0dnadg*rhoval)*fac
1021     &             + (dH1dg + d2H1dnadg*rhoval)*fac
1022            Cmat2(n,D2_RA_GAB) = Cmat2(n,D2_RA_GAB)
1023     &             + 2.0d0*(dH0dg + d2H0dnadg*rhoval)*fac
1024     &             + 2.0d0*(dH1dg + d2H1dnadg*rhoval)*fac
1025            Cmat2(n,D2_RA_GBB) = Cmat2(n,D2_RA_GBB)
1026     &             + (dH0dg + d2H0dnadg*rhoval)*fac
1027     &             + (dH1dg + d2H1dnadg*rhoval)*fac
1028            Cmat2(n,D2_RB_GAA) = Cmat2(n,D2_RB_GAA)
1029     &             + (dH0dg + d2H0dnbdg*rhoval)*fac
1030     &             + (dH1dg + d2H1dnbdg*rhoval)*fac
1031            Cmat2(n,D2_RB_GAB) = Cmat2(n,D2_RB_GAB)
1032     &             + 2.0d0*(dH0dg + d2H0dnbdg*rhoval)*fac
1033     &             + 2.0d0*(dH1dg + d2H1dnbdg*rhoval)*fac
1034            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB)
1035     &             + (dH0dg + d2H0dnbdg*rhoval)*fac
1036     &             + (dH1dg + d2H1dnbdg*rhoval)*fac
1037            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA)
1038     &             + d2H0dg2*rhoval*fac
1039     &             + d2H1dg2*rhoval*fac
1040            Cmat2(n,D2_GAA_GAB) = Cmat2(n,D2_GAA_GAB)
1041     &             + 2.0d0*d2H0dg2*rhoval*fac
1042     &             + 2.0d0*d2H1dg2*rhoval*fac
1043            Cmat2(n,D2_GAA_GBB) = Cmat2(n,D2_GAA_GBB)
1044     &             + d2H0dg2*rhoval*fac
1045     &             + d2H1dg2*rhoval*fac
1046            Cmat2(n,D2_GAB_GAB) = Cmat2(n,D2_GAB_GAB)
1047     &             + 4.0d0*d2H0dg2*rhoval*fac
1048     &             + 4.0d0*d2H1dg2*rhoval*fac
1049            Cmat2(n,D2_GAB_GBB) = Cmat2(n,D2_GAB_GBB)
1050     &             + 2.0d0*d2H0dg2*rhoval*fac
1051     &             + 2.0d0*d2H1dg2*rhoval*fac
1052            Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB)
1053     &             + d2H0dg2*rhoval*fac
1054     &             + d2H1dg2*rhoval*fac
1055#endif
1056         endif
1057   20 continue
1058c
1059      return
1060      end
1061c
1062#ifndef SECOND_DERIV
1063#define SECOND_DERIV
1064c
1065c     Compile source again for the 2nd derivative case
1066c
1067#include "xc_perdew91.F"
1068#endif
1069