1c 2c LDA screened exchange functional 3c 4c BGJ - 2/01 5c 6c Reference: 7c 8c D. M. Bylander and L. Kleinman, Phys. Rev. B 41, 7868 (1990). 9c 10#ifndef SECOND_DERIV 11 Subroutine xc_dirac_screened(tol_rho, fac, lfac, nlfac, rho, 12 & Amat, nq, ipol, Ex, Ks, qwght, ldew, func) 13#else 14#include "dft2drv.fh" 15 Subroutine xc_dirac_screened_d2(tol_rho, fac, lfac, nlfac, rho, 16 & Amat, Amat2, nq, ipol, Ex, Ks, qwght, ldew, func) 17#endif 18c 19C$Id$ 20c 21 implicit none 22#include "errquit.fh" 23c 24c 25 double precision tol_rho, fac, Ex 26 integer nq, ipol 27 logical lfac, nlfac,ldew 28 double precision Ks ! Thomas-Fermi screening constant [input] 29 double precision func(*) ! value of the functional [output] 30c 31c Charge Density 32c 33 double precision rho(nq,ipol*(ipol+1)/2) 34c 35c Quadrature Weights 36c 37 double precision qwght(nq) 38c 39c Sampling Matrices for the XC Potential 40c 41 double precision Amat(nq,ipol) 42#ifdef SECOND_DERIV 43 double precision Amat2(nq,NCOL_AMAT2) 44#endif 45c 46 double precision ONE3, FOUR3 47 Parameter (ONE3 = 1d0/3d0, FOUR3 = 4d0/3d0) 48 integer n 49 double precision pi, C, TPP, rhoval, rho13, rho43, z, d1z, 50 & f, d1f 51 double precision kf, d1kf, f1, f2, d1f1, d1f2 52 double precision rz, Tworz, atan2rz, datan2rz, rdatan, lnrdatan, 53 & z243, Eightrz2, Eightrz3, dzdp 54#ifdef SECOND_DERIV 55 double precision d2z, d2f 56c 57 call errquit('Must implement 2nd derivs in xc_dirac_screened',0, 58 & CAPMIS_ERR) 59#endif 60c 61 pi = acos(-1d0) 62 C = -1.5d0*(0.75d0/pi)**ONE3 63 TPP = (3.d0*pi*pi)**ONE3 64c 65 if (ipol.eq.1) then 66c 67c ======> SPIN-RESTRICTED <====== 68c 69 do 10 n = 1, nq 70 if (rho(n,1).lt.tol_rho) goto 10 71 rhoval = 0.5d0*rho(n,1) 72c 73c Spin alpha: 74c 75 rho13 = rhoval**ONE3 76 rho43 = rho13*rhoval 77c 78 kf = TPP*rho13 79 z = Ks/kf 80 rz = 1d0/z 81 Tworz = 2d0*rz 82 Eightrz2 = 4d0*Tworz*rz 83 Eightrz3 = Eightrz2*rz 84 rdatan = 1d0 + Tworz*Tworz 85 datan2rz = 1d0/rdatan 86 lnrdatan = log(rdatan) 87 atan2rz = atan(Tworz) 88 z243 = z*z/4d0 + 3d0 89 f1 = -z*z/6d0 90 f2 = 1d0 - z243*lnrdatan + 8d0*rz*atan2rz 91 d1f1 = -z/3d0 92 d1f2 = -z*0.5d0*lnrdatan - Eightrz2*atan2rz 93 & - Eightrz3*datan2rz*(2d0 - z243) 94 f = 1d0 + f1*f2 95 dzdp = -ONE3*z/rhoval 96 d1f = (d1f1*f2 + f1*d1f2)*dzdp 97c 98 if (lfac) then 99 Ex = Ex + 2d0*rho43*C*f*qwght(n)*fac 100 if(ldew)func(n) = func(n) + 2.d0*rho43*C*f*fac 101 Amat(n,1) = Amat(n,1) + (FOUR3*rho13*f+rho43*d1f)*C*fac 102#ifdef SECOND_DERIV 103 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 104 & + ONE3*(FOUR3*rho13/rhoval)*C*fac 105#endif 106 endif 107c 108 10 continue 109c 110 else 111c 112c ======> SPIN-UNRESTRICTED <====== 113c 114 do 20 n = 1, nq 115 if (rho(n,1).lt.tol_rho) goto 20 116 if (rho(n,2).lt.tol_rho) goto 25 117c 118c Spin alpha: 119c 120 rhoval = rho(n,2) 121 rho13 = rhoval**ONE3 122 rho43 = rho13*rhoval 123c 124 kf = TPP*rho13 125 z = Ks/kf 126 rz = 1d0/z 127 Tworz = 2d0*rz 128 Eightrz2 = 4d0*Tworz*rz 129 Eightrz3 = Eightrz2*rz 130 rdatan = 1d0 + Tworz*Tworz 131 datan2rz = 1d0/rdatan 132 lnrdatan = log(rdatan) 133 atan2rz = atan(Tworz) 134 z243 = z*z/4d0 + 3d0 135 f1 = -z*z/6d0 136 f2 = 1d0 - z243*lnrdatan + 8d0*rz*atan2rz 137 d1f1 = -z/3d0 138 d1f2 = -z*0.5d0*lnrdatan - Eightrz2*atan2rz 139 & - Eightrz3*datan2rz*(2d0 - z243) 140 f = 1d0 + f1*f2 141 dzdp = -ONE3*z/rhoval 142 d1f = (d1f1*f2 + f1*d1f2)*dzdp 143c 144 if (lfac) then 145 Ex = Ex + rho43*C*f*qwght(n)*fac 146 if(ldew)func(n) = func(n) + rho43*C*f*fac 147 Amat(n,1) = Amat(n,1) + (FOUR3*rho13*f+rho43*d1f)*C*fac 148#ifdef SECOND_DERIV 149 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) 150 & + ONE3*(FOUR3*rho13/rhoval)*C*fac 151#endif 152 endif 153c 154 25 continue 155c 156c Spin beta: 157c 158 if (rho(n,3).lt.tol_rho) goto 20 159c 160 rhoval = rho(n,3) 161 rho13 = rhoval**ONE3 162 rho43 = rho13*rhoval 163c 164 kf = TPP*rho13 165 z = Ks/kf 166 rz = 1d0/z 167 Tworz = 2d0*rz 168 Eightrz2 = 4d0*Tworz*rz 169 Eightrz3 = Eightrz2*rz 170 rdatan = 1d0 + Tworz*Tworz 171 datan2rz = 1d0/rdatan 172 lnrdatan = log(rdatan) 173 atan2rz = atan(Tworz) 174 z243 = z*z/4d0 + 3d0 175 f1 = -z*z/6d0 176 f2 = 1d0 - z243*lnrdatan + 8d0*rz*atan2rz 177c f2 = 1d0 - z243*log(1d0+rz) + 8d0*rz*atan2rz 178 d1f1 = -z/3d0 179 d1f2 = -z*0.5d0*lnrdatan - Eightrz2*atan2rz 180 & - Eightrz3*datan2rz*(2d0 - z243) 181 f = 1d0 + f1*f2 182 dzdp = -ONE3*z/rhoval 183 d1f = (d1f1*f2 + f1*d1f2)*dzdp 184c 185 if (lfac) then 186 Ex = Ex + rho43*C*f*qwght(n)*fac 187 if(ldew)func(n) = func(n) + rho43*C*f*fac 188 Amat(n,2) = Amat(n,2) + (FOUR3*rho13*f+rho43*d1f)*C*fac 189#ifdef SECOND_DERIV 190 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) 191 & + ONE3*(FOUR3*rho13/rhoval)*C*fac 192#endif 193 endif 194c 195 20 continue 196c 197 endif 198c 199 return 200 end 201#ifndef SECOND_DERIV 202#define SECOND_DERIV 203c 204c Compile source again for the 2nd derivative case 205c 206#include "xc_dirac_screened.F" 207#endif 208