1c 2c Modified to handle second derivatives while reusing code 3c 4c BGJ - 8/98 5c 6C> \ingroup nwdft_xc 7C> @{ 8C> 9C> \brief Do clever stuff to the partial derivatives 10C> 11C> The density functionals return plain partial derivatives of the 12C> of the functionals in addition to the energy. This routine modifies 13C> those quantities to get them in a form that is more "convenient" 14C> in the construction of various matrices. 15C> 16C> In all cases this modification includes multiplying the partial 17C> derivatives for each grid point with the corresponding quadrature 18C> weight. 19C> 20C> If `xform_Cmat` is set to `.TRUE.` then the 1st derivatives with 21C> respect to the norm of the electron density gradient are replaced 22C> in the closed shell case with 23C> \f{eqnarray*}{ 24C> Cmat(i,1,1) &=& \frac{\partial\rho(r)}{\partial x} \left( 25C> \frac{\partial f}{\partial\gamma^{\alpha\alpha}} + 26C> \frac{1}{2}\frac{\partial f}{\partial\gamma^{\alpha\beta}} 27C> \right) w_i\\\\ 28C> Cmat(i,2,1) &=& \frac{\partial\rho(r)}{\partial y} \left( 29C> \frac{\partial f}{\partial\gamma^{\alpha\alpha}} + 30C> \frac{1}{2}\frac{\partial f}{\partial\gamma^{\alpha\beta}} 31C> \right) w_i\\\\ 32C> Cmat(i,3,1) &=& \frac{\partial\rho(r)}{\partial z} \left( 33C> \frac{\partial f}{\partial\gamma^{\alpha\alpha}} + 34C> \frac{1}{2}\frac{\partial f}{\partial\gamma^{\alpha\beta}} 35C> \right) w_i\\\\ 36C> \f} 37C> and in the open shell case with 38C> \f{eqnarray*}{ 39C> Cmat(i,1,1) &=& \left(2\frac{\partial\rho^\alpha(r)}{\partial x} 40C> \frac{\partial f}{\partial\gamma^{\alpha\alpha}} + 41C> \frac{\partial\rho^\beta(r)}{\partial x} 42C> \frac{\partial f}{\partial\gamma^{\alpha\beta}}\right)w_i \\\\ 43C> Cmat(i,2,1) &=& \left(2\frac{\partial\rho^\alpha(r)}{\partial y} 44C> \frac{\partial f}{\partial\gamma^{\alpha\alpha}} + 45C> \frac{\partial\rho^\beta(r)}{\partial y} 46C> \frac{\partial f}{\partial\gamma^{\alpha\beta}}\right)w_i \\\\ 47C> Cmat(i,3,1) &=& \left(2\frac{\partial\rho^\alpha(r)}{\partial z} 48C> \frac{\partial f}{\partial\gamma^{\alpha\alpha}} + 49C> \frac{\partial\rho^\beta(r)}{\partial z} 50C> \frac{\partial f}{\partial\gamma^{\alpha\beta}}\right)w_i \\\\ 51C> Cmat(i,1,2) &=& \left(2\frac{\partial\rho^\beta(r)}{\partial x} 52C> \frac{\partial f}{\partial\gamma^{\beta\beta}} + 53C> \frac{\partial\rho^\alpha(r)}{\partial x} 54C> \frac{\partial f}{\partial\gamma^{\alpha\beta}}\right)w_i \\\\ 55C> Cmat(i,2,2) &=& \left(2\frac{\partial\rho^\beta(r)}{\partial y} 56C> \frac{\partial f}{\partial\gamma^{\beta\beta}} + 57C> \frac{\partial\rho^\alpha(r)}{\partial y} 58C> \frac{\partial f}{\partial\gamma^{\alpha\beta}}\right)w_i \\\\ 59C> Cmat(i,3,2) &=& \left(2\frac{\partial\rho^\beta(r)}{\partial z} 60C> \frac{\partial f}{\partial\gamma^{\beta\beta}} + 61C> \frac{\partial\rho^\alpha(r)}{\partial z} 62C> \frac{\partial f}{\partial\gamma^{\alpha\beta}}\right)w_i \\\\ 63C> \f} 64C> The grid point is indexed with `i`, the density functional is 65C> referred to with \f$f\f$, and \f$w\f$ is the quadrature weight. 66C> 67#if !defined SECOND_DERIV && !defined THIRD_DERIV 68 Subroutine setACmat(delrho, Amat, Cmat, qwght, ipol, nq, GRAD, 69 & xform_Cmat, kske, Mmat, kslap, Lmat) 70#elif defined(SECOND_DERIV) && !defined THIRD_DERIV 71c Second derivatives are not yet implemented for the meta-GGAs 72 Subroutine setACmat_d2(delrho, Amat, Amat2, Cmat, Cmat2, qwght, 73 & ipol, nq, GRAD, xform_Cmat, kske, Mmat, Mmat2, 74 & kslap) 75#else 76c Third derivatives are not yet implemented for the meta-GGAs 77 Subroutine setACmat_d3(delrho, Amat, Amat2, Amat3, Cmat, Cmat2, 78 & Cmat3, qwght, ipol, nq, GRAD, xform_Cmat, kske, kslap) 79#endif 80c 81C$Id$ 82c 83 implicit none 84c 85#include "dft2drv.fh" 86#include "dft3drv.fh" 87#include "stdio.fh" 88c !!! BGJ test 89#include "bgj.fh" 90c !!! BGJ test 91c 92 integer ipol !< [Input] The number of spin channels 93 integer nq !< [Input] The number of grid points 94 Logical GRAD !< [Input] If `.TRUE.` the density functional is 95 !< a GGA 96 logical xform_Cmat !< [Input] If `.TRUE.` do the elaborate 97 !< transformation 98c 99c Density gradients - used for transforming fnl gamma derivatives 100c 101 double precision delrho(nq,3,ipol) !< [Input] The gradient of the 102 !< density 103c 104c Sampling Matrices for the XC Potential & Energy 105c 106 double precision Amat(nq,ipol) !< [In/Output] \f$\partial f/\partial\rho^\sigma\f$ 107 double precision Cmat(nq,3,ipol) !< [In/Output] \f$\partial f/\partial\gamma^{\sigma\sigma'}\f$ 108 double precision Mmat(nq,ipol) !< [In/Output] \f$\partial f/\partial\tau^\sigma\f$ 109 logical kske !< [Input] If `.TRUE.` the density functional is 110 !< a meta-GGA 111c 112 logical kslap !< [Input] If `.TRUE.` the density functional is 113 !< a laplacian-dependent meta-GGA 114 double precision Lmat(nq,ipol) !< [In/Output] \f$\partial f/\partial\lap^\sigma\f$ 115 116c 117#ifdef SECOND_DERIV 118 double precision Amat2(nq,NCOL_AMAT2) !< [In/Output] Similar to 119 !< `Amat` but 2nd order derivatives 120 double precision Cmat2(nq,NCOL_CMAT2) !< [In/Output] Similar to 121 !< `Cmat` but 2nd order derivatives 122 double precision Mmat2(nq,NCOL_MMAT2) !< [In/Output] Similar to 123 !< `Mmat` but 2nd order derivatives 124#endif 125c 126#ifdef THIRD_DERIV 127 double precision Amat3(nq,NCOL_AMAT3) 128 double precision Cmat3(nq,NCOL_CMAT3) 129#endif 130c 131c Quadrature Weights 132c 133 double precision qwght(nq) !< [Input] The quadrature weights 134c 135 integer ii, jj 136c 137c 138c Transform derivatives of functional with respect to gammas 139c to derivatives of functional with respect to density gradient 140c 141 if (GRAD .and. xform_Cmat) then 142 call transform_Cmat(delrho, Cmat, ipol, nq) 143 endif 144c 145c Combine derivatives of functional with quadrature weights 146c 147 if (GRAD)then 148 if (xform_Cmat) then 149 do ii = 1, ipol 150 do jj = 1, nq 151 Amat(jj,ii) = Amat(jj,ii)*qwght(jj) 152 Cmat(jj,1,ii) = Cmat(jj,1,ii)*qwght(jj) 153 Cmat(jj,2,ii) = Cmat(jj,2,ii)*qwght(jj) 154 Cmat(jj,3,ii) = Cmat(jj,3,ii)*qwght(jj) 155 enddo 156 enddo 157 else 158 if (ipol .eq. 1) then 159 do jj = 1, nq 160 Amat(jj,1) = Amat(jj,1)*qwght(jj) 161 Cmat(jj,D1_GAA,1) = Cmat(jj,D1_GAA,1)*qwght(jj) 162 Cmat(jj,D1_GAB,1) = Cmat(jj,D1_GAB,1)*qwght(jj) 163 enddo 164 else 165 do jj = 1, nq 166 Amat(jj,1) = Amat(jj,1)*qwght(jj) 167 Amat(jj,2) = Amat(jj,2)*qwght(jj) 168 Cmat(jj,D1_GAA,1) = Cmat(jj,D1_GAA,1)*qwght(jj) 169 Cmat(jj,D1_GAB,1) = Cmat(jj,D1_GAB,1)*qwght(jj) 170 Cmat(jj,D1_GBB,1) = Cmat(jj,D1_GBB,1)*qwght(jj) 171 enddo 172 endif 173 endif 174 else 175 do ii = 1, ipol 176 do jj = 1, nq 177 Amat(jj,ii) = Amat(jj,ii)*qwght(jj) 178 enddo 179 enddo 180 endif 181 if (kske) then 182 do ii = 1, ipol 183 do jj = 1, nq 184 Mmat(jj,ii) = Mmat(jj,ii)*qwght(jj) 185 enddo 186 enddo 187 endif 188 if (kslap) then 189 do ii = 1, ipol 190 do jj = 1, nq 191 Lmat(jj,ii) = Lmat(jj,ii)*qwght(jj) 192 enddo 193 enddo 194 endif 195c 196#if 0 197 if (bgj_print() .gt. 1) then 198 write(LuOut,*) ' setACmat: AMAT out' 199 call output(amat, 1, nq, 1, ipol, nq, ipol, 1) 200 if (grad) then 201 write(LuOut,*) ' setACmat: CMAT out ',xform_Cmat 202 call output(cmat, 1, nq, 1, 3*ipol, nq, 3*ipol, 1) 203 endif 204 if (kske) then 205 write(LuOut,*) ' setACmat: MMAT out' 206 call output(mmat, 1, nq, 1, ipol, nq, ipol, 1) 207 endif 208 endif 209#endif 210c 211#if defined(SECOND_DERIV) || defined(THIRD_DERIV) 212 do ii = 1, NCOL_AMAT2 213 do jj = 1, nq 214 Amat2(jj,ii) = Amat2(jj,ii)*qwght(jj) 215 enddo 216 enddo 217 if (GRAD)then 218 do ii = 1, NCOL_CMAT2 219 do jj = 1, nq 220 Cmat2(jj,ii) = Cmat2(jj,ii)*qwght(jj) 221 enddo 222 enddo 223 endif 224 if (kske)then 225 do ii = 1, NCOL_MMAT2 226 do jj = 1, nq 227 Mmat2(jj,ii) = Mmat2(jj,ii)*qwght(jj) 228 enddo 229 enddo 230 endif 231#if 0 232 if (bgj_print() .gt. 1) then 233 write(LuOut,*) ' setACmat_d2: AMAT2 out' 234 call output(amat2, 1, nq, 1, NCOL_AMAT2, nq, NCOL_AMAT2, 1) 235 if (grad) then 236 write(LuOut,*) ' setACmat_d2: CMAT2 out' 237 call output(cmat2, 1, nq, 1, NCOL_CMAT2, nq, NCOL_CMAT2, 1) 238 endif 239 endif 240#endif 241#endif 242c 243#ifdef THIRD_DERIV 244 do ii = 1, NCOL_AMAT3 245 do jj = 1, nq 246 Amat3(jj,ii) = Amat3(jj,ii)*qwght(jj) 247 enddo 248 enddo 249 if (GRAD) then 250 do ii = 1, NCOL_CMAT3 251 do jj = 1, nq 252 Cmat3(jj,ii) = Cmat3(jj,ii)*qwght(jj) 253 enddo 254 enddo 255 endif 256#endif 257 return 258 end 259c 260#if !defined(SECOND_DERIV) && !defined(THIRD_DERIV) 261c 262c Transform Cmat from gamma form to delrho form 263c 264c BGJ - 8/98 265c 266 Subroutine transform_Cmat(delrho, Cmat, ipol, nq) 267c 268 implicit none 269c 270#include "dft2drv.fh" 271c 272 integer ipol, nq 273 double precision delrho(nq,3,ipol), Cmat(nq,3,ipol) 274c 275 integer n 276 double precision gaa, gab, gbb 277c 278 if (ipol .eq. 1) then 279 do n = 1, nq 280c Must account for delrho being total density gradient, not alpha 281 gaa = Cmat(n,D1_GAA,1) + Cmat(n,D1_GAB,1)*0.5d0 282 Cmat(n,1,1) = delrho(n,1,1)*gaa 283 Cmat(n,2,1) = delrho(n,2,1)*gaa 284 Cmat(n,3,1) = delrho(n,3,1)*gaa 285 enddo 286 else 287 do n = 1, nq 288 gaa = Cmat(n,D1_GAA,1) 289 gab = Cmat(n,D1_GAB,1) 290 gbb = Cmat(n,D1_GBB,1) 291 Cmat(n,1,1) = 2d0*delrho(n,1,1)*gaa + delrho(n,1,2)*gab 292 Cmat(n,2,1) = 2d0*delrho(n,2,1)*gaa + delrho(n,2,2)*gab 293 Cmat(n,3,1) = 2d0*delrho(n,3,1)*gaa + delrho(n,3,2)*gab 294 Cmat(n,1,2) = 2d0*delrho(n,1,2)*gbb + delrho(n,1,1)*gab 295 Cmat(n,2,2) = 2d0*delrho(n,2,2)*gbb + delrho(n,2,1)*gab 296 Cmat(n,3,2) = 2d0*delrho(n,3,2)*gbb + delrho(n,3,1)*gab 297 enddo 298 endif 299c 300 return 301 end 302#endif 303c 304#ifndef SECOND_DERIV 305#define SECOND_DERIV 306c 307c Compile source again for the 2nd derivative case 308c 309#include "setACmat.F" 310#endif 311#ifndef THIRD_DERIV 312#define THIRD_DERIV 313c 314c Compile source again for the 3rd derivative case 315c 316#include "setACmat.F" 317#endif 318C> @} 319