1cc AJL/Begin/FDE 2c 3 Subroutine ts_vw(tol_rho, fac, lfac, nlfac, rho, delrho, 4 & Amat, Cmat, nq, ipol, Ts, qwght, ldew, func) 5c 6 Implicit none 7#include "errquit.fh" 8c 9#include "stdio.fh" 10c 11#include "dft2drv.fh" 12c#include "dft3drv.fh" 13c 14 integer nq, ipol 15 double precision fac, Ts 16 logical ldew, lfac, nlfac 17 double precision func(*) ! value of the functional [output] 18c 19c Charge Density 20c 21 double precision rho(nq,(ipol*(ipol+1))/2) 22c 23c Charge Density Gradient 24c 25 double precision delrho(nq,3,ipol) 26c 27c Charge Density Laplacian 28c 29c double precision del2rho(nq,ipol*(ipol+1)/2) 30c 31c Quadrature Weights 32c 33 double precision qwght(nq) 34c 35c Partial First Derivatives of the Kinetic Energy Functional 36c 37 double precision Amat(nq,ipol) 38c double precision Cmat(nq,3,ipol) 39 double precision Cmat(nq,*) 40c 41c Compute the partial derivatives of the Kinetic functional of Von Weiszacker 42c 43 double precision tol_rho 44 double precision one_over_rho 45c 46 double precision delrho_2 47 double precision temp_energy, temp_potential, term 48 double precision one_eighth, one_quarter 49 Parameter (one_eighth = 1.d0/8.d0) 50 Parameter (one_quarter = 2.d0*one_eighth) 51 integer n 52c Spin-polarised variables 53 double precision delrho32_2 54 double precision delrho33_2 55 double precision temp_energy1, temp_energy2 56c double precision temp_potential1, temp_potential2 57 58 double precision one_over_rho32, one_over_rho33 59c 60 term = 0.d0 61c 62 if (ipol.eq.1)then 63c 64c ======> SPIN-RESTRICTED <====== 65c 66c Tested against values from Acharya et al. (PNAS, 1980), 67c doi: 10.1073/pnas.77.12.6978 68 69c Reported Ts for He: 2.86 (No basis details, HF orbitals) 70c Our calculated value: 2.8599 (6-311+G, uncontracted, xc = HFexch) 71c 72c Reported Ts for Ne: 90.63 73c Our calculated value: 90.6048 (6-311+G, uncontracted, and xc = HFexch) 74c 75c For spin-polarised implementation we got identical results. All good. 76c 77 78 do 10 n = 1, nq 79c write(luout,*) rho(n,1), delrho(n,1,1), del2rho(n,1) 80 if (rho(n,1).gt.tol_rho)then 81 one_over_rho = 1.d0/rho(n,1) 82 delrho_2= (delrho(n,1,1)*delrho(n,1,1)) 83 & + (delrho(n,2,1)*delrho(n,2,1)) 84 & + (delrho(n,3,1)*delrho(n,3,1)) 85c delrho_2 = max(delrho_2,tol_rho_2) 86c I think this is just in case of numerical instability 87 temp_energy = delrho_2*one_over_rho 88 term = term + temp_energy*qwght(n) 89c 90cc Original derivation of potential, using laplacian if density 91c temp_potential = 2*max(del2rho(n,1),tol_rho) 92c I think this is just in case of numerical instability also! 93c temp_potential = 2*del2rho(n,1) 94c temp_potential = temp_energy - temp_potential 95c \frac{\partial Tw[p]}{\partial \rho} 96c = \frac{\partial tw}{\partial \rho} - 97c \triangledown . \frac{\partial tw}{\partial \triangledown \rho} 98c Amat(n,1) =Amat(n,1) + 99c & one_eighth*temp_potential*one_over_rho*fac 100 101c 102ccc Updated version, using method of Pople and Johnson (1992, CPL) 103cc \frac{\partial tw}{\partial \rho} = 104 Amat(n,1) = Amat(n,1) - 105 & one_eighth*temp_energy*one_over_rho*fac 106cc \frac{\partial tw}{\partial |\triangledown \rho|^{2}} = \frac{\partial tw}{\partial\gamma^{\alpha\alpha}} 107cc = \frac{1}{8 \rho} 108 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 109cc A factor of two is missing when converted from being \partial 110cc \gamma^{\alpha\alpha} to \partial \triangledown \rho. 111cc Therefore we account this here in to our numbers, and it seems to 112cc produce the correct results. 113c & one_eighth*one_over_rho*fac 114 & one_quarter*one_over_rho*fac 115 endif 116 10 continue 117c 118 Ts = Ts + one_eighth*term*fac 119 120 else 121c 122c ======> SPIN-UNRESTRICTED <====== 123c 124c Tested against values from Iyengar, 125c Ernzerhof, Maximoff and Scuseria, DOI: 10.1103/PhysRevA.63.052508, 126 127c Reported Ts for H: 0.500 (6-311+G, uncontracted, xc = B88-PW91) 128c Our calculated value: 0.5088 (6-311+G, uncontracted, xc = B88-PW91) 129c 130c Reported Ts for C: 32.197 (triplet, same setup as above) 131c Our calculated value: 34.0596 (triplet, same setup as above). 132c 133c As with the TF functional, I put these differences down to numerical 134c implementation seeing as the cited reference is ~15 years old. 135c 136c For MPI implementation we got identical results. All good. 137c 138 do 20 n = 1, nq 139c write(luout,*) rho(n,1), delrho(n,1,1)+delrho(n,1,2), 140c & del2rho(n,1), del2rho(n,2), del2rho(n,3) 141 if (rho(n,1).gt.tol_rho) then 142 143 temp_energy1 = 0.0d0 144 temp_energy2 = 0.0d0 145 146 one_over_rho32 = 1.d0/rho(n,2) 147 one_over_rho33 = 1.d0/rho(n,3) 148 149 delrho32_2 = (delrho(n,1,1)*delrho(n,1,1)) 150 & + (delrho(n,2,1)*delrho(n,2,1)) 151 & + (delrho(n,3,1)*delrho(n,3,1)) 152 153 delrho33_2 = (delrho(n,1,2)*delrho(n,1,2)) 154 & + (delrho(n,2,2)*delrho(n,2,2)) 155 & + (delrho(n,3,2)*delrho(n,3,2)) 156 157 if (rho(n,2).gt.tol_rho) 158 & temp_energy1 = delrho32_2*one_over_rho32 159 if (rho(n,3).gt.tol_rho) 160 & temp_energy2 = delrho33_2*one_over_rho33 161 temp_energy = temp_energy1 + temp_energy2 162 term = term + temp_energy*qwght(n) 163 164c temp_potential1 = 0.5d0*temp_energy1 165c temp_potential2 = 0.5d0*temp_energy2 166c 167c Needs correcting to use Cmat 168c temp_potential1 = temp_potential1-del2rho(n,1) 169c temp_potential2 = temp_potential2-del2rho(n,2) 170 171 if (rho(n,2).gt.tol_rho) then 172 Amat(n,1) =Amat(n,1) - 173c & one_quarter*temp_potential1*one_over_rho32*fac 174 & one_eighth*temp_energy1*one_over_rho32*fac 175 Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + 176 & one_eighth*one_over_rho32*fac 177 endif 178c 179 if (rho(n,3).gt.tol_rho) then 180 Amat(n,2) =Amat(n,2) - 181c & one_quarter*temp_potential2*one_over_rho33*fac 182 & one_eighth*temp_energy2*one_over_rho33*fac 183 Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + 184 & one_eighth*one_over_rho33*fac 185 endif 186c 187 endif 188 20 continue 189c 190 Ts = Ts + one_eighth*term*fac 191 192 endif 193c 194 return 195 end 196c 197cc AJL/End 198