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