1c
2c     this is the gradient correction term by Keal and Tozer,
3c     which can be used on its own (as in the KT1 functional),
4c     or forms part of the SSB-D functional
5c
6c     it never includes LDA !!
7c     (this is taken care of somewhere else)
8c
9c     note that even though the energy is assigned here to the
10c     exchange, this gradient correction is strictly speaking
11c     NOT an exchange functional !
12c
13c     KT1/KT2 reference:
14c     T.W. Keal, D.J. Tozer, JCP 2006, 119, 3015
15c     SSB-D reference:
16c     M. Swart, M. Sola, F.M. Bickelhaupt, JCP 2009, 131, 094103
17c
18#ifndef SECOND_DERIV
19      Subroutine xc_kt1(tol_rho, fac, rho, delrho,
20     &                     Amat, Cmat, nq, ipol, Ex, qwght,ldew,func)
21#else
22      Subroutine xc_kt1_d2(tol_rho, fac, rho, delrho,
23     &                     Amat, Amat2, Cmat, Cmat2, nq, ipol,
24     &                     Ex, qwght,ldew,func)
25#endif
26c
27C$Id$
28c
29      implicit none
30c
31#include "dft2drv.fh"
32c
33      double precision tol_rho, fac, Ex
34      integer nq, ipol
35      logical ldew
36      double precision func(*)  ! value of the functional [output]
37c
38c     Charge Density
39c
40      double precision rho(nq,ipol*(ipol+1)/2)
41c
42c     Charge Density Gradient
43c
44      double precision delrho(nq,3,ipol)
45c
46c     Quadrature Weights
47c
48      double precision qwght(nq)
49c
50c     Sampling Matrices for the XC Potential
51c
52      double precision Amat(nq,ipol), Cmat(nq,*)
53c
54#ifdef SECOND_DERIV
55c
56c     Second Derivatives of the Exchange Energy Functional
57c
58      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
59#endif
60c
61      double precision DELTA, GAMKT
62      Parameter (DELTA = 0.1D0, GAMKT= -0.006d0)
63c
64c References:
65c
66c    Keal, Tozer, JCP 119, 3015 (2003), JCP 121, 5654 (2004)
67c    Swart, Sola, Bickelhaupt, JCP 131, XXXX (2009)
68c    Johnson, Gill & Pople, J. Chem. Phys. 98, 5612 (1993)
69c
70c***************************************************************************
71c
72      integer n
73      double precision hrho
74      double precision rho13, rho43, gamma, g, gdenom, gdenom2
75#ifdef SECOND_DERIV
76      double precision rho23, rhom23, gdenom3
77#endif
78c
79c     NOTE: the gamma from the KT1 formulation is here called
80c           gamkt, gamma is in NWChem reserved for grad**2
81c
82      if (ipol.eq.1) then
83c
84c        ======> SPIN-RESTRICTED <======
85c
86         do 10 n = 1, nq
87            if (rho(n,1).lt.tol_rho) goto 10
88c
89c           Spin alpha:
90c
91            hrho  = 0.5d0*rho(n,1)
92            rho13 = hrho**(1.d0/3.d0)
93            rho43 = rho13*hrho
94            gamma = delrho(n,1,1)*delrho(n,1,1) +
95     &              delrho(n,2,1)*delrho(n,2,1) +
96     &              delrho(n,3,1)*delrho(n,3,1)
97            if (dsqrt(gamma).gt.tol_rho) then
98               gamma = 0.25d0 * gamma
99            else
100               goto 10
101            endif
102c
103            gdenom = 1.d0 / (rho43 + DELTA)
104            gdenom2 = gdenom*gdenom
105            g = GAMKT * gamma * gdenom
106c
107            Ex = Ex + 2.d0*g*qwght(n)*fac
108            if (ldew)  func(n) = func(n) + 2.d0*g*fac
109            Amat(n,1) = Amat(n,1) - (4d0/3d0)*GAMKT*gamma*rho13*
110     &                  fac*gdenom2
111            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + GAMKT*gdenom*fac
112c
113#ifdef SECOND_DERIV
114            rho23 = rho13*rho13
115            rhom23 = rho13 / (0.5d0*rho(n,1))
116            gdenom3 = gdenom2*gdenom
117c
118            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + (4d0/3d0)*GAMKT*
119     &          gamma*(rho23*7d0/3d0 - DELTA*rhom23/3d0)*gdenom3*fac
120            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA)
121     &           - (4d0/3d0)*GAMKT*rho13*gdenom2*fac
122c
123c      second derivative w.r.t. gamma is zero !
124c      (by construction)
125c      therefore, nothing added to Cmat2(n,D2_GAA_GAA)
126c
127#endif
128c
129 10      continue
130c
131      else
132c
133c        ======> SPIN-UNRESTRICTED <======
134c
135         do 20 n = 1, nq
136            if (rho(n,1).lt.tol_rho) goto 20
137            if (rho(n,2).lt.tol_rho) goto 25
138c
139c           Spin alpha:
140c
141            rho13 = rho(n,2)**(1.d0/3.d0)
142            rho43 = rho13*rho(n,2)
143            gamma = delrho(n,1,1)*delrho(n,1,1) +
144     &              delrho(n,2,1)*delrho(n,2,1) +
145     &              delrho(n,3,1)*delrho(n,3,1)
146            if (dsqrt(gamma).lt.tol_rho) then
147               goto 25
148            endif
149c
150            gdenom = 1d0 / (rho43 + DELTA)
151            gdenom2 = gdenom*gdenom
152            g = GAMKT * gamma * gdenom
153c
154            Ex = Ex + g*qwght(n)*fac
155            if (ldew)  func(n) = func(n) + g*fac
156            Amat(n,1) = Amat(n,1) - (4d0/3d0)*GAMKT*gamma*rho13*
157     &                 gdenom2*fac
158            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + GAMKT*gdenom*fac
159c
160#ifdef SECOND_DERIV
161            rho23 = rho13*rho13
162            rhom23 = rho13 / rho(n,2)
163            gdenom3 = gdenom2*gdenom
164c
165            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + (4d0/3d0)*GAMKT*
166     &          gamma*(rho23*7d0/3d0 - DELTA*rhom23/3d0)*gdenom3*fac
167            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA)
168     &           - (4d0/3d0)*GAMKT*rho13*gdenom2*fac
169c
170c      second derivative w.r.t. gamma is zero !
171c      (by construction)
172c      therefore, nothing added to Cmat2(n,D2_GAA_GAA)
173#endif
174c
175 25         continue
176c
177c           Spin beta:
178c
179            if (rho(n,3).lt.tol_rho) goto 20
180c
181            rho13 = rho(n,3)**(1.d0/3.d0)
182            rho43 = rho13*rho(n,3)
183            gamma = delrho(n,1,2)*delrho(n,1,2) +
184     &              delrho(n,2,2)*delrho(n,2,2) +
185     &              delrho(n,3,2)*delrho(n,3,2)
186            if (dsqrt(gamma).lt.tol_rho) then
187               goto 20
188            endif
189c
190            gdenom = 1d0 / (rho43 + DELTA)
191            gdenom2 = gdenom*gdenom
192            g = GAMKT * gamma * gdenom
193c
194            Ex = Ex + g*qwght(n)*fac
195            if (ldew)  func(n) = func(n) + g*fac
196            Amat(n,2) = Amat(n,2) - (4d0/3d0)*GAMKT*gamma*rho13*
197     &            fac*gdenom2
198            Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + GAMKT*gdenom*fac
199c
200#ifdef SECOND_DERIV
201            rho23 = rho13*rho13
202            rhom23 = rho13 / rho(n,3)
203            gdenom3 = gdenom2*gdenom
204c
205            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + (4d0/3d0)*GAMKT*
206     &          gamma*(rho23*7d0/3d0 - DELTA*rhom23/3d0)*gdenom3*fac
207            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB)
208     &           - (4d0/3d0)*GAMKT*rho13*gdenom2*fac
209c
210c      second derivative w.r.t. gamma is zero !
211c      (by construction)
212c      therefore, nothing added to Cmat2(n,D2_GAA_GAA)
213c
214#endif
215c
216 20      continue
217c
218      endif
219c
220      return
221      end
222#ifndef SECOND_DERIV
223#define SECOND_DERIV
224c
225c     Compile source again for the 2nd derivative case
226c
227#include "xc_kt1.F"
228#endif
229