1cc AJL/Begin/FDE 2 Subroutine ts_tf(tol_rho, fac, lfac, nlfac, rho, Amat, nq, 3 & ipol, Ts, qwght, ldew, func) 4c 5 Implicit none 6#include "errquit.fh" 7c 8#include "stdio.fh" 9c 10 integer nq, ipol 11 double precision fac, Ts 12 logical ldew, lfac, nlfac 13 double precision func(*) ! value of the functional [output] 14 double precision two_rho32, two_rho33 15c 16c Charge Density 17c 18 double precision rho(nq,(ipol*(ipol+1))/2) 19c 20c Quadrature Weights 21c 22 double precision qwght(nq) 23c 24c Partial First Derivatives of the Kinetic Energy Functional 25c 26 double precision Amat(nq,ipol) 27c 28c Compute the partial derivatives of the Kinetic functional of Dirac. 29c 30 double precision P1, P2, P3, P4, tol_rho 31c 32c Daniel (5-23-12): P1 and P2 are for restricted calculations, P3 33c and P4 are for unrestricted calculations. 34c 35cc AJL: I made an error in these definitions. Updating on 23/09/2015 36cc P1 = (1/2) *(3*PI^2)**(1/3) 37cc P2 = (3/10) *(3*PI^2)**(1/3) 38cc P3 = 2**(2/3)* (1/2) *(3*PI^2)**(1/3) 39cc P4 = 2**(2/3)* (3/10) *(3*PI^2)**(1/3) 40cc 41cc Calculated 14/10/2014 with Wolfram Alpha 42c Parameter (P1 = 0.15468338631400680D+01) 43c Parameter (P2 = 0.92810031788404078D+00) 44c Parameter (P3 = 0.24554457015685778D+01) 45c Parameter (P4 = 0.14732674209411467D+01) 46 47c New values 48c P1 = (1/2) *(3*PI^2)**(2/3) 49c P2 = (3/10) *(3*PI^2)**(2/3) 50c P3 = 2**(2/3)* (1/2) *(3*PI^2)**(2/3) 51c P4 = 2**(2/3)* (3/10) *(3*PI^2)**(2/3) 52c 53c Calculated 23/09/2015 with Wolfram Alpha 54 Parameter (P1 = 0.47853900003136530D+01) 55 Parameter (P2 = 0.28712340001881918D+01) 56 Parameter (P3 = 0.75963331205759952D+01) 57 Parameter (P4 = 0.45577998723455971D+01) 58 double precision rho13, rhom23, two_third, term 59 double precision rho32, rho33 60 Parameter (two_third = 2.d0/3.d0) 61 integer n 62c 63 term = 0.d0 64c 65 if (ipol.eq.1)then 66c 67c ======> SPIN-RESTRICTED <====== 68c 69c Tested against values from Jones and Gunnarsson, 70c doi: 10.1103/RevModPhys.61.689 71 72c Reported Ts for He: 2.5605 (No basis details, HF orbitals) 73c Our calculated value: 2.5600 (6-311+G, uncontracted, xc = HFexch) 74c 75c Reported Ts for Ne: 117.78 76c Our calculated value: 117.7226 (6-311+G, uncontracted, and xc = HFexch) 77c 78c For spin-polarised implementation we got identical results. All good. 79c 80 do n = 1, nq 81c write(luout,*) rho(n,1) 82 if (rho(n,1).gt.tol_rho)then 83 rho13=rho(n,1)**two_third 84 Amat(n,1) = Amat(n,1) + rho13 85 term = term + rho(n,1)*rho13*qwght(n) 86 if(ldew) func(n) = func(n) + rho(n,1)*rho13*fac*P2 87 endif 88 enddo 89 call dscal(nq, P1*fac, Amat(1,1), 1) 90c 91 Ts = Ts + P2*term*fac 92c 93 else 94c 95c ======> SPIN-UNRESTRICTED <====== 96c 97c The only useful initial test for this I could find was for a hydrogen atom, 98c as listed in Parr & Yang's book "Density Functional Theory of Atoms and 99c Molecules". On page 175 they show the analytical values for Ts for a 100c hydrogen atom as being 0.2891 in spin-paired (with 0.5 e occupancy in 101c each channel) and 0.4590 for the spin-unpaired (1 e in alpha channel). 102c Due to the implementation of spin-paired calculations in NWChem I can't 103c test the former, but for the latter the below implementation results in 104c a value of 0.4540 (xc = LDA). Of course, our result is dependent on this 105c functional but the similarity of the results leaves me in good 106c confidence that this has been correctly implemented now. 107c 108c Additional note for spin-unpaired (23/09/2015): 109c I have found I can also reference against the work of Iyengar, 110c Ernzerhof, Maximoff and Scuseria, DOI: 10.1103/PhysRevA.63.052508 111c Here they list the energies for a whole range of gas phase atoms, with 112c KS orbitals and densities obtained with B88-PW91. For H, they achieved 113c Ts of 0.456, compared to 0.4631 (6-311+G uncontracted, same basis as 114c them) for our implementations and so I am confident this is correct. 115c 116c For an additional test, I did C(triplet) with the same basis and XC. 117c In the published case, the value of Ts is 33.983, whereas for our 118c implementation we got 34.0596. Again, I can only believe the 119c difference is in the numerical implementation that is used in 120c solving the KS equations and the accuracy of the settings for each 121c type of calculation. 122c 123 do n = 1,nq 124c write(luout,*) rho(n,1), rho(n,2), rho(n,3) 125 if (rho(n,1).gt.tol_rho)then 126 127 rho32 = 0.d0 128 rho33 = 0.d0 129 130 if (rho(n,2).gt.tol_rho) rho32 = rho(n,2)**two_third 131 if (rho(n,3).gt.tol_rho) rho33 = rho(n,3)**two_third 132c Debug: Testing the functionality 133c rho32 = rho(n,2)**two_third 134c rho33 = rho(n,3)**two_third 135 136 Amat(n,1) = Amat(n,1) + P3*rho32*fac 137 Amat(n,2) = Amat(n,2) + P3*rho33*fac 138 139 if (ldew)func(n) = func(n) + ((rho(n,2)*rho32) 140 & + (rho(n,3)*rho33))*P4*fac 141 142 term = term+((rho(n,2)*rho32) 143 & + (rho(n,3)*rho33))*qwght(n) 144 145c write(LuOut,*)'from ts_tf ; rho32, rho33, qwght(n): ', 146c & rho32, rho33, qwght(n) 147 endif 148 enddo 149c 150 Ts = Ts + P4*term*fac 151 endif 152c write(LuOut,*)'from ts_tf ; p4, term, fac, Ex: ', 153c & P4, term, fac, Ts 154 return 155 end 156c 157cc AJL/End 158