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