1c Perdew-Wang 86 exchange functional 2c 3c References: 4c [a] J. P. Perdew and W. Yue (actually, Y. Wang) 5c Phys. Rev. B 33 (1986) 8800 6c 7#ifndef SECOND_DERIV 8 Subroutine xc_xpw86(tol_rho, fac, lfac, nlfac, rho, delrho, 9 & Amat, Cmat, nq, ipol, Ex, qwght,ldew,func) 10#else 11 Subroutine xc_xpw86_d2(tol_rho, fac, lfac, nlfac, rho, delrho, 12 & Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex, 13 & qwght,ldew,func) 14#endif 15c 16c$Id$ 17c 18 implicit none 19c 20#include "dft2drv.fh" 21c 22 double precision tol_rho ! zero-density crit. (input) 23 double precision fac ! energy multiplier (input) 24 logical lfac, nlfac ! local and non-local labels (input) 25 logical ldew ! compute funct(n) (input) 26 double precision Ex ! exchange energy (output) 27 integer nq ! number of grid points (input) 28 integer ipol ! spinpolarization flag (input) 29 double precision func(*) ! value of the functional [output] 30 double precision rho(nq,ipol*(ipol+1)/2) ! charge density (in) 31 double precision delrho(nq,3,ipol) ! gradient of rho (in) 32 double precision qwght(nq) ! quadrature weights (in) 33 double precision amat(nq,ipol) ! dex/drho matrix 34 double precision cmat(nq,*) ! 2*dex/(dgrho) 35#ifdef SECOND_DERIV 36c second derivatives? dkdc for now. 37 double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2) 38#endif 39 40c constants for PW86 41 double precision a, b, c, m1, Ax, F43, s_pfac 42 parameter(a = 1.296d0, b = 14d0, c = 0.2d0) 43 parameter(m1 = 1d0/15d0) 44 parameter(s_pfac = 6.18733545256027d0) 45 parameter(Ax = -0.738558766382022d0) 46 parameter(F43 = 4d0/3d0) 47 48 double precision pi, axr43, rh 49 double precision s, s2, s3, s4, s5, s6 50 double precision fs, sx, drho, drho2, dfds,oneminus 51 integer n 52c 53 54 55 pi = acos(-1.d0) 56 57c spin-restricted 58 if (ipol.eq.1) then 59#ifdef IFCV81 60CDEC$ NOSWP 61#endif 62 do 10 n = 1, nq 63 if (rho(n,1).lt.tol_rho) goto 10 64 65 rh = rho(n,1) 66 drho2 = (delrho(n,1,1)*delrho(n,1,1) + 67 & delrho(n,2,1)*delrho(n,2,1) + 68 & delrho(n,3,1)*delrho(n,3,1)) 69 drho = dsqrt(drho2) 70 if (drho.gt.tol_rho) then 71 s = drho / (s_pfac*rh**F43) 72 s2 = s * s 73 s3 = s2 * s 74 s4 = s3 * s 75 s5 = s4 * s 76 s6 = s5 * s 77 fs = (1 + a*s2 + b*s4 + c*s6)**m1 78 dfds = (2*a*s + 4*b*s3 + 6*c*s5)/(15.d0*fs**14) 79 oneminus=1d0-dfds/fs*s 80 else 81 s=0d0 82 fs=1d0 83 dfds=0d0 84 oneminus=1d0 85 endif 86 87 axr43 = Ax * rh**F43 88 sx = axr43 * fs 89 Ex = Ex + sx * qwght(n) * fac 90 if (ldew) func(n) = func(n) + sx * fac 91 92c GC: dex/drho 93 Amat(n,1) = Amat(n,1) + F43*sx/rh * oneminus 94 if(drho.gt.tol_rho) then 95c GC: 2*dex/d(grho2) 96 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + sx/fs*dfds*s/drho2 97 endif 98 99c xxxx missing 100#ifdef SECOND_DERIV 101 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 0 102 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + 0 103 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + 0 104#endif 105 10 continue 106c 107 else 108c spin unrestricted 109#ifdef IFCV81 110CDEC$ NOSWP 111#endif 112 do 100 n = 1, nq 113 if (rho(n,1).lt.tol_rho) goto 100 114 115c alpha 116 if (rho(n,2).lt.tol_rho) goto 150 117 rh = rho(n,2) * 2 118 drho2 = (delrho(n,1,1)*delrho(n,1,1) + 119 & delrho(n,2,1)*delrho(n,2,1) + 120 & delrho(n,3,1)*delrho(n,3,1)) * 4 121 drho = dsqrt(drho2) 122 s = drho / (s_pfac*rh**F43) 123 s2 = s * s 124 s3 = s2 * s 125 s4 = s3 * s 126 s5 = s4 * s 127 s6 = s5 * s 128 129 axr43 = Ax * rh**F43 130 fs = (1 + a*s2 + b*s4 + c*s6)**m1 131 sx = axr43 * fs 132 Ex = Ex + 0.5d0 * sx * qwght(n) * fac 133 if (ldew) func(n) = func(n) + 0.5d0 * sx * fac 134 135 dfds = (2*a*s + 4*b*s3 + 6*c*s5)/(15.d0*fs**14) 136c GC: dex/drho 137 Amat(n,1) = Amat(n,1) + F43*sx/rh * (1d0-dfds/fs*s) 138c GC: 2*dex/d(grho2) 139 if(drho2.gt.tol_rho) then 140 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + sx/fs*dfds*s/drho2 141 endif 142 143c xxxx missing 144#ifdef SECOND_DERIV 145 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 0 146 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + 0 147 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + 0 148#endif 149 150 150 continue 151c beta 152 if (rho(n,3).lt.tol_rho) goto 100 153 rh = rho(n,3) * 2 154 drho2 = (delrho(n,1,2)*delrho(n,1,2) + 155 & delrho(n,2,2)*delrho(n,2,2) + 156 & delrho(n,3,2)*delrho(n,3,2)) * 4 157 drho = dsqrt(drho2) 158 s = drho / (s_pfac*rh**F43) 159 s2 = s * s 160 s3 = s2 * s 161 s4 = s3 * s 162 s5 = s4 * s 163 s6 = s5 * s 164 165 axr43 = Ax * rh**F43 166 fs = (1 + a*s2 + b*s4 + c*s6)**m1 167 sx = axr43 * fs 168 Ex = Ex + 0.5d0 * sx * qwght(n) * fac 169 if (ldew) func(n) = func(n) + 0.5d0 * sx * fac 170 171 dfds = (2*a*s + 4*b*s3 + 6*c*s5)/(15.d0*fs**14) 172c GC: dex/drho 173 Amat(n,2) = Amat(n,2) + F43*sx/rh * (1d0-dfds/fs*s) 174c GC: 2*dex/d(grho2) 175 if(drho2.gt.tol_rho) then 176 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + sx/fs*dfds*s/drho2 177 endif 178 179c xxxx missing 180#ifdef SECOND_DERIV 181 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 0 182 Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + 0 183 Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + 0 184#endif 185 186 100 continue 187 endif 188c 189 return 190 end 191#ifndef SECOND_DERIV 192#define SECOND_DERIV 193c 194c Compile source again for the 2nd derivative case 195c 196#include "xc_xpw86.F" 197#endif 198