1#if !defined SECOND_DERIV 2 subroutine xc_lb94_e(tol_rho, xfac, rho, delrho, Amat, nq, 3 I ipol, Ex, qwght,ldew,func,delrhoq) 4#elif defined(SECOND_DERIV) 5c For locations of 2nd derivatives of functionals in array 6 subroutine xc_lb94_e_d2(tol_rho, xfac, rho, delrho, Amat, 7 I Amat2, nq, ipol, Ex, qwght,ldew,func,delrhoq) 8#include "dft2drv.fh" 9 implicit none 10#endif 11c 12 double precision tol_rho, xfac 13 double precision Ex 14 integer nq, ipol 15 logical ldew 16 double precision func(*),qwght(*) 17 double precision rho(nq,*) 18 double precision delrho(nq,3,*) 19 double precision Amat(nq,*) 20#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 21c 22c Partial Second Derivatives of the Exchange Energy Functional 23c 24 double precision Amat2(nq,*) 25#endif 26 double precision delrhoq(nq,*) ! [in] 3*rho+ddot(delrho,qxyz) 27c 28 integer n,ii 29 double precision ex_lda ! ignored 30c 31 32 ex_lda=0d0 33#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 34 call xc_dirac(tol_rho, 1d0, .true.,.false.,rho, 35 & amat, nq, ipol,ex_lda,qwght,.false.,func) 36 37 call xc_lb94(tol_rho, 0d0, rho, delrho, 38 & amat, nq, ipol) 39#endif 40#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 41 call xc_dirac_d2(tol_rho, 1d0, .true., .false., rho, Amat, 42 & Amat2, nq, ipol, ex_lda, qwght, .false., 43 & func) 44 45 call xc_lb94_d2(tol_rho, 0d0, rho, delrho, 46 & amat, amat2, nq, ipol) 47#endif 48 49c 50c eq. 12 of DOI:10.1080/00268979600100011 51c or 52c eq. 35 of DOI:10.1103/PhysRevA.51.170 53c delrhoq is computed in xc_delrhodotr 54c 55#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 56 do ii=1,ipol 57 do n=1,nq 58 ex=ex+delrhoq(n,ii)*amat(n,ii)*qwght(n) 59 enddo 60 enddo 61#endif 62 return 63 end 64#if !defined SECOND_DERIV 65 Subroutine xc_lb94(tol_rho, fac, rho, delrho, 66 & Amat, nq, ipol) 67#elif defined(SECOND_DERIV) 68 Subroutine xc_lb94_d2(tol_rho, fac, rho, delrho, 69 & Amat, Amat2,nq, ipol) 70#include "dft2drv.fh" 71#endif 72c 73C$Id$ 74c 75 implicit none 76c 77c 78 double precision tol_rho, fac 79 integer nq, ipol 80c 81c Charge Density 82c 83 double precision rho(nq,*) 84c 85c Charge Density Gradient 86c 87 double precision delrho(nq,3,*) 88c 89c Sampling Matrices for the XC Potential 90c 91 double precision Amat(nq,*) 92c 93#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 94c 95c Partial Second Derivatives of the Exchange Energy Functional 96c 97 double precision Amat2(nq,*) 98#endif 99c 100 double precision BETA 101 Parameter (BETA = 0.05D0) 102c 103c References: 104c 105c R. van Leeuwen & E. J. Baerends, Phys. Rev. A 49, 2421 (1994). 106c 107c*************************************************************************** 108c 109 integer n 110 double precision arcsinh,darcsinh 111 double precision hrho 112 double precision rho13, rho43, gamma, x, g, gdenom 113c double precisiobn rhom23 114#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 115 double precision dg,dgdenom 116#endif 117c 118 arcsinh(x)=log(x+dsqrt(1d0+x*x)) 119 darcsinh(x)=1d0/dsqrt(1d0+x*x) 120c 121 if (ipol.eq.1) then 122c 123c ======> SPIN-RESTRICTED <====== 124c 125 do 10 n = 1, nq 126 if (rho(n,1).lt.tol_rho) goto 10 127c 128c Spin alpha: 129c 130 hrho = 0.5d0*rho(n,1) 131 rho13 = hrho**(1.d0/3.d0) 132 rho43 = rho13*hrho 133c rhom23=rho13/hrho 134 gamma = delrho(n,1,1)*delrho(n,1,1) + 135 & delrho(n,2,1)*delrho(n,2,1) + 136 & delrho(n,3,1)*delrho(n,3,1) 137 if (dsqrt(gamma).gt.tol_rho)then 138 gamma = 0.25d0 * gamma 139 x = dsqrt(gamma) / rho43 140 else 141 x = 0d0 142 endif 143c 144 gdenom = 1d0 + 3d0*BETA*x*arcsinh(x) 145 g = -BETA*x*x / gdenom 146c 147 Amat(n,1) = Amat(n,1) + rho13*g*(1d0-fac) 148#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 149 dgdenom = 3d0*BETA*(arcsinh(x) + x*darcsinh(x)) 150 dg = BETA*x*(x*dgdenom - 2d0*gdenom) / 151 G (gdenom*gdenom) 152 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 153 R BETA/(3d0*rho13*rho13)*(4d0*x*dg - g) 154#endif 155c 156 10 continue 157c 158 else 159c 160c ======> SPIN-UNRESTRICTED <====== 161c 162 do 20 n = 1, nq 163 if (rho(n,1).lt.tol_rho) goto 20 164 if (rho(n,2).lt.tol_rho) goto 25 165c 166c Spin alpha: 167c 168 rho13 = rho(n,2)**(1.d0/3.d0) 169 rho43 = rho13*rho(n,2) 170 gamma = delrho(n,1,1)*delrho(n,1,1) + 171 & delrho(n,2,1)*delrho(n,2,1) + 172 & delrho(n,3,1)*delrho(n,3,1) 173 if (dsqrt(gamma).gt.tol_rho)then 174 x = dsqrt(gamma) / rho43 175 else 176 x = 0d0 177 endif 178c 179 gdenom = 1d0 + 3d0*BETA*x*arcsinh(x) 180 g = -BETA*x*x / gdenom 181c 182 Amat(n,1) = Amat(n,1) + rho13*g*(1d0-fac) 183c 184#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 185 dgdenom = 3d0*BETA*(arcsinh(x) + x*darcsinh(x)) 186 dg = BETA*x*(x*dgdenom - 2d0*gdenom) / 187 G (gdenom*gdenom) 188 Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + 189 R BETA/(3d0*rho13*rho13)*(4d0*x*dg - g) 190#endif 191 25 continue 192c 193c Spin beta: 194c 195 if (rho(n,3).lt.tol_rho) goto 20 196c 197 rho13 = rho(n,3)**(1.d0/3.d0) 198 rho43 = rho13*rho(n,3) 199 gamma = delrho(n,1,2)*delrho(n,1,2) + 200 & delrho(n,2,2)*delrho(n,2,2) + 201 & delrho(n,3,2)*delrho(n,3,2) 202 if (dsqrt(gamma).gt.tol_rho)then 203 x = dsqrt(gamma) / rho43 204 else 205 x = 0d0 206 endif 207c 208 gdenom = 1d0 + 3d0*BETA*x*arcsinh(x) 209 g = -BETA*x*x / gdenom 210c 211 Amat(n,2) = Amat(n,2) + rho13*g*(1d0-fac) 212#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 213 dgdenom = 3d0*BETA*(arcsinh(x) + x*darcsinh(x)) 214 dg = BETA*x*(x*dgdenom - 2d0*gdenom) / 215 G (gdenom*gdenom) 216 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + 217 R BETA/(3d0*rho13*rho13)*(4d0*x*dg - g) 218#endif 219c 220 20 continue 221c 222 endif 223c 224 return 225 end 226#ifndef SECOND_DERIV 227#define SECOND_DERIV 228c 229c Compile source again for the 2nd derivative case 230c 231#include "xc_lb94.F" 232#endif 233