1      subroutine tddft_grad_gxc(ipol,nroot,nbf,ibf,nq,chi,dchi,xpy_rho,
2     &                          amat3,scr,scrmat)
3      implicit none
4c
5c     $Id$
6c
7c     Compute the DFT contributions to R and W in AO basis. This routine
8c     will be called from within the DFT drivers, so it has to be
9c     adapted to the idiosyncracies thereof.
10c
11c     The general expression to evaluate is
12c
13c        Rijs = Rijs + Sum_klt Sum_mnu g_ijskltmnu*(X+Y)_klt*(X+Y)_mnu
14c
15c     where i,j,k,l,m and n are AO labels and s,t and u are spin labels.
16c     The quantity g is the third derivative of the functional with
17c     respect to the density matrix, i.e.:
18c
19c                      d3 f(ra,rb,gaa,gab,gbb,ta,tb)
20c        g_ijskltmnu = -----------------------------
21c                      d D_ijs * d D_klt * d D_mnu
22c
23c     Grouping factors in the expression at the top it can be reduced
24c     to an equation involving 3 different kinds of factors:
25c
26c       1. partial derivatives of the functional
27c       2. "basis function" factors
28c       3. "density" (like) factors
29c
30c     The first kind of factors are obvious, the basis function factors
31c     are:
32c
33c       X_i(r)*X_j(r)                                          (101)
34c       ... more for (meta) GGAs ...
35c
36c     The density factors are:
37c
38c       Sum_ij X_i(r)*X_j(r)*(X+Y)_ijs                         (201)
39c       ... more for (meta) GGAs ...
40c
41c     In particular computing the density factors up front helps
42c     eliminating indeces early and is therefore recommended to
43c     enhance the efficiency. As these expressions are identical to
44c     the normal density expressions the same machinery can be used
45c     to evaluate them. Note that typical density evaluation routines
46c     will assume that the density presented to them is symmetric.
47c     Due to the form of the expressions this can be dealt with by
48c     symmetrizing (X+Y) before passing it into the density evaluators.
49c
50c     Written by Huub van Dam, November 2006.
51c
52#include "dft3drv.fh"
53#include "global.fh"
54c
55c     Input:
56c
57      integer ipol     ! =1 (restricted), =2 (unrestricted)
58      integer nroot    ! the number of roots
59      integer nbf      ! the number of active basis functions
60      integer ibf(nbf) ! the table giving the true basis function number
61                       ! for every element of compressed list
62      integer nq       ! the number of grid points
63c
64      double precision chi(nq,nbf)            ! basis function values
65      double precision dchi(nq,3,nbf)         ! basis function gradients
66      double precision xpy_rho(nq,ipol,nroot) ! density factor eq. (201)
67c
68c     3rd order functional derivatives wrt rho
69      double precision amat3(nq,NCOL_AMAT3)
70c
71c     Input/Output:
72c
73c     buffer to accumulate R matrix contributions into:
74c
75      double precision scrmat(ipol*nroot*nbf*nbf)
76c
77c     integer g_r ! the global array holding the TDDFT gradient
78c                 ! right-hand-side. The contribitions are added onto
79c                 ! this.
80c
81c     Workspace:
82c
83      double precision scr(nq) ! buffer to minimize multiplications
84c
85c     Local:
86c
87      integer kbf, lbf ! counters over basis functions
88      integer iq       ! counter over quadrature points
89      integer ir       ! counter over roots
90c
91      double precision t ! a term
92c
93c     Code:
94c
95cDEBUG
96      logical oroot
97      oroot = ga_nodeid().eq.0
98c     if (oroot) write(*,*)'*** tddft_grad_gxc begin'
99c     if (oroot) write(*,*)'*** basis functions'
100c     do kbf=1,nbf
101c       do iq=1,nq
102c         if (oroot) write(*,*)'*** chi=',kbf,iq,chi(iq,kbf)
103c       enddo
104c     enddo
105c     if (oroot) write(*,*)'*** amat3'
106c     do iq=1,nq
107c       if (oroot) write(*,*)'*** amat3=',amat3(iq,D3_RA_RA_RA)
108c     enddo
109c     if (oroot) write(*,*)'*** xpy_rho'
110c     do ir=1,nroot
111c       do iq=1,nq
112c         if (oroot) write(*,*)'*** xpy_rho=',ir,iq,xpy_rho(iq,1,ir)
113c       enddo
114c     enddo
115c     if (oroot) write(*,*)'*** tddft_grad_gxc end'
116cDEBUG
117      if (ipol.eq.1) then
118c
119c       Only one spin component to do
120c
121        do ir = 1, nroot
122          do kbf = 1, nbf
123            do iq = 1, nq
124              scr(iq) =
125     &            chi(iq,kbf)*amat3(iq,D3_RA_RA_RA)* ! gaaa
126     &            xpy_rho(iq,1,ir)*xpy_rho(iq,1,ir)
127     &        + 2*chi(iq,kbf)*amat3(iq,D3_RA_RA_RB)* ! gaab+gaba
128     &            xpy_rho(iq,1,ir)*xpy_rho(iq,1,ir)
129     &        +   chi(iq,kbf)*amat3(iq,D3_RA_RA_RB)* ! gabb
130     &            xpy_rho(iq,1,ir)*xpy_rho(iq,1,ir)
131              scr(iq) = 0.25d0*scr(iq) ! xpy_rho: total rho -> alpha rho
132              scr(iq) = 2.0d0*scr(iq)  ! 2*gxc
133cDEBUG
134c     scr(iq) = 0.d0
135cDEBUG
136            enddo
137            do lbf = 1, kbf-1
138              t = 0.0d0
139              do iq = 1, nq
140                t = t + chi(iq,lbf)*scr(iq)
141              enddo
142              scrmat(ir+(kbf-1)*nroot+(lbf-1)*nroot*nbf) = t
143              scrmat(ir+(lbf-1)*nroot+(kbf-1)*nroot*nbf) = t
144            enddo
145            t = 0.0d0
146            do iq = 1, nq
147              t = t + chi(iq,kbf)*scr(iq)
148            enddo
149            scrmat(ir+(kbf-1)*nroot+(kbf-1)*nroot*nbf) = t
150          enddo
151        enddo
152c
153      else if (ipol.eq.2) then
154c
155c       Alpa density component first
156c
157        do ir = 1, nroot
158          do kbf = 1, nbf
159            do iq = 1, nq
160              scr(iq) = chi(iq,kbf)*(
161     &                    amat3(iq,D3_RA_RA_RA)*
162     &                    xpy_rho(iq,1,ir)*xpy_rho(iq,1,ir)
163     &                 +2*amat3(iq,D3_RA_RA_RB)*
164     &                    xpy_rho(iq,1,ir)*xpy_rho(iq,2,ir)
165     &                 +  amat3(iq,D3_RA_RB_RB)*
166     &                    xpy_rho(iq,2,ir)*xpy_rho(iq,2,ir)
167     &                  )
168            enddo
169            do lbf = 1, kbf-1
170              t = 0.0d0
171              do iq = 1, nq
172                t = t + chi(iq,lbf)*scr(iq)
173              enddo
174              scrmat(ir+(lbf-1)*nroot*ipol+(kbf-1)*nroot*ipol*nbf) = t
175              scrmat(ir+(kbf-1)*nroot*ipol+(lbf-1)*nroot*ipol*nbf) = t
176            enddo
177            t = 0.0d0
178            do iq = 1, nq
179              t = t + chi(iq,kbf)*scr(iq)
180            enddo
181            scrmat(ir+(kbf-1)*nroot*ipol+(kbf-1)*nroot*ipol*nbf) = t
182          enddo
183        enddo
184c
185c       Beta density component next
186c
187        do ir = 1, nroot
188          do kbf = 1, nbf
189            do iq = 1, nq
190              scr(iq) = chi(iq,kbf)*(
191     &                    amat3(iq,D3_RA_RA_RB)*
192     &                    xpy_rho(iq,1,ir)*xpy_rho(iq,1,ir)
193     &                 +2*amat3(iq,D3_RA_RB_RB)*
194     &                    xpy_rho(iq,1,ir)*xpy_rho(iq,2,ir)
195     &                 +  amat3(iq,D3_RB_RB_RB)*
196     &                    xpy_rho(iq,2,ir)*xpy_rho(iq,2,ir)
197     &                  )
198            enddo
199            do lbf = 1, kbf-1
200              t = 0.0d0
201              do iq = 1, nq
202                t = t + chi(iq,lbf)*scr(iq)
203              enddo
204              scrmat(ir+(ipol-1)*nroot+(lbf-1)*nroot*ipol
205     &               +(kbf-1)*nroot*ipol*nbf) = t
206              scrmat(ir+(ipol-1)*nroot+(kbf-1)*nroot*ipol
207     &               +(lbf-1)*nroot*ipol*nbf) = t
208            enddo
209            t = 0.0d0
210            do iq = 1, nq
211              t = t + chi(iq,kbf)*scr(iq)
212            enddo
213            scrmat(ir+(ipol-1)*nroot+(kbf-1)*nroot*ipol
214     &             +(kbf-1)*nroot*ipol*nbf) = t
215          enddo
216        enddo
217      endif
218c
219      end
220