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