1c Copyright 2019 (C) Orbital-free DFT group at University of Florida
2c Licensed under the Educational Community License, Version 2.0
3c (the "License"); you may not use this file except in compliance with
4c the License. You may obtain a copy of the License at
5c
6c    http://www.osedu.org/licenses/ECL-2.0
7c
8c Unless required by applicable law or agreed to in writing,
9c software distributed under the License is distributed on an "AS IS"
10c BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
11c or implied. See the License for the specific language governing
12c permissions and limitations under the License.
13c
14#include "dft2drv.fh"
15#include "dft3drv.fh"
16c     Nearly correct asymptotic potential (NCAP) functional
17c     (Exchange part only)
18c           GGA
19C         utilizes ingredients:
20c                              rho   -  density
21c                              delrho - gradient of density
22c
23c     Written by:
24c     Daniel Mejia-Rodriguez
25c     QTP, Department of Physics, University of Florida
26c
27c     References:
28c     J. Carmona-Espindola, J.L. Gazquez, A. Vela, S.B. Trickey
29c     JCTC 15, 303 (2019).
30c     DOI: 10.1021/acs.jctc.8b00998
31c
32c
33#if !defined SECOND_DERIV && !defined THIRD_DERIV
34      Subroutine xc_xncap(tol_rho, fac, rho, delrho,
35     &                    Amat, Cmat, nq, ipol, Ex, qwght,ldew,func)
36#elif defined(SECOND_DERIV) && !defined THIRD_DERIV
37      Subroutine xc_xncap_d2(tol_rho, fac, rho, delrho,
38     &                       Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
39     &                       qwght,ldew,func)
40#else
41      Subroutine xc_xncap_d3(tol_rho, fac, rho, delrho,
42     &                       Amat, Amat2, Amat3, Cmat, Cmat2, Cmat3,
43     &                       nq, ipol, Ex, qwght,ldew,func)
44#endif
45c
46      implicit none
47c
48      double precision fac, Ex
49      integer nq, ipol
50      logical ldew
51      double precision func(*)  ! value of the functional [output]
52c
53c     Charge Density & Its Cube Root
54c
55      double precision rho(nq,ipol*(ipol+1)/2)
56c
57c     Charge Density Gradient
58c
59      double precision delrho(nq,3,ipol)
60c
61c     Quadrature Weights
62c
63      double precision qwght(nq)
64c
65c     Sampling Matrices for the XC Potential & Energy
66c
67      double precision amat(nq,ipol), cmat(nq,*)
68c
69#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
70      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
71#endif
72#ifdef THIRD_DERIV
73      double precision Amat3(nq,NCOL_AMAT3), Cmat3(nq,NCOL_CMAT3)
74#endif
75c
76      double precision tol_rho, pi, mu, alpha, beta, zeta
77      double precision ckf, Ax
78      double precision F43, F13, F23, F49
79c
80      parameter(mu=0.2195149727645171d0, beta=0.01808569669d0)
81      parameter(zeta=0.30412141859531383d0)
82
83      parameter (F43=4.d0/3.d0, F13=1.d0/3.d0, F23=2.0d0/3.0d0)
84      parameter (F49=4d0/9d0)
85c
86      integer n
87      double precision rrho, rho43, rho13, gamma, gam12, s, s2, d1s
88      double precision arcsinh, darcsinh, d2arcsinh, d3arcsinh
89      double precision Fx, dFx, d2Fx, d3Fx
90      double precision tansin, dtansin
91      double precision factor, dfactor
92      double precision asymp, dasymp
93      double precision denom, ddenom
94
95#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
96      double precision d2denom, d2asymp, d2tansin, d2factor, d2s, rhom23
97#endif
98#if defined(THIRD_DERIV)
99      double precision d3denom, d3asymp, d3tansin, d3factor, d3s, rhom53
100#endif
101      double precision ops, omz, logops
102
103      arcsinh(s)=dlog(s+dsqrt(1d0+s*s))
104      darcsinh(s)=1d0/dsqrt(1d0+s*s)
105      d2arcsinh(s) = -s/dsqrt(1d0+s*s)**3
106      d3arcsinh(s) = (2d0*s*s - 1d0)/dsqrt(1d0+s*s)**5
107c
108      pi = acos(-1.d0)
109      ckf = (3d0*pi*pi)**F13
110      Ax = -3d0/(4d0*pi)*ckf
111      alpha = 4d0*pi/3d0*beta/mu
112c
113      if (ipol.eq.1 )then
114c
115c        ======> SPIN-RESTRICTED <======
116c
117#ifdef IFCV81
118CDEC$ NOSWP
119#endif
120         do 10 n = 1, nq
121            if (rho(n,1).lt.tol_rho) goto 10
122            rho43 = rho(n,1)**F43
123            rrho = 1d0/rho(n,1)
124            rho13 = rho43*rrho
125            gamma = delrho(n,1,1)*delrho(n,1,1) +
126     &              delrho(n,2,1)*delrho(n,2,1) +
127     &              delrho(n,3,1)*delrho(n,3,1)
128            if (dsqrt(gamma).le.tol_rho**2) goto 10
129
130            s = dsqrt(gamma)/(2d0*ckf*rho43)
131            s2 = s*s
132
133            if (s.lt.1d-8) then
134               Fx = 1d0 + mu*s2 + alpha*zeta*mu*s2*s +
135     &              mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s2
136               dFx = 2d0*mu*s + 3d0*alpha*zeta*mu*s2 +
137     &               4d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s
138               d2Fx = 2d0*mu + 6d0*alpha*zeta*mu*s +
139     &                12d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2
140               d3Fx = 6d0*alpha*zeta*mu +
141     &                24d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s
142            else
143              ops = 1d0 + s
144              omz = 1d0 - zeta
145              logops = dlog(ops)
146
147              asymp = 1d0 + alpha*(omz*s*logops + zeta*s)
148              dasymp = alpha*( (s + zeta)/ops + omz*logops )
149
150              if (s.gt.30d0) then
151                tansin = arcsinh(s)
152                dtansin = darcsinh(s)
153              else
154                tansin = dtanh(s)*arcsinh(s)
155                dtansin = arcsinh(s)/dcosh(s)**2 + dtanh(s)*darcsinh(s)
156              endif
157
158              denom = 1d0 + beta*tansin
159              ddenom = beta*dtansin
160
161              factor = tansin/denom
162              dfactor = (dtansin - factor*ddenom)/denom
163
164              Fx = 1d0 + mu*factor*asymp
165              dFx = mu*(dfactor*asymp + factor*dasymp)
166
167#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
168              if (s.gt.30d0) then
169                d2tansin = d2arcsinh(s)
170              else
171                d2tansin = 2d0*darcsinh(s)/dcosh(s)**2 +
172     &                     dtanh(s)*d2arcsinh(s) -
173     &                     2d0*tansin/dcosh(s)**2
174              endif
175              d2asymp = alpha*omz*(2d0 + s)/ops**2
176              d2denom = beta*d2tansin
177              d2factor = (d2tansin - 2d0*dfactor*ddenom -
178     &                    factor*d2denom)/denom
179              d2Fx = mu*(d2factor*asymp + 2d0*dfactor*dasymp +
180     &                   factor*d2asymp)
181
182#endif
183#if defined(THIRD_DERIV)
184              if (s.gt.30d0) then
185                d3tansin = d3arcsinh(s)
186              else
187                d3tansin = 3d0*d2arcsinh(s)/dcosh(s)**2 -
188     &                     4d0*darcsinh(s)*dtanh(s)/dcosh(s)**2 +
189     &                     dtanh(s)*d3arcsinh(s) -
190     &                     2d0*dtansin/dcosh(s)**2 +
191     &                     4d0*tansin*dtanh(s)/dcosh(s)**2
192              endif
193              d3asymp = -alpha*omz*(3d0 + s)/ops**3
194              d3denom = beta*d3tansin
195              d3factor = (d3tansin - 3d0*d2factor*ddenom -
196     &                    3d0*dfactor*d2denom - factor*d3denom)/denom
197              d3Fx = mu*(d3factor*asymp + 3d0*d2factor*dasymp +
198     &                   3d0*dfactor*d2asymp + factor*d3asymp)
199#endif
200            endif
201
202            Ex = Ex + Ax*rho43*Fx*qwght(n)*fac
203            if (ldew) func(n) = func(n) + Ax*rho43*Fx*fac
204
205            d1s = 0.5d0*s/gamma
206            Amat(n,1) = Amat(n,1) + F43*Ax*rho13*(Fx - s*dFx)*fac
207            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 2d0*Ax*rho43*dFx*d1s*fac
208
209#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
210
211            d2s = -0.5d0*d1s/gamma
212            rhom23 = rho13*rrho
213            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 2d0*F49*Ax*rhom23*
214     &                          (Fx - s*dFx + 4d0*s2*d2Fx)*fac
215
216            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) -
217     &                           (F43*Ax*rho13*d2Fx*d1s*s)*4d0*fac
218
219            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + Ax*rho43*
220     &                            (d2Fx*d1s**2 + dFx*d2s)*8d0*fac
221#endif
222#ifdef THIRD_DERIV
223            rhom53 = rhom23*rrho
224            d3s = -1.5d0*d2s/gamma
225
226            Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) - 4d0*fac*
227     &      F23*F49*Ax*rhom53*(Fx - s*dFx + 18d0*s2*d2Fx +8d0*s2*s*d3Fx)
228
229            Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) +
230     &      F49*Ax*rhom23*d1s*s*(7d0*d2Fx + 4d0*d3Fx*s)*8d0*fac
231
232            Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) -
233     &      F43*Ax*rho13*(d2Fx*d1s**2 + d3Fx*d1s**2*s +
234     &                     d2Fx*d2s*s)*16d0*fac
235
236            Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) +
237     &      Ax*rho43*(d3Fx*d1s**3 + 3d0*d2Fx*d1s*d2s +
238     &                dFx*d3s)*32d0*fac
239
240#endif
241 10      continue
242c
243      else
244c
245c        ======> SPIN-UNRESTRICTED <======
246c
247#ifdef IFCV81
248CDEC$ NOSWP
249#endif
250         do 20 n = 1, nq
251            if (rho(n,1).lt.tol_rho) goto 20
252c
253c     Alpha
254c
255            if (rho(n,2).lt.tol_rho) goto 25
256            rho43 = (2d0*rho(n,2))**F43
257            rrho = 0.5d0/rho(n,2)
258            rho13 = rho43*rrho
259c
260            gamma = delrho(n,1,1)*delrho(n,1,1) +
261     &              delrho(n,2,1)*delrho(n,2,1) +
262     &              delrho(n,3,1)*delrho(n,3,1)
263            gam12 = 2d0*dsqrt(gamma)
264            if (gam12.le.tol_rho**2) goto 25
265c
266            s = gam12/(2d0*ckf*rho43)
267            s2 = s*s
268
269            d1s = 0.5d0*s/gamma
270
271            if (s.lt.1d-8) then
272               Fx = 1d0 + mu*s2 + alpha*zeta*mu*s2*s +
273     &              mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s2
274               dFx = 2d0*mu*s + 3d0*alpha*zeta*mu*s2 +
275     &               4d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s
276               d2Fx = 2d0*mu + 6d0*alpha*zeta*mu*s +
277     &                12d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2
278               d3Fx = 6d0*alpha*zeta*mu +
279     &                24d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s
280            else
281              ops = 1d0 + s
282              omz = 1d0 - zeta
283              logops = dlog(ops)
284
285              asymp = 1d0 + alpha*(omz*s*logops + zeta*s)
286              dasymp = alpha*( (s + zeta)/ops + omz*logops )
287
288              if (s.gt.30d0) then
289                tansin = arcsinh(s)
290                dtansin = darcsinh(s)
291              else
292                tansin = dtanh(s)*arcsinh(s)
293                dtansin = arcsinh(s)/dcosh(s)**2 + dtanh(s)*darcsinh(s)
294              endif
295
296              denom = 1d0 + beta*tansin
297              ddenom = beta*dtansin
298
299              factor = tansin/denom
300              dfactor = (dtansin - factor*ddenom)/denom
301
302              Fx = 1d0 + mu*factor*asymp
303              dFx = mu*(dfactor*asymp + factor*dasymp)
304#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
305              if (s.gt.30d0) then
306                d2tansin = d2arcsinh(s)
307              else
308                d2tansin = 2d0*darcsinh(s)/dcosh(s)**2 +
309     &                     dtanh(s)*d2arcsinh(s) -
310     &                     2d0*tansin/dcosh(s)**2
311              endif
312              d2asymp = alpha*omz*(2d0 + s)/ops**2
313              d2denom = beta*d2tansin
314              d2factor = (d2tansin - 2d0*dfactor*ddenom -
315     &                    factor*d2denom)/denom
316              d2Fx = mu*(d2factor*asymp + 2d0*dfactor*dasymp +
317     &                   factor*d2asymp)
318#endif
319#if defined(THIRD_DERIV)
320              if (s.gt.30d0) then
321                d3tansin = d3arcsinh(s)
322              else
323                d3tansin = 3d0*d2arcsinh(s)/dcosh(s)**2 -
324     &                     4d0*darcsinh(s)*dtanh(s)/dcosh(s)**2 +
325     &                     dtanh(s)*d3arcsinh(s) -
326     &                     2d0*dtansin/dcosh(s)**2 +
327     &                     4d0*tansin*dtanh(s)/dcosh(s)**2
328              endif
329              d3asymp = -alpha*omz*(3d0 + s)/ops**3
330              d3denom = beta*d3tansin
331              d3factor = (d3tansin - 3d0*d2factor*ddenom -
332     &                    3d0*dfactor*d2denom - factor*d3denom)/denom
333              d3Fx = mu*(d3factor*asymp + 3d0*d2factor*dasymp +
334     &                   3d0*dfactor*d2asymp + factor*d3asymp)
335#endif
336            endif
337
338            Ex = Ex + 0.5d0*Ax*rho43*Fx*qwght(n)*fac
339            if (ldew) func(n) = func(n) + 0.5d0*Ax*rho43*Fx*fac
340
341            Amat(n,1) = Amat(n,1) + F43*Ax*rho13*(Fx - s*dFx)*fac
342            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) +
343     &                       0.5d0*Ax*rho43*dFx*d1s*fac
344
345#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
346            d2s = -0.5d0*d1s/gamma
347            rhom23 = rho13*rrho
348            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 2d0*F49*Ax*rhom23*
349     &                          (Fx - s*dFx + 4d0*s2*d2Fx)*fac
350
351            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) -
352     &                           (F43*Ax*rho13*d2Fx*d1s*s)*fac
353
354            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + fac*Ax*rho43*
355     &                            (d2Fx*d1s**2 + dFx*d2s)*0.5d0
356#endif
357#ifdef THIRD_DERIV
358            rhom53 = rhom23*rrho
359            d3s = -1.5d0*d2s/gamma
360
361            Amat3(n,D3_RA_RA_RA) = Amat3(n,D3_RA_RA_RA) - 4d0*fac*
362     &      F23*F49*Ax*rhom53*(Fx - s*dFx + 18d0*s2*d2Fx +8d0*s2*s*d3Fx)
363
364            Cmat3(n,D3_RA_RA_GAA) = Cmat3(n,D3_RA_RA_GAA) +
365     &      F49*Ax*rhom23*d1s*s*(7d0*d2Fx + 4d0*d3Fx*s)*2d0*fac
366
367            Cmat3(n,D3_RA_GAA_GAA) = Cmat3(n,D3_RA_GAA_GAA) -
368     &      F43*Ax*rho13*(d2Fx*d1s**2 + d3Fx*d1s**2*s +
369     &                     d2Fx*d2s*s)*fac
370
371            Cmat3(n,D3_GAA_GAA_GAA) = Cmat3(n,D3_GAA_GAA_GAA) +
372     &      Ax*rho43*(d3Fx*d1s**3 + 3d0*d2Fx*d1s*d2s +
373     &                dFx*d3s)*0.5d0*fac
374
375#endif
376c
377c     Beta
378c
379 25         continue
380            if (rho(n,3).lt.tol_rho) goto 20
381            rho43 = (2d0*rho(n,3))**F43
382            rrho = 0.5d0/rho(n,3)
383            rho13 = rho43*rrho
384c
385            gamma = delrho(n,1,2)*delrho(n,1,2) +
386     &              delrho(n,2,2)*delrho(n,2,2) +
387     &              delrho(n,3,2)*delrho(n,3,2)
388            gam12 = 2d0*dsqrt(gamma)
389            if (gam12.le.tol_rho**2) goto 20
390c
391            s = gam12/(2d0*ckf*rho43)
392            s2 = s*s
393
394            d1s = 0.5d0*s/gamma
395
396            if (s.lt.1d-8) then
397               Fx = 1d0 + mu*s2 + alpha*zeta*mu*s2*s +
398     &              mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s2
399               dFx = 2d0*mu*s + 3d0*alpha*zeta*mu*s2 +
400     &               4d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2*s
401               d2Fx = 2d0*mu + 6d0*alpha*zeta*mu*s +
402     &                12d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s2
403               d3Fx = 6d0*alpha*zeta*mu +
404     &                24d0*mu*(alpha - beta - alpha*zeta - 0.5d0)*s
405            else
406              ops = 1d0 + s
407              omz = 1d0 - zeta
408              logops = dlog(ops)
409
410              asymp = 1d0 + alpha*(omz*s*logops + zeta*s)
411              dasymp = alpha*( (s + zeta)/ops + omz*logops )
412
413              if (s.gt.30d0) then
414                tansin = arcsinh(s)
415                dtansin = darcsinh(s)
416              else
417                tansin = dtanh(s)*arcsinh(s)
418                dtansin = arcsinh(s)/dcosh(s)**2 + dtanh(s)*darcsinh(s)
419              endif
420
421              denom = 1d0 + beta*tansin
422              ddenom = beta*dtansin
423
424              factor = tansin/denom
425              dfactor = (dtansin - factor*ddenom)/denom
426
427              Fx = 1d0 + mu*factor*asymp
428              dFx = mu*(dfactor*asymp + factor*dasymp)
429
430#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
431              if (s.gt.30d0) then
432                d2tansin = d2arcsinh(s)
433              else
434                d2tansin = 2d0*darcsinh(s)/dcosh(s)**2 +
435     &                     dtanh(s)*d2arcsinh(s) -
436     &                     2d0*tansin/dcosh(s)**2
437              endif
438              d2asymp = alpha*omz*(2d0 + s)/ops**2
439              d2denom = beta*d2tansin
440              d2factor = (d2tansin - 2d0*dfactor*ddenom -
441     &                    factor*d2denom)/denom
442              d2Fx = mu*(d2factor*asymp + 2d0*dfactor*dasymp +
443     &                   factor*d2asymp)
444#endif
445#if defined(THIRD_DERIV)
446
447              if (s.gt.30d0) then
448                d3tansin = d3arcsinh(s)
449              else
450                d3tansin = 3d0*d2arcsinh(s)/dcosh(s)**2 -
451     &                     4d0*darcsinh(s)*dtanh(s)/dcosh(s)**2 +
452     &                     dtanh(s)*d3arcsinh(s) -
453     &                     2d0*dtansin/dcosh(s)**2 +
454     &                     4d0*tansin*dtanh(s)/dcosh(s)**2
455              endif
456              d3asymp = -alpha*omz*(3d0 + s)/ops**3
457              d3denom = beta*d3tansin
458              d3factor = (d3tansin - 3d0*d2factor*ddenom -
459     &                    3d0*dfactor*d2denom - factor*d3denom)/denom
460              d3Fx = mu*(d3factor*asymp + 3d0*d2factor*dasymp +
461     &                   3d0*dfactor*d2asymp + factor*d3asymp)
462#endif
463            endif
464
465            Ex = Ex + 0.5d0*Ax*rho43*Fx*qwght(n)*fac
466            if (ldew) func(n) = func(n) + 0.5d0*Ax*rho43*Fx*fac
467
468            Amat(n,2) = Amat(n,2) + F43*Ax*rho13*(Fx - s*dFx)*fac
469            Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + 0.5d0*Ax*rho43*dFx*d1s*fac
470
471#if defined(SECOND_DERIV) || defined(THIRD_DERIV)
472            d2s = -0.5d0*d1s/gamma
473            rhom23 = rho13*rrho
474            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + 2d0*F49*Ax*rhom23*
475     &                          (Fx - s*dFx + 4d0*s2*d2Fx)*fac
476
477            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) -
478     &                           (F43*Ax*rho13*d2Fx*d1s*s)*fac
479
480            Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + fac*Ax*rho43*
481     &                            (d2Fx*d1s**2 + dFx*d2s)*0.5d0
482#endif
483#ifdef THIRD_DERIV
484            rhom53 = rhom23*rrho
485            d3s = -1.5d0*d2s/gamma
486
487            Amat3(n,D3_RB_RB_RB) = Amat3(n,D3_RB_RB_RB) - 4d0*fac*
488     &      F23*F49*Ax*rhom53*(Fx - s*dFx + 18d0*s2*d2Fx +8d0*s2*s*d3Fx)
489
490            Cmat3(n,D3_RB_RB_GBB) = Cmat3(n,D3_RB_RB_GBB) +
491     &      F49*Ax*rhom23*d1s*s*(7d0*d2Fx + 4d0*d3Fx*s)*2d0*fac
492
493            Cmat3(n,D3_RB_GBB_GBB) = Cmat3(n,D3_RB_GBB_GBB) -
494     &      F43*Ax*rho13*(d2Fx*d1s**2 + d3Fx*d1s**2*s +
495     &                     d2Fx*d2s*s)*fac
496
497            Cmat3(n,D3_GBB_GBB_GBB) = Cmat3(n,D3_GBB_GBB_GBB) +
498     &      Ax*rho43*(d3Fx*d1s**3 + 3d0*d2Fx*d1s*d2s +
499     &                dFx*d3s)*0.5d0*fac
500#endif
501c
502 20      continue
503      endif
504c
505      return
506      end
507#ifndef SECOND_DERIV
508#define SECOND_DERIV
509c
510c     Compile source again for the 2nd derivative case
511c
512#include "xc_xncap.F"
513#endif
514#ifndef THIRD_DERIV
515#define THIRD_DERIV
516c
517c     Compile source again for the 3rd derivative case
518c
519#include "xc_xncap.F"
520#endif
521c $Id$
522