1#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
2#if !defined(NWAD_PRINT)
3C> \ingroup nwxc
4C> @{
5C>
6C> \file nwxc_x_tpss03.F
7C> The TPSS exchange functional
8C>
9C> @}
10#endif
11#endif
12C>
13C> \ingroup nwxc_priv
14C> @{
15C>
16C> \brief Evaluate the TPSS exchange functional
17C>
18C> Evaluate the TPSS meta-GGA exchange functional [1,2].
19C>
20C> Due to the form of the meta-GGAs we need to screen on the kinetic
21C> energy density to ensure that LDA will be obtained when the kinetic
22C> energy density goes to zero [3].
23C>
24C> ### References ###
25C>
26C> [1] J. Tao, J.P. Perdew, V.N. Staveroverov, G.E. Scuseria,
27C> "Climbing the density functional ladder: Nonemperical meta-
28C> generalized gradient approximation designed for molecules
29C> and solids",
30C> Phys. Rev. Lett. <b>91</b>, 146401-146404 (2003), DOI:
31C> <a href="https://doi.org/10.1103/PhysRevLett.91.146401">
32C> 10.1103/PhysRevLett.91.146401</a>.
33C>
34C> [2] J.P. Perdew, J. Tao, V.N. Staveroverov, G.E. Scuseria,
35C> "Meta-generalized gradient approximation: Explanation of a
36C> realistic nonempirical density functional",
37C> J. Chem. Phys. <b>120</b>, 6898-6911 (2004), DOI:
38C> <a href="https://doi.org/10.1063/1.1665298">
39C> 10.1103/1.1665298</a>.
40C>
41C> [3] J. Gr&auml;fenstein, D. Izotov, D. Cremer,
42C>     "Avoiding singularity problems associated with meta-GGA exchange
43C>     and correlation functionals containing the kinetic energy
44C>     density", J. Chem. Phys. <b>127</b>, 214103 (2007), DOI:
45C>     <a href="https://doi.org/10.1063/1.2800011">
46C>     10.1063/1.2800011</a>.
47C>
48c
49c$Id$
50c
51c    Tao,Perdew, Staroverov, Scuseria exchange functional
52c           META GGA
53C         utilizes ingredients:
54c                              rho   -  density
55c                              delrho - gradient of density
56c                              tau - K.S kinetic energy density
57c                              tauW - von Weiszacker kinetic energy density
58c                              tauU - uniform-gas KE density
59c     References:
60c     [a] J. Tao, J.P. Perdew, V.N.Staroverov, G. Scuseria
61c         PRL 91, 146401 (2003).
62c     [b] J. Tao, J.P. Perdew, V.N.Staroverov, G. Scuseria
63c         JCP 120, 6898 (2004).
64#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
65#if defined(NWAD_PRINT)
66      Subroutine nwxc_x_tpss03_p(tol_rho, ipol, nq, wght,
67     &                           rho, rgamma, tau, func)
68#else
69      Subroutine nwxc_x_tpss03(tol_rho, ipol, nq, wght,
70     &                         rho, rgamma, tau, func)
71#endif
72#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
73      Subroutine nwxc_x_tpss03_d2(tol_rho, ipol, nq, wght,
74     &                            rho, rgamma, tau, func)
75#else
76      Subroutine nwxc_x_tpss03_d3(tol_rho, ipol, nq, wght,
77     &                            rho, rgamma, tau, func)
78#endif
79c
80#include "nwad.fh"
81c
82      implicit none
83c
84#include "nwxc_param.fh"
85c
86c     Input and other parameters
87c
88      double precision tol_rho !< [Input] The lower limit on the density
89      integer nq               !< [Input] The number of points
90      integer ipol             !< [Input] The number of spin channels
91      double precision wght    !< [Input] The weight of the functional
92c
93c     Charge Density
94c
95      type(nwad_dble)::rho(nq,*) !< [Input] The density
96c
97c     Charge Density Gradient
98c
99      type(nwad_dble)::rgamma(nq,*) !< [Input] The norm of the density gradients
100c
101c     Kinetic Energy Density
102c
103      type(nwad_dble)::tau(nq,*) !< [Input] The kinetic energy density
104c
105c     The functional
106c
107      type(nwad_dble)::func(*)  !< [Output] The value of the functional
108c
109c     Sampling Matrices for the XC Potential & Energy
110c
111c     double precision Amat(nq,*) !< [Output] The derivative wrt rho
112c     double precision Cmat(nq,*) !< [Output] The derivative wrt rgamma
113c     double precision Mmat(nq,*) !< [Output] The derivative wrt tau
114c
115      integer ispin,cmatpos
116c
117      if (ipol.eq.1 )then
118c
119c     SPIN-RESTRICTED
120c     Ex = Ex[n]
121c
122#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
123#if defined(NWAD_PRINT)
124         call nwxc_x_tpss03_cs_p(1.0d0, tol_rho, ipol, nq, wght,
125     &                      rho, rgamma, tau, func)
126#else
127         call nwxc_x_tpss03_cs(1.0d0, tol_rho, ipol, nq, wght,
128     &                      rho, rgamma, tau, func)
129#endif
130#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
131         call nwxc_x_tpss03_cs_d2(1.0d0, tol_rho, ipol, nq, wght,
132     &                      rho, rgamma, tau, func)
133#else
134         call nwxc_x_tpss03_cs_d3(1.0d0, tol_rho, ipol, nq, wght,
135     &                      rho, rgamma, tau, func)
136#endif
137      else
138c
139c     SPIN-UNRESTRICTED
140c     Ex = Ex[2n_up]/2 + Ex[2n_down]/2
141
142#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
143#if defined(NWAD_PRINT)
144         call nwxc_x_tpss03_cs_p(2.0d0, tol_rho, ipol, nq, wght,
145     &                      rho(1,R_A), rgamma(1,G_AA), tau(1,T_A),
146     &                      func)
147         call nwxc_x_tpss03_cs_p(2.0d0, tol_rho, ipol, nq, wght,
148     &                      rho(1,R_B), rgamma(1,G_BB), tau(1,T_B),
149     &                      func)
150#else
151         call nwxc_x_tpss03_cs(2.0d0, tol_rho, ipol, nq, wght,
152     &                      rho(1,R_A), rgamma(1,G_AA), tau(1,T_A),
153     &                      func)
154         call nwxc_x_tpss03_cs(2.0d0, tol_rho, ipol, nq, wght,
155     &                      rho(1,R_B), rgamma(1,G_BB), tau(1,T_B),
156     &                      func)
157#endif
158#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
159         call nwxc_x_tpss03_cs_d2(2.0d0, tol_rho, ipol, nq, wght,
160     &                      rho(1,R_A), rgamma(1,G_AA), tau(1,T_A),
161     &                      func)
162         call nwxc_x_tpss03_cs_d2(2.0d0, tol_rho, ipol, nq, wght,
163     &                      rho(1,R_B), rgamma(1,G_BB), tau(1,T_B),
164     &                      func)
165#else
166         call nwxc_x_tpss03_cs_d3(2.0d0, tol_rho, ipol, nq, wght,
167     &                      rho(1,R_A), rgamma(1,G_AA), tau(1,T_A),
168     &                      func)
169         call nwxc_x_tpss03_cs_d3(2.0d0, tol_rho, ipol, nq, wght,
170     &                      rho(1,R_B), rgamma(1,G_BB), tau(1,T_B),
171     &                      func)
172#endif
173      endif
174      return
175      end
176#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV)
177#if defined(NWAD_PRINT)
178      Subroutine nwxc_x_tpss03_cs_p(facttwo, tol_rho, ipol, nq, wght,
179     &                      rho, rgamma, tau, func)
180#else
181      Subroutine nwxc_x_tpss03_cs(facttwo, tol_rho, ipol, nq, wght,
182     &                      rho, rgamma, tau, func)
183#endif
184#elif defined(SECOND_DERIV) && !defined(THIRD_DERIV)
185      Subroutine nwxc_x_tpss03_cs_d2(facttwo, tol_rho, ipol, nq, wght,
186     &                      rho, rgamma, tau, func)
187#else
188      Subroutine nwxc_x_tpss03_cs_d3(facttwo, tol_rho, ipol, nq, wght,
189     &                      rho, rgamma, tau, func)
190#endif
191#include "nwad.fh"
192      implicit none
193c
194c     Input and other parameters
195c
196      double precision facttwo !< [Input] Scale factor
197                               !< - 1 for closed shell calculations
198                               !< - 2 for open shell calculations
199      double precision tol_rho !< [Input] The lower limit on the density
200      integer nq               !< [Input] The number of points
201      integer ipol             !< [Input] The number of spin channels
202      double precision wght    !< [Input] The weight of the functional
203c
204c     Charge Density
205c
206      type(nwad_dble)::rho(nq) !< [Input] The density
207c
208c     Charge Density Gradient
209c
210      type(nwad_dble)::rgamma(nq) !< [Input] The norm of the density gradients
211c
212c     Kinetic Energy Density
213c
214      type(nwad_dble)::tau(nq) !< [Input] The kinetic energy density
215c
216c     The functional
217c
218      type(nwad_dble)::func(*)  !< [Output] The value of the functional
219c
220c     Sampling Matrices for the XC Potential & Energy
221c
222c     double precision Amat(nq) !< [Output] The derivative wrt rho
223c     double precision Cmat(nq) !< [Output] The derivative wrt rgamma
224c     double precision Mmat(nq) !< [Output] The derivative wrt tau
225c
226      double precision pi
227      integer n
228      type(nwad_dble)::rrho, rho43, gamma
229      type(nwad_dble)::tauN, tauW, tauU
230
231      type(nwad_dble):: p, qtil, x,  al, mt, z
232      double precision   F83, F23, F53, F43, F13
233      double precision   G920
234      double precision  b,c,e,es
235      double precision    C1, C2, C3
236      double precision    kap, mu
237      type(nwad_dble)::xb,xc,xd
238      type(nwad_dble)::x1,x2,x3,x4,x5,x6,x7
239      double precision   P32, Ax
240c     functional derivatives below FFFFFFFFFFFF
241      double precision dzdn, dpdn, daldn, dqtdn
242      double precision til1, til2
243      double precision dtil2dn, dtil1dn
244      double precision ax1, bx1, dx1dn
245      double precision dx2dn
246      double precision dxddn, dxcdn, dx3dn
247      double precision dx4dn, dx5dn, dx6dn, dx7dn
248      double precision  dxdn
249      double precision xmany, dxmanydn
250      double precision dmtdn, derivn
251
252      double precision dzdg, dpdg, daldg, dqtdg
253      double precision dtil2dg
254      double precision dx1dg, dx2dg
255      double precision dxcdg, dxddg,dx3dg
256      double precision dx4dg, dx5dg, dx6dg, dx7dg
257      double precision dxmanydg, dxdg, dmtdg, derivg
258
259      double precision dzdt, daldt, dqtdt
260      double precision dx1dt, dx2dt, dx3dt
261      double precision dx5dt
262      double precision dxmanydt, dxdt, dmtdt, derivt
263      double precision afact2
264      type(nwad_dble)::rhoval
265
266c     functional derivatives above FFFFFFFFFFFF
267
268      parameter(kap=0.8040d0, mu=0.21951d0)
269      parameter (F43=4.d0/3.d0, F13=1.d0/3.d0)
270      parameter (F83=8.d0/3.d0, F23=2.d0/3.d0, F53=5.d0/3.d0)
271      parameter (G920  =9.d0/20.d0 )
272
273      parameter(b=0.40d0, c=1.59096d0, e=1.537d0)
274      parameter (C1  =  10.d0/81.d0,
275     &     C2  = 146.d0/2025.d0,
276     &     C3  = -73.d0/405.d0 )
277c
278      pi=acos(-1d0)
279      Ax = (-0.75d0)*(3d0/pi)**F13
280      P32 = (3.d0*pi**2)**F23
281      es=dsqrt(e)
282      afact2=1d0/facttwo
283c
284      do n = 1, nq
285         rhoval=rho(n)*facttwo
286         if (rhoval.ge.tol_rho) then
287
288c     rho43= n*e_x^unif=exchange energy per electron for uniform electron gas
289c     = n* Ax*n^(1/3)   or n*C*n^(1/3)
290
291            rho43 = Ax*rhoval**F43 ! Ax*n^4/3
292            rrho = 1d0/rhoval   ! reciprocal of rho
293c           rho13 = rho43*rrho
294
295C     Below we just sum up the LDA contribution to the functional
296            func(n)= func(n) + rho43*wght*afact2
297c           Amat(n) = Amat(n) + F43*rho13*wght
298
299c
300c           gamma = delrho(n,1)*delrho(n,1) +
301c    &           delrho(n,2)*delrho(n,2) +
302c    &           delrho(n,3)*delrho(n,3)
303            gamma=rgamma(n)
304            gamma=gamma*facttwo*facttwo
305            tauN = tau(n)*facttwo
306            tauW=0.125d0*gamma*rrho
307            tauU=0.3d0*P32*rhoval**F53
308c
309c     Evaluate the Fx, i.e. mt(x) = Fx - 1 (LDA bit already there)
310c
311              p=gamma/(rhoval**F83*P32*4.d0)
312              if (tauN.ge.tol_rho) then
313                z=tauW/tauN
314              else
315                z=0.0d0
316              endif
317              al=(tauN-tauW)/tauU
318c     al=dabs(al)
319              if(al.lt.0d0)  al=0d0
320
321              qtil=(G920*(al-1.d0)/((1.d0+b*al*(al-1.d0))**.5d0)) +
322     +             F23*p
323
324              xb=(c*z**2)/( (1.0d0+z**2)**2)
325              x1=(C1 + xb)*p
326              x2=C2*qtil*qtil
327              xc=C3*qtil
328              if (gamma.lt.tol_rho) then
329                xd = 0.0d0
330              else
331                xd=(0.5d0*(.6d0*z)**2 + .5d0*p**2)**.5d0
332              endif
333              x3=xc*xd
334              x4=C1*C1*p**2/kap
335              x5=2.d0*es*C1*(.6d0*z)**2
336              x6= e*mu*p**3
337              x7 = (1.d0+es*p)**(-2)
338
339              x=(x1+x2+x3 +x4 +x5+x6)*x7
340
341              if (abs(x).lt.tol_rho) write(0,*) ' x for fx ',x
342
343c     functional derivatives FFFFFFFFFFFFFFFFFFFFFFFFFFFF
344
345C     Derivatives wrt n, density   below
346c           dzdn=-z*rrho
347c           dpdn = -p*rrho*F83
348c           daldn=F53*( -p*dzdn/z**2 +dpdn*(-1.d0+1.d0/z) )
349c
350c           til1=al-1.d0
351c           til2=(1.d0+b*al*(al-1.d0))**(-0.5d0)
352c           dtil1dn=daldn
353c           dtil2dn=b*daldn*(2.d0*al-1d0)*
354c    &           (-.5d0)*(til2**3)
355c           dqtdn = G920*(til2*dtil1dn+til1*dtil2dn)+F23*dpdn
356c
357c           ax1=c*p*z*z
358c           bx1=(1+z*z)**(-2.d0)
359c           dx1dn=(x1/p)*dpdn + 2d0*c*p*z*dzdn/((1d0+z*z)**3)*(1d0-z*z)
360c           dx2dn=2.d0*C2*qtil*dqtdn
361c
362c           dxddn=.5d0/xd*( (3d0/5d0)**2*z*dzdn +
363c    +           p*dpdn)
364c           dxcdn=C3*dqtdn
365c           dx3dn=xc*dxddn+xd*dxcdn
366c
367c           dx4dn=(2.d0*x4/p)*dpdn
368c           dx5dn=(2.d0*x5/z)*dzdn
369c           dx6dn=(3.d0*x6/p)*dpdn
370c           dx7dn=-2.d0*es*dpdn/(1.d0+es*p)**3
371c
372c           xmany=x1+x2+x3 +x4 +x5+x6
373c           dxmanydn= dx1dn+dx2dn+dx3dn+dx4dn+dx5dn+dx6dn
374c           dxdn=x7*dxmanydn+xmany*dx7dn
375C     Derivatives wrt n, density   above
376
377C     Derivatives wrt gamma,    below
378
379c           dpdg=p/gamma
380c           dzdg=z/gamma
381c           daldg=(al/p)*dpdg-F53*(p/(z*z))*dzdg
382c
383c           dtil2dg=-0.5d0*daldg*b*(2.d0*al-1d0)*til2**3.d0
384c           dqtdg=G920*(til1*dtil2dg+til2*daldg)+F23*dpdg
385c           dx1dg=(x1/p)*dpdg + 2d0*c*p*z*dzdg/((1d0+z*z)**3)*(1d0-z*z)
386c
387c           dx2dg=C2*2.d0*qtil*dqtdg
388c
389c           dxcdg=C3*dqtdg
390c           dxddg=.5d0/xd*( (3d0/5d0)**2*z*dzdg +
391c    +           p*dpdg)
392c           dx3dg=xc*dxddg+xd*dxcdg
393c
394c           dx4dg=(2.d0*x4/p)*dpdg
395c           dx5dg=(2.d0*x5/z)*dzdg
396c           dx6dg=(3.d0*x6/p)*dpdg
397c
398c           dx7dg=-2.d0*es*dpdg*(1.d0+p*es)**(-3.d0)
399c
400c           dxmanydg= dx1dg+dx2dg+dx3dg+dx4dg+dx5dg+dx6dg
401c           dxdg=x7*dxmanydg+xmany*dx7dg
402
403C     Derivatives wrt tau,    below
404c     ttttttttttttttttttttttttttttttttttttttttttttttttt
405c           dzdt= -z/tauN
406c           daldt=1.d0/tauU
407c
408c           dqtdt=g920*daldt*til2*(1d0-
409c    -           0.5d0*b*til1*til2*til2*(2d0*al-1d0))
410c
411c           dx1dt=c*p*dzdt*2d0*z*(1d0-z*z)/((1.d0+z*z)**3)
412c           dx2dt=2*c2*qtil*dqtdt
413c           dx3dt=x3*(dqtdt/qtil +
414c    &           0.5d0*(3d0/5d0)**2*z*dzdt/(xd*xd))
415c           dx5dt=2d0*(x5/z)*dzdt
416c
417c           dxmanydt= dx1dt+dx2dt+dx3dt+dx5dt
418c           dxdt=x7*dxmanydt
419c     ttttttttttttttttttttttttttttttttttttttttttttttttttt
420
421              mt = kap - kap/(1.d0 + x/kap)
422
423              func(n)= func(n) + mt*rho43*wght*afact2
424
425c           dmtdn=dxdn/(1.d0+x/kap)**2
426c           derivn=mt*F43*rho13+rho43*dmtdn
427
428c           dmtdg=dxdg/(1.d0+x/kap)**2
429c           derivg = rho43*dmtdg
430c
431c           dmtdt=dxdt/(1.d0+x/kap)**2
432c           derivt = rho43*dmtdt
433c           Amat(n) = Amat(n) + derivn*wght
434c
435c     4x factor comes from gamma_aa = gamma_total/4
436c
437c           Cmat(n)=  Cmat(n) + 2d0*derivg*wght
438c           Mmat(n)=  Mmat(n) +0.5d0*derivt*wght
439         endif
440      enddo
441      return
442      end
443#ifndef NWAD_PRINT
444#define NWAD_PRINT
445c
446c     Compile source again for the 2nd derivative case
447c
448#include "nwxc_x_tpss03.F"
449#endif
450#ifndef SECOND_DERIV
451#define SECOND_DERIV
452c
453c     Compile source again for the 2nd derivative case
454c
455#include "nwxc_x_tpss03.F"
456#endif
457#ifndef THIRD_DERIV
458#define THIRD_DERIV
459c
460c     Compile source again for the 3rd derivative case
461c
462#include "nwxc_x_tpss03.F"
463#endif
464#undef NWAD_PRINT
465C>
466C> @}
467