1c $Id: xc_eval_fnl.F 27518 2015-09-17 09:34:57Z alogsdail $ 2 subroutine xc_xn12(tol_rho, fac,lfac,nlfac, rho, delrho, 3 & amat, cmat, nq, ipol, ex, 4 & qwght, ldew, func) 5 implicit none 6************************************************************************ 7* * 8* N12 evaluates the exchange part of the N12 and N12-SX * 9* functionals on the grid. * 10* * 11* OUTPUT: * 12* F - Functional values * 13* D1F - First derivatives with respect to RA, RB, GA, GB * 14* * 15* INPUT: * 16* ijzy - 1 N12 * 17* ijzy - 2 N12-SX * 18* * 19* RA,B - Spin densities * 20* D1RA,B - Spin density gradients * 21* NGrid - number of grid points * 22* * 23* Analytic second derivatives are not available yet!!! * 24* * 25* RP (09/12) * 26* * 27************************************************************************ 28#include "errquit.fh" 29#include "dft2drv.fh" 30c 31 double precision fac, Ex 32 integer nq, ipol 33 logical lfac, nlfac,ldew, uselc 34 double precision func(*) ! value of the functional [output] 35c 36c Charge Density & Its Cube Root 37c 38 double precision rho(nq,ipol*(ipol+1)/2) 39c 40c Charge Density Gradient 41c 42 double precision delrho(nq,3,ipol) 43c 44c Quadrature Weights 45c 46 double precision qwght(nq) 47c 48c Sampling Matrices for the XC Potential & Energy 49c 50 double precision Amat(nq,ipol), Cmat(nq,*) 51 double precision tol_rho 52 external xc_xn12_os 53 call xc_os2cs(xc_xn12_os, 54 & tol_rho, fac,lfac,nlfac, rho, delrho, 55 & amat, cmat, nq, ipol, ex, 56 & qwght, ldew, func) 57 return 58 end 59 subroutine xc_xn12_os(tol_rho, fac,lfac,nlfac, rho, delrho, 60 & amat, cmat, nq, ex, 61 & qwght, ldew, func, fact_cs) 62 Implicit Real*8(A-H,O-Z) 63#include "errquit.fh" 64#include "dft2drv.fh" 65c 66c wrapper for n12 open-shell 67 double precision fac, Ex 68 double precision fact_cs 69 integer nq, ipol 70 logical lfac, nlfac,ldew, uselc 71 double precision func(*) ! value of the functional [output] 72c 73c Charge Density & Its Cube Root 74c 75 double precision rho(nq) 76c 77c Charge Density Gradient 78c 79 double precision delrho(nq,3) 80c 81c Quadrature Weights 82c 83 double precision qwght(nq) 84c 85c Sampling Matrices for the XC Potential & Energy 86c 87 double precision Amat(nq), Cmat(nq) 88 double precision tol_rho 89 90 call xc_xn12_0(tol_rho, fac,lfac,nlfac, rho, delrho, 91 & amat, cmat, nq, ex, 92 & qwght, ldew, func, 1,fact_cs) 93 return 94 end 95 subroutine xc_xn12sx(tol_rho, fac,lfac,nlfac, rho, delrho, 96 & amat, cmat, nq, ipol, ex, 97 & qwght, ldew, func) 98 implicit none 99#include "errquit.fh" 100#include "dft2drv.fh" 101c 102 double precision fac, Ex 103 integer nq, ipol 104 logical lfac, nlfac,ldew, uselc 105 double precision func(*) ! value of the functional [output] 106 double precision rho(nq,ipol*(ipol+1)/2) 107 double precision delrho(nq,3,ipol) 108 double precision qwght(nq) 109 double precision Amat(nq,ipol), Cmat(nq,*) 110 double precision tol_rho 111 external xc_xn12sx_os 112 call xc_os2cs(xc_xn12sx_os, 113 & tol_rho, fac,lfac,nlfac, rho, delrho, 114 & amat, cmat, nq, ipol, ex, 115 & qwght, ldew, func) 116 return 117 end 118 subroutine xc_xn12sx_os(tol_rho, fac,lfac,nlfac, rho, delrho, 119 & amat, cmat, nq, ex, 120 & qwght, ldew, func, fact_cs) 121 Implicit Real*8(A-H,O-Z) 122#include "errquit.fh" 123#include "dft2drv.fh" 124c 125c wrapper for n12 open-shell 126 double precision fac, Ex 127 double precision fact_cs 128 integer nq, ipol 129 logical lfac, nlfac,ldew, uselc 130 double precision func(*) ! value of the functional [output] 131c 132c Charge Density & Its Cube Root 133c 134 double precision rho(nq) 135c 136c Charge Density Gradient 137c 138 double precision delrho(nq,3) 139c 140c Quadrature Weights 141c 142 double precision qwght(nq) 143c 144c Sampling Matrices for the XC Potential & Energy 145c 146 double precision Amat(nq), Cmat(nq) 147 double precision tol_rho 148 149 call xc_xn12_0(tol_rho, fac,lfac,nlfac, rho, delrho, 150 & amat, cmat, nq, ex, 151 & qwght, ldew, func, 2,fact_cs) 152 return 153 end 154 subroutine xc_xn12_0(tol_rho, fac,lfac,nlfac, rho, delrho, 155 & amat, cmat, nq, ex, 156 & qwght, ldew, func, ijzy,fact_cs) 157 Implicit Real*8(A-H,O-Z) 158#include "errquit.fh" 159#include "dft2drv.fh" 160c 161 integer ijzy 162 double precision fac, Ex 163 double precision fact_cs 164 integer nq, ipol 165 logical lfac, nlfac,ldew, uselc 166 double precision func(*) ! value of the functional [output] 167c 168c Charge Density & Its Cube Root 169c 170 double precision rho(nq) 171c 172c Charge Density Gradient 173c 174 double precision delrho(nq,3) 175c 176c Quadrature Weights 177c 178 double precision qwght(nq) 179c 180c Sampling Matrices for the XC Potential & Energy 181c 182 double precision Amat(nq), Cmat(nq) 183 double precision tol_rho, pi 184c Real*8 RA(*),RB(*),D1RA(NGrid,3),D1RB(NGrid,3),F(*),D1F(NGrid,*) 185 Save One,Two,Three,Four,Five,Six,Seven,Eight,Nine 186 Data One/1.0d0/,Two/2.0d0/,Three/3.0d0/,Four/4.0d0/,Five/5.0d0/, 187 $ Six/6.0d0/,Seven/7.0d0/,Eight/8.0d0/,Nine/9.0d0/ 188C 189C 190 G = 0.004d+0 191 ome = 2.5d+0 192 If (ijzy.eq.1) then 193c 194c N12 195c 196 CC00 = 1.00000D+00 197 CC01 = 5.07880D-01 198 CC02 = 1.68233D-01 199 CC03 = 1.28887D-01 200 CC10 = 8.60211D-02 201 CC11 = -1.71008D+01 202 CC12 = 6.50814D+01 203 CC13 = -7.01726D+01 204 CC20 = -3.90755D-01 205 CC21 = 5.13392D+01 206 CC22 = -1.66220D+02 207 CC23 = 1.42738D+02 208 CC30 = 4.03611D-01 209 CC31 = -3.44631D+01 210 CC32 = 7.61661D+01 211 CC33 = -2.41834D+00 212 Else If (ijzy.eq.2) then 213c 214c N12-SX 215c 216 CC00 = 6.81116D-01 217 CC01 = 1.88858D+00 218 CC02 = 1.78590D+00 219 CC03 = 8.79456D-01 220 CC10 = -8.12270D-02 221 CC11 = -1.08723D+00 222 CC12 = -4.18682D+00 223 CC13 = -3.00000D+01 224 CC20 = 5.36236D-01 225 CC21 = -5.45678D+00 226 CC22 = 3.00000D+01 227 CC23 = 5.51105D+01 228 CC30 = -7.09913D-01 229 CC31 = 1.30001D+01 230 CC32 = -7.24877D+01 231 CC33 = 2.98363D+01 232 End If 233 F12 = Two * Six 234 F24 = Four * Six 235 F28 = Four * Seven 236 F2o3 = Two / Three 237 F3o2 = Three / Two 238 F1o3 = One / Three 239 F4o3 = Four / Three 240 F7o3 = Seven / Three 241 F8o3 = Eight / Three 242 F10o3 = F2o3 * Five 243 F28o9 = F28 / Nine 244 Pi = ACos(-One) 245C 246 Ax = -F3o2*(F4o3*PI)**(-F1o3) 247C 248C Alpha contributions ... 249C 250 Do 10 n = 1, nq 251 If(rho(n).gt.tol_rho) then 252 pX = rho(n) 253 GamX2 =(delrho(n,1)*delrho(n,1) + 254 & delrho(n,2)*delrho(n,2) + 255 & delrho(n,3)*delrho(n,3)) 256 S2 = GamX2*pX**(-F8o3) 257 U = G*S2/(One+G*S2) 258 E = Ax*pX**F4o3 259 FV = ome*pX**F1o3/(One+ome*pX**F1o3) 260c 261 FN12 =CC00 +CC01 *U+CC02 *U**2+CC03 *U**3+ 262 $ CC10*FV +CC11*FV *U+CC12*FV *U**2+CC13*FV *U**3+ 263 $ CC20*FV**2+CC21*FV**2*U+CC22*FV**2*U**2+CC23*FV**2*U**3+ 264 $ CC30*FV**3+CC31*FV**3*U+CC32*FV**3*U**2+CC33*FV**3*U**3 265c 266 ex = ex + E*FN12*qwght(n)*fac*fact_cs 267 if(ldew) func(n)=func(n) + E*FN12*fac*fact_cs 268c 269c First Derivative 270c 271 ER = F4o3*E/pX 272 S = Sqrt(S2) 273 GamX = Sqrt(GamX2) 274 SR = -F4o3*S/pX 275 SG = S/GamX 276 US = Two*G*S/((One+G*S*S)**2) 277c 278 dFVdR = (ome/(Three*pX**F2o3)) 279 $ *(One+ome*pX**F1o3)**(-Two) 280c 281 dFN12dU=CC01+2.0d0*CC02 *U+3.0d0*CC03 *U**2+ 282 $ CC11*FV +2.0d0*CC12*FV *U+3.0d0*CC13*FV *U**2+ 283 $ CC21*FV**2+2.0d0*CC22*FV**2*U+3.0d0*CC23*FV**2*U**2+ 284 $ CC31*FV**3+2.0d0*CC32*FV**3*U+3.0d0*CC33*FV**3*U**2 285 dFN12dV= 286 $ CC10+CC11 *U+CC12 *U**2+CC13 *U**3+ 287 $ 2.0d0*CC20*FV +2.0d0*CC21*FV*U+ 288 $ 2.0d0*CC22*FV*U**2+2.0d0*CC23*FV*U**3+ 289 $ 3.0d0*CC30*FV**2+3.0d0*CC31*FV**2*U+ 290 $ 3.0d0*CC32*FV**2*U**2+3.0d0*CC33*FV**2*U**3 291c 292 dFN12dR = dFN12dU*US*SR+dFN12dV*dFVdR 293 dFN12dG = dFN12dU*US*SG 294c 295 amat(n) = amat(n) + (ER*FN12 296 $ + E*dFN12dR)*fac 297 cmat(n) = cmat(n) + E*dFN12dG/(Two*GamX)*fac 298 endIf 299 10 Continue 300 Return 301 End 302 Subroutine xc_xn12_d2() 303 call errquit(' not coded ',0,0) 304 return 305 end 306