1c 2C$Id$ 3c 4 subroutine xc_optx(tol_rho, fac, p, rho, delrho, 5 & Amat, Cmat, nq, ipol, Ex, qwght,ldew,func) 6 implicit none 7c 8#include "dft2drv.fh" 9c 10 double precision tol_rho, fac, Ex 11 integer nq, ipol 12 logical ldew 13 double precision func(*) ! value of the functional [output] 14c 15c Charge Density 16c 17 double precision rho(nq,ipol*(ipol+1)/2) 18c 19c Charge Density Gradient 20c 21 double precision delrho(nq,3,ipol) 22c 23c Quadrature Weights 24c 25 double precision qwght(nq) 26c 27c Sampling Matrices for the XC Potential 28c 29 double precision Amat(nq,ipol), Cmat(nq,*) 30c 31 integer p ! [in] 32c 33c 34c References: 35c 36c Becke, (1986) 37c Handy NC, Cohen AJ, Mol Phys 99 (5): 403-412 MAR 2001 38c idem, Mol Phys 99 (7); 607-615 2001 39c 40c*************************************************************************** 41c 42 integer n 43 double precision rho13, rho43, gamma, x, g, dg, 44 & t, hrho 45 double precision gamma86 46 Parameter (gamma86=0.006d0) 47 double precision ux,uxp,gx 48 ux(x,gx)=gx*x*x/(1d0+gx*x*x) 49 uxp(x,gx)=gx*x*2d0/(1d0+gx*x*x)**2 50c 51 if (ipol.eq.1) then 52c 53c ======> SPIN-RESTRICTED <====== 54c 55 do 10 n = 1, nq 56 if (rho(n,1).lt.tol_rho) goto 10 57c 58c Spin alpha: 59c 60 hrho = 0.5d0*rho(n,1) 61 rho13 = hrho**(1.d0/3.d0) 62 rho43 = rho13*hrho 63 gamma = delrho(n,1,1)*delrho(n,1,1) + 64 & delrho(n,2,1)*delrho(n,2,1) + 65 & delrho(n,3,1)*delrho(n,3,1) 66 if (dsqrt(gamma).gt.tol_rho)then 67 gamma = 0.25d0 * gamma 68 x = sqrt(gamma) / rho43 69 else 70 x = 0d0 71 endif 72c 73 g = -ux(x,gamma86)**p 74 dg = -p*ux(x,gamma86)**(p-1)*uxp(x,gamma86) 75c 76c 77 Ex = Ex + 2d0*rho43*g*qwght(n)*fac 78 if(ldew)func(n) = func(n) + 2.d0*rho43*g*fac 79 Amat(n,1) = Amat(n,1) + (4d0/3d0)*rho13*(g-x*dg)*fac 80c 81 if (x.gt.tol_rho) then 82 t = 0.5d0 * dg / sqrt(gamma) * fac 83 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + t 84 endif 85c 86c 87 10 continue 88c 89 else 90c 91c ======> SPIN-UNRESTRICTED <====== 92c 93 do 20 n = 1, nq 94 if (rho(n,1).lt.tol_rho) goto 20 95 if (rho(n,2).lt.tol_rho) goto 25 96c 97c Spin alpha: 98c 99 rho13 = rho(n,2)**(1.d0/3.d0) 100 rho43 = rho13*rho(n,2) 101 gamma = delrho(n,1,1)*delrho(n,1,1) + 102 & delrho(n,2,1)*delrho(n,2,1) + 103 & delrho(n,3,1)*delrho(n,3,1) 104 if (dsqrt(gamma).gt.tol_rho)then 105 x = dsqrt(gamma) / rho43 106 else 107 x = 0d0 108 endif 109c 110 g = -ux(x,gamma86)**p 111 dg = -p*ux(x,gamma86)**(p-1)*uxp(x,gamma86) 112c 113c 114 Ex = Ex + rho43*g*qwght(n)*fac 115 if (ldew)func(n) = func(n) + rho43*g*fac 116 Amat(n,1) = Amat(n,1) + (4d0/3d0)*rho13*(g-x*dg)*fac 117c 118 if (x.gt.tol_rho) then 119 t = dg / sqrt(gamma) * fac 120 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + t * 0.5d0 121 endif 122c 123c 124 25 continue 125c 126c Spin beta: 127c 128 if (rho(n,3).lt.tol_rho) goto 20 129c 130 rho13 = rho(n,3)**(1.d0/3.d0) 131 rho43 = rho13*rho(n,3) 132 gamma = delrho(n,1,2)*delrho(n,1,2) + 133 & delrho(n,2,2)*delrho(n,2,2) + 134 & delrho(n,3,2)*delrho(n,3,2) 135 if (dsqrt(gamma).gt.tol_rho)then 136 x = dsqrt(gamma) / rho43 137 else 138 x = 0d0 139 endif 140 g = -ux(x,gamma86)**p 141 dg = -p*ux(x,gamma86)**(p-1)*uxp(x,gamma86) 142c 143c 144 145 Ex = Ex + rho43*g*qwght(n)*fac 146 if (ldew)func(n) = func(n) +rho43*g*fac 147 Amat(n,2) = Amat(n,2) + (4d0/3d0)*rho13*(g-x*dg)*fac 148c 149 if (x.gt.tol_rho) then 150 t = dg / sqrt(gamma) * fac 151 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + t * 0.5d0 152 endif 153c 154c 155 20 continue 156c 157 endif 158c 159 return 160 end 161