1#ifndef SECOND_DERIV 2 Subroutine xc_gill96(tol_rho, fac, lfac, nlfac, rho, delrho, 3 & Amat, Cmat, nq, ipol, Ex, qwght,ldew,func) 4#else 5#include "dft2drv.fh" 6 Subroutine xc_gill96_d2(tol_rho, fac, lfac, nlfac, rho, delrho, 7 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 8 & qwght,ldew,func) 9#endif 10c 11C$Id$ 12c 13 implicit none 14c 15c 16 double precision tol_rho, fac, Ex 17 integer nq, ipol 18 logical lfac, nlfac,ldew 19 double precision func(*) ! value of the functional [output] 20c 21c Charge Density 22c 23 double precision rho(nq,ipol*(ipol+1)/2) 24c 25c Charge Density Gradient 26c 27 double precision delrho(nq,3,ipol) 28c 29c Quadrature Weights 30c 31 double precision qwght(nq) 32c 33c Sampling Matrices for the XC Potential 34c 35 double precision Amat(nq,ipol), Cmat(nq,*) 36#ifdef SECOND_DERIV 37 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 38#endif 39c 40 double precision BETA, ONE3, FOUR3 41#ifdef SECOND_DERIV 42 double precision SEVEN3 43#endif 44 Parameter (BETA = 1d0/137d0) 45 Parameter (ONE3 = 1d0/3d0, FOUR3 = 4d0/3d0) 46#ifdef SECOND_DERIV 47 Parameter (SEVEN3 = 7d0/3d0) 48#endif 49c 50c References: 51c 52c Gill , Mol. Phys. 89, 433 (1996) 53c 54c*************************************************************************** 55c 56 integer n 57 double precision C, rhoval, rho13, rho43, gamma, x, d1x(2), 58 & g, d1g 59#ifdef SECOND_DERIV 60 double precision d2x(3), d2g 61#endif 62c 63c 64c Uniform electron gas constant 65c 66 C = -(1.5d0)*(0.75d0/acos(-1d0))**(ONE3) 67c 68 if (ipol.eq.1) then 69c 70c ======> SPIN-RESTRICTED <====== 71c 72 do 10 n = 1, nq 73 if (rho(n,1).lt.tol_rho) goto 10 74 rhoval = 0.5d0*rho(n,1) 75c 76c Spin alpha: 77c 78 rho13 = rhoval**ONE3 79 rho43 = rho13*rhoval 80c Include factor of 4/3 in rho13 since it always appears with it 81 rho13 = FOUR3*rho13 82c 83 if (lfac) then 84 Ex = Ex + 2d0*rho43*C*qwght(n)*fac 85 if(ldew)func(n) = func(n) + 2.d0*rho43*C*fac 86 Amat(n,1) = Amat(n,1) + rho13*C*fac 87#ifdef SECOND_DERIV 88 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 89 & + ONE3*(rho13/rhoval)*C*fac 90#endif 91 endif 92c 93 gamma = delrho(n,1,1)*delrho(n,1,1) + 94 & delrho(n,2,1)*delrho(n,2,1) + 95 & delrho(n,3,1)*delrho(n,3,1) 96 if (nlfac.and.dsqrt(gamma).gt.tol_rho)then 97 gamma = 0.25d0*gamma 98 x = dsqrt(gamma)/rho43 99 d1x(1) = -FOUR3*x/rhoval 100 d1x(2) = 0.5d0*x/gamma 101 g = -BETA*x*sqrt(x) 102 d1g = -1.5d0*BETA*sqrt(x) 103c 104 Ex = Ex + 2d0*rho43*g*qwght(n)*fac 105 if(ldew)func(n) = func(n) + 2.d0*rho43*g*fac 106 Amat(n,1) = Amat(n,1) + (rho13*g+rho43*d1g*d1x(1))*fac 107 Cmat(n,1) = Cmat(n,1) + rho43*d1g*d1x(2)*fac 108#ifdef SECOND_DERIV 109 d2g = 0.5d0*d1g/x 110 d2x(1) = -SEVEN3*d1x(1)/rhoval 111 d2x(2) = -FOUR3*d1x(2)/rhoval 112 d2x(3) = -0.5d0*d1x(2)/gamma 113 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 114 & + ONE3*(rho13/rhoval)*g*fac 115 & + 2.d0*rho13*d1g*d1x(1)*fac 116 & + rho43*(d2g*d1x(1)*d1x(1)+d1g*d2x(1))*fac 117 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 118 & + rho13*d1g*d1x(2)*fac 119 & + rho43*(d2g*d1x(1)*d1x(2)+d1g*d2x(2))*fac 120 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 121 & + rho43*(d2g*d1x(2)*d1x(2)+d1g*d2x(3))*fac 122#endif 123 endif 124c 125 10 continue 126c 127 else 128c 129c ======> SPIN-UNRESTRICTED <====== 130c 131 do 20 n = 1, nq 132 if (rho(n,1).lt.tol_rho) goto 20 133 if (rho(n,2).lt.tol_rho) goto 25 134c 135c Spin alpha: 136c 137 rhoval = rho(n,2) 138 rho13 = rhoval**ONE3 139 rho43 = rho13*rhoval 140c Include factor of 4/3 in rho13 since it always appears with it 141 rho13 = FOUR3*rho13 142c 143 if (lfac) then 144 Ex = Ex + rho43*C*qwght(n)*fac 145 if(ldew)func(n) = func(n) + rho43*C*fac 146 Amat(n,1) = Amat(n,1) + rho13*C*fac 147#ifdef SECOND_DERIV 148 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 149 & + ONE3*(rho13/rhoval)*C*fac 150#endif 151 endif 152c 153 gamma = delrho(n,1,1)*delrho(n,1,1) + 154 & delrho(n,2,1)*delrho(n,2,1) + 155 & delrho(n,3,1)*delrho(n,3,1) 156 if (nlfac.and.dsqrt(gamma).gt.tol_rho)then 157 x = dsqrt(gamma)/rho43 158 d1x(1) = -FOUR3*x/rhoval 159 d1x(2) = 0.5d0*x/gamma 160 g = -BETA*x*sqrt(x) 161 d1g = -1.5d0*BETA*sqrt(x) 162c 163 Ex = Ex + rho43*g*qwght(n)*fac 164 if(ldew)func(n) = func(n) + rho43*g*fac 165 Amat(n,1) = Amat(n,1) + (rho13*g+rho43*d1g*d1x(1))*fac 166 Cmat(n,1) = Cmat(n,1) + rho43*d1g*d1x(2)*fac 167#ifdef SECOND_DERIV 168 d2g = 0.5d0*d1g/x 169 d2x(1) = -SEVEN3*d1x(1)/rhoval 170 d2x(2) = -FOUR3*d1x(2)/rhoval 171 d2x(3) = -0.5d0*d1x(2)/gamma 172 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 173 & + ONE3*(rho13/rhoval)*g*fac 174 & + 2.d0*rho13*d1g*d1x(1)*fac 175 & + rho43*(d2g*d1x(1)*d1x(1)+d1g*d2x(1))*fac 176 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) 177 & + rho13*d1g*d1x(2)*fac 178 & + rho43*(d2g*d1x(1)*d1x(2)+d1g*d2x(2))*fac 179 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) 180 & + rho43*(d2g*d1x(2)*d1x(2)+d1g*d2x(3))*fac 181#endif 182 endif 183c 184c 185 25 continue 186c 187c Spin beta: 188c 189 if (rho(n,3).lt.tol_rho) goto 20 190c 191 rhoval = rho(n,3) 192 rho13 = rhoval**ONE3 193 rho43 = rho13*rhoval 194c Include factor of 4/3 in rho13 since it always appears with it 195 rho13 = FOUR3*rho13 196c 197 if (lfac) then 198 Ex = Ex + rho43*C*qwght(n)*fac 199 if(ldew)func(n) = func(n) + rho43*C*fac 200 Amat(n,2) = Amat(n,2) + rho13*C*fac 201#ifdef SECOND_DERIV 202 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 203 & + ONE3*(rho13/rhoval)*C*fac 204#endif 205 endif 206c 207 gamma = delrho(n,1,2)*delrho(n,1,2) + 208 & delrho(n,2,2)*delrho(n,2,2) + 209 & delrho(n,3,2)*delrho(n,3,2) 210 if (nlfac.and.dsqrt(gamma).gt.tol_rho)then 211 x = dsqrt(gamma)/rho43 212 d1x(1) = -FOUR3*x/rhoval 213 d1x(2) = 0.5d0*x/gamma 214 g = -BETA*x*sqrt(x) 215 d1g = -1.5d0*BETA*sqrt(x) 216c 217 Ex = Ex + rho43*g*qwght(n)*fac 218 if(ldew)func(n) = func(n) + rho43*g*fac 219 Amat(n,2) = Amat(n,2) + (rho13*g+rho43*d1g*d1x(1))*fac 220 Cmat(n,3) = Cmat(n,3) + rho43*d1g*d1x(2)*fac 221#ifdef SECOND_DERIV 222 d2g = 0.5d0*d1g/x 223 d2x(1) = -SEVEN3*d1x(1)/rhoval 224 d2x(2) = -FOUR3*d1x(2)/rhoval 225 d2x(3) = -0.5d0*d1x(2)/gamma 226 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 227 & + ONE3*(rho13/rhoval)*g*fac 228 & + 2.d0*rho13*d1g*d1x(1)*fac 229 & + rho43*(d2g*d1x(1)*d1x(1)+d1g*d2x(1))*fac 230 Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) 231 & + rho13*d1g*d1x(2)*fac 232 & + rho43*(d2g*d1x(1)*d1x(2)+d1g*d2x(2))*fac 233 Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) 234 & + rho43*(d2g*d1x(2)*d1x(2)+d1g*d2x(3))*fac 235#endif 236 endif 237c 238c 239 20 continue 240c 241 endif 242c 243 return 244 end 245#ifndef SECOND_DERIV 246#define SECOND_DERIV 247c 248c Compile source again for the 2nd derivative case 249c 250#include "xc_gill96.F" 251#endif 252