1cc AJL/Begin/FDE
2c
3c     Wrapper routine for evaluating Ts functional and derivatives,
4c     and combining them with quadrature weights.
5c     Copied from xc_eval_fnl, with minor adaptation
6c
7c     AJL (10/2014)
8c
9C> \ingroup nwdft_xc
10C> @{
11C>
12C> \brief the functional evaluation routine
13C>
14C> This routine iterates over the components of the current density
15C> functional and evaluates and accumulates all density dependent
16C> terms.
17C> At present it will evaluate
18C>
19C> - the kinetic energy energy
20C>
21C> - the kinetic energy 1st order partial derivatives wrt density
22C>
23C> as desired.
24C>
25      Subroutine ts_eval_fnl(rho, delrho, Amat, Amat2, Cmat, Cmat2,
26     &                       nq, Ts, qwght, grad, ldew, func,
27     &                       do_2nd, ttau, kske, Mmat, Mmat2,
28     &                       do_3rd, Amat3, Cmat3, ldew2)
29c
30      implicit none
31c
32#include "cdft.fh"
33c xc-second derivative header file
34#include "dft2drv.fh"
35c xc-third derivative header file
36#include "dft3drv.fh"
37#include "stdio.fh"
38c#include "steric.fh"
39c
40cc AJL/Unused
41c      integer nq, is12x
42      integer nq
43      double precision rho(*)    !< [Input] The electron density
44                                 !< \f$\rho^\alpha\f$ and
45                                 !< \f$\rho^\beta\f$.
46      double precision delrho(*) !< [Input] The electron density gradient
47                                 !< \f$\partial\rho^\alpha/\partial r\f$
48                                 !< and
49                                 !< \f$\partial\rho^\beta/\partial r\f$
50c
51      double precision Amat(*), Cmat(*)
52      double precision Amat2(*),Cmat2(*)
53      double precision Amat3(*),Cmat3(*)
54      double precision ttau(*) !< [Input] The kinetic energy density
55      double precision Mmat(*), Mmat2(*)
56c
57      double precision Ts !< [Output] The exchange energy
58      double precision qwght(nq), func(nq)
59cc AJL/Unused
60c      logical grad, kske, dohcth
61      logical grad, kske
62c
63c      character*4 whichf
64
65      logical do_2nd, do_3rd
66      logical ldew, ldew2, ldew3
67c
68      double precision eps
69c AJL/Unused
70c      integer nx,nc,dumi
71      parameter (eps=1.e-8)
72c
73c     Initialize the Ts potential and energy sampling matrices.
74c
75      call dfill(ipol*nq, 0.d0, Amat, 1)
76      if(grad) call dfill(3*nq*ipol, 0.d0, Cmat, 1)
77      if(kske) call dfill(nq*ipol, 0.d0, Mmat, 1)
78      if (do_2nd) then
79         call dfill(nq*NCOL_AMAT2, 0.d0, Amat2, 1)
80         if (grad) call dfill(nq*NCOL_CMAT2, 0.d0, Cmat2, 1)
81         if (kske) call dfill(nq*NCOL_MMAT2, 0.d0, Mmat2, 1)
82      endif
83c
84c     Initialize the 3rd derivatives.  Note that we may need
85c     the 2nd derivative matrices in some cases as well, so
86c     those are initialized also.
87      if (do_3rd) then
88         call dfill(nq*NCOL_AMAT2, 0.0d0, Amat2, 1)
89         call dfill(nq*NCOL_AMAT3, 0.0d0, Amat3, 1)
90         if (grad) call dfill(nq*NCOL_CMAT2, 0.0d0, Cmat2, 1)
91         if (grad) call dfill(nq*NCOL_CMAT3, 0.0d0, Cmat3, 1)
92      endif
93c
94c     This prevents the Ts-kernel from being multiplied
95c     with the quadrature weights in TDDFT gradients.
96c
97      ldew3 = (do_3rd.and.ldew)
98      if (ldew.or.ldew2.or.ldew3) call dfill(nq, 0.d0, func, 1)
99c      if (ldew.or.ldew3) call dfill(nq, 0.d0, func, 1)
100c      if (ldew) call dfill(nq, 0.d0, func, 1)
101c
102      if (tsfac(1).gt.eps) then
103         if (.not. do_2nd .and. .not. do_3rd) then
104            call ts_tf(tol_rho, tsfac(1), 0.0d0, 0.0d0, rho,
105     &           Amat, nq, ipol, Ts, qwght,
106     &           ldew, func)
107c         else if (.not. do_3rd) then
108c            call ts_tf_d2(tol_rho, 1.d0, 0.d0, 0.d0, rho,
109c     &           Amat, Amat2, nq, ipol, Ts, qwght,
110c     &           .false., func)
111c         else
112c            call ts_tf_d3(tol_rho, 1.d0, 0.d0, 0.d0, rho,
113c     &           Amat, Amat2, Amat3, nq, ipol, Ts, qwght,
114c     &           .false., func)
115         endif
116      endif
117
118      if(tsfac(2).gt.eps) then
119         if (.not. do_2nd .and. .not. do_3rd) then
120            call ts_vw(tol_rho, tsfac(2), 0.0d0, 0.0d0, rho,
121     &           delrho, Amat, Cmat, nq, ipol,
122     &           Ts, qwght, ldew, func)
123!         else if (.not. do_3rd) then
124!            call ts_vw_d2(tol_rho, 1.d0, 0.d0, 0.d0, rho,
125!     &           Amat, Amat2, nq, ipol, Ts, qwght,
126!     &           .false., func)
127!         else
128!            call ts_vw_d3(tol_rho, 1.d0, 0.d0, 0.d0, rho,
129!     &           Amat, Amat2, Amat3, nq, ipol, Ts, qwght,
130!     &           .false., func)
131         endif
132      endif
133c
134c     Calculate the steric energy
135c
136c      if (lsteric) then
137c        StericEnergy = 0.d0
138c        call steric_energy(tol_rho,xfac(1),rho,delrho,nq,
139c     &   qwght,ipol,StericEnergy)
140c      endif
141c
142c     Combine with quadrature weights
143c
144      if (.not.(ldew2.or.ldew3)) then
145
146        if (.not. do_2nd .and. .not. do_3rd) then
147           call setACmat(delrho, Amat, Cmat, qwght, ipol, nq, grad,
148     &           (.not. do_2nd), kske, Mmat, .false., 0d0)
149c        else if (.not. do_3rd) then
150c           call setACmat_d2(delrho, Amat, Amat2, Cmat, Cmat2, qwght,
151c     &           ipol, nq, grad, (.not. do_2nd), kske, Mmat, Mmat2,
152c     &           .false.)
153c        else
154c           call setACmat_d3(delrho, Amat, Amat2, Amat3, Cmat, Cmat2,
155c     &           Cmat3, qwght, ipol, nq, grad, (.not. do_3rd), .false.,
156c     &.          false.)
157        endif
158      endif
159c
160      return
161      end
162
163cC>
164cC> @}
165cc AJL/End
166