1 subroutine xc_3rd_deriv(prho, pdelrho, 2 1 Amat3, Cmat3, Cmat2, delrho, npert, ipol, nq, grad, triplet) 3c $Id$ 4c 5c Calculates the 3rd order derivatives of the XC-functional, which are 6c required for evaluating quadratic response functions with TDDFT. 7c This is not yet implemented for meta-GGA functionals. Explicitly, we 8c begin construction of: 9c 10c G^{xc} = g^{xc}*(X+Y)*(X+Y) 11c 12c This is quantity is later assembled in the AO basis by xc_tabcd. 13c 14c Called from: grid_quadv0b 15c 16 implicit none 17#include "rtdb.fh" 18#include "dft2drv.fh" 19#include "dft3drv.fh" 20#include "stdio.fh" 21c !!! BGJ test 22#include "bgj.fh" 23c !!! BGJ test 24#include "errquit.fh" 25c 26 integer npert ! Number CPKS perturbations [input] 27 integer ipol ! Restricted or unrestricted [input] 28 integer nq ! Number of grid points [input] 29c 30 logical grad ! Whether gradient-corrected [input] 31 logical triplet ! Whether a triplet calculation is being done 32c 33c Current approximate perturbed spin densities and density gradients 34c These are overwritten with the XC matrix coefficients to save space 35c 36 double precision prho(nq,ipol,npert) 37 double precision pdelrho(nq,3,ipol,npert) 38c 39c Third derivatives of XC functional [input] 40c 41 double precision Amat3(nq,NCOL_AMAT3) 42 double precision Cmat3(nq,NCOL_CMAT3) 43c 44c GC second derivatives of XC functional [input] 45c 46 double precision Cmat2(nq,NCOL_CMAT2) 47c 48c Gradients of spin densities [input] 49c 50 double precision delrho(nq,3,ipol) 51c 52 integer ipert, n ! Loop indices 53c 54 double precision ptmp(10) 55c double precision t(2) 56 double precision pdra(3), pdrb(3) 57 double precision term_rrr, term_rrg, term_rgg, term_ggg 58 double precision term_rg,term_gg 59 double precision term_prho, term_pdelrho 60c 61 double precision term_rrg1, term_rrg2, term_rrg3 62 double precision term_rgg1, term_rgg2, term_rgg3 63 double precision term_rg1, term_rg2, term_rg3 64 double precision term_gg1, term_gg2, term_gg3 65c 66c Unrestricted specific variables 67c 68 double precision term_rrr1a, term_rrr1b 69 double precision term_rrg1a, term_rrg2a, term_rrg3a, term_rrg4a 70 double precision term_rrg1b, term_rrg2b, term_rrg3b, term_rrg4b 71 double precision term_rgg1a, term_rgg2a, term_rgg3a, term_rgg4a, 72 1 term_rgg5a, term_rgg6a 73 double precision term_rgg1b, term_rgg2b, term_rgg3b, term_rgg4b, 74 1 term_rgg5b, term_rgg6b 75 double precision term_ggg1a, term_ggg2a, term_ggg3a, term_ggg4a, 76 1 term_ggg5a, term_ggg6a, term_ggg7a, term_ggg8a 77 double precision term_ggg1b, term_ggg2b, term_ggg3b, term_ggg4b, 78 1 term_ggg5b, term_ggg6b, term_ggg7b, term_ggg8b 79 double precision term_rg1a, term_rg2a, term_rg3a 80 double precision term_rg1b, term_rg2b, term_rg3b 81 double precision term_gg1a, term_gg2a, term_gg3a, term_gg4a 82 double precision term_gg1b, term_gg2b, term_gg3b, term_gg4b 83 double precision term_pdelrho1a, term_pdelrho1b 84 double precision term_pdelrho2a, term_pdelrho2b 85c 86 character*32 pname 87c 88 pname = "xc_3rd_deriv: " 89c 90c preliminaries 91c 92 do n = 1,10 93 ptmp(n) = 0.d0 94 end do 95c 96 if (ipol.eq.1) then 97c 98c Since the total densities are evaluated in the restricted case, 99c scale them by a factor of 0.5 so that the correct CPKS matrices 100c will be produced. 101c 102 call dscal(nq*ipol*npert,0.5d0,prho,1) 103 if (grad) then 104 call dscal(nq*3*ipol*npert,0.5d0,pdelrho,1) 105 call dscal(nq*3*ipol,0.5d0,delrho,1) 106 endif 107 endif 108c 109 do ipert = 1, npert 110c 111c !!! Put in cutoffs here similar to xc_tabcd? !!! 112c 113 if (ipol.eq.2) then ! unrestricted 114c -------------------------------------------------------------- 115c Unrestricted case 116c (TDDFT excitation energy gradients) 117c (More stuff later) 118c -------------------------------------------------------------- 119c Daniel (3-26-13): For this case, the user needs to be aware that 120c spin contamination can result in very different gradients. This 121c results from the significant changes in both the excitation energy 122c and excitation vectors. 123 if (.not. grad) then ! local functionals 124 do n = 1, nq 125 ptmp(1) = prho(n,1,ipert)*prho(n,1,ipert) 126 ptmp(2) = prho(n,1,ipert)*prho(n,2,ipert) 127 ptmp(3) = prho(n,2,ipert)*prho(n,2,ipert) 128c Here, we multiply by the perturbed density twice, mimicing what was 129c done for the XC-kernel where the perturbed density is multiplied in 130c once. 131 prho(n,1,ipert) = Amat3(n,D3_RA_RA_RA)*ptmp(1) 132 1 + 2.0d0*Amat3(n,D3_RA_RA_RB)*ptmp(2) 133 2 + Amat3(n,D3_RA_RB_RB)*ptmp(3) 134c 135 prho(n,2,ipert) = Amat3(n,D3_RA_RA_RB)*ptmp(1) 136 1 + 2.0d0*Amat3(n,D3_RA_RB_RB)*ptmp(2) 137 2 + Amat3(n,D3_RB_RB_RB)*ptmp(3) 138 enddo 139 else ! gradient dependent functionals 140 do n = 1, nq 141c 142c Perturbed densities 143c 144 ptmp(1) = prho(n,1,ipert) ! prho_alpha 145c 146 ptmp(2) = prho(n,2,ipert) ! prho_beta 147c 148 ptmp(3) = prho(n,1,ipert) + prho(n,2,ipert) ! prho_alpha + prho_beta 149c 150c Products of the ground state density gradient and perturbed density 151c gradient 152c 153 ptmp(4) = delrho(n,1,1)*pdelrho(n,1,1,ipert) ! delrho_alpha * pdelrho_alpha 154 1 + delrho(n,2,1)*pdelrho(n,2,1,ipert) 155 2 + delrho(n,3,1)*pdelrho(n,3,1,ipert) 156c 157 ptmp(5) = delrho(n,1,1)*pdelrho(n,1,2,ipert) ! delrho_alpha * pdelrho_beta 158 1 + delrho(n,2,1)*pdelrho(n,2,2,ipert) 159 2 + delrho(n,3,1)*pdelrho(n,3,2,ipert) 160c 161 ptmp(6) = delrho(n,1,2)*pdelrho(n,1,1,ipert) ! delrho_beta * pdelrho_alpha 162 1 + delrho(n,2,2)*pdelrho(n,2,1,ipert) 163 2 + delrho(n,3,2)*pdelrho(n,3,1,ipert) 164c 165 ptmp(7) = delrho(n,1,2)*pdelrho(n,1,2,ipert) ! delrho_beta * pdelrho_beta 166 1 + delrho(n,2,2)*pdelrho(n,2,2,ipert) 167 2 + delrho(n,3,2)*pdelrho(n,3,2,ipert) 168c 169c Products of the perturbed density gradients 170c 171 ptmp(8) = pdelrho(n,1,1,ipert)*pdelrho(n,1,1,ipert) ! pdelrho_alpha * pdelrho_alpha 172 1 + pdelrho(n,2,1,ipert)*pdelrho(n,2,1,ipert) 173 2 + pdelrho(n,3,1,ipert)*pdelrho(n,3,1,ipert) 174c 175 ptmp(9) = pdelrho(n,1,1,ipert)*pdelrho(n,1,2,ipert) ! pdelrho_alpha * pdelrho_beta 176 1 + pdelrho(n,2,1,ipert)*pdelrho(n,2,2,ipert) ! OR 177 2 + pdelrho(n,3,1,ipert)*pdelrho(n,3,2,ipert) ! pdelrho_beta * pdelrho_alpha 178c 179 ptmp(10) = pdelrho(n,1,2,ipert)*pdelrho(n,1,2,ipert) ! pdelrho_beta * pdelrho_beta 180 1 + pdelrho(n,2,2,ipert)*pdelrho(n,2,2,ipert) 181 2 + pdelrho(n,3,2,ipert)*pdelrho(n,3,2,ipert) 182c ------------------------------------------------------------------- 183c Third derivative contributions: drdrdr 184c ------------------------------------------------------------------- 185 term_rrr1a = Amat3(n,D3_RA_RA_RA)*ptmp(1)*ptmp(1) 186 1 + 2.0d0*Amat3(n,D3_RA_RA_RB)*ptmp(1)*ptmp(2) 187 2 + Amat3(n,D3_RA_RB_RB)*ptmp(2)*ptmp(2) 188c 189 term_rrr1b = Amat3(n,D3_RA_RA_RB)*ptmp(1)*ptmp(1) 190 1 + 2.0d0*Amat3(n,D3_RA_RB_RB)*ptmp(1)*ptmp(2) 191 2 + Amat3(n,D3_RB_RB_RB)*ptmp(2)*ptmp(2) 192c ------------------------------------------------------------------- 193c Third derivative contributions: drdrdg 194c ------------------------------------------------------------------- 195c Alpha spin terms 196 term_rrg1a = 2.0d0*( 2.0d0*Cmat3(n,D3_RA_RA_GAA)*ptmp(4) 197 1 + Cmat3(n,D3_RA_RA_GAB) 198 2 *( ptmp(6) + ptmp(5) ) 199 3 + 2.0d0*Cmat3(n,D3_RA_RA_GBB) 200 4 *ptmp(7) ) 201c 202 term_rrg2a = 2.0d0*( 2.0d0*Cmat3(n,D3_RA_RB_GAA)*ptmp(4) 203 1 + Cmat3(n,D3_RA_RB_GAB) 204 2 *( ptmp(6) + ptmp(5) ) 205 3 + 2.0d0*Cmat3(n,D3_RA_RB_GBB) 206 4 *ptmp(7) ) 207c 208 term_rrg3a = 2.0d0*( Cmat3(n,D3_RA_RA_GAA) 209 1 *ptmp(1)*ptmp(1) 210 2 + 2.0d0*Cmat3(n,D3_RA_RB_GAA) 211 3 *ptmp(2)*ptmp(1) 212 4 + Cmat3(n,D3_RB_RB_GAA) 213 5 *ptmp(2)*ptmp(2) ) 214c 215 term_rrg4a = Cmat3(n,D3_RA_RA_GAB)*ptmp(1)*ptmp(1) 216 1 + 2.0d0*Cmat3(n,D3_RA_RB_GAB)*ptmp(2)*ptmp(1) 217 2 + Cmat3(n,D3_RB_RB_GAB)*ptmp(2)*ptmp(2) 218c 219c Beta spin terms 220 term_rrg1b = 2.0d0*( 2.0d0*Cmat3(n,D3_RA_RB_GAA)*ptmp(4) 221 1 + Cmat3(n,D3_RA_RB_GAB) 222 2 *( ptmp(6) + ptmp(5) ) 223 3 + 2.0d0*Cmat3(n,D3_RA_RB_GBB)*ptmp(7) ) 224c 225 term_rrg2b = 2.0d0*( 2.0d0*Cmat3(n,D3_RB_RB_GAA)*ptmp(4) 226 1 + Cmat3(n,D3_RB_RB_GAB) 227 2 *( ptmp(6) + ptmp(5) ) 228 3 + 2.0d0*Cmat3(n,D3_RB_RB_GBB)*ptmp(7) ) 229c 230 term_rrg3b = Cmat3(n,D3_RA_RA_GAB)*ptmp(1)*ptmp(1) 231 1 + 2.0d0*Cmat3(n,D3_RA_RB_GAB)*ptmp(2)*ptmp(1) 232 2 + Cmat3(n,D3_RB_RB_GAB)*ptmp(2)*ptmp(2) 233c 234 term_rrg4b = 2.0d0*( Cmat3(n,D3_RA_RA_GBB) 235 1 *ptmp(1)*ptmp(1) 236 2 + 2.0d0*Cmat3(n,D3_RA_RB_GBB) 237 3 *ptmp(2)*ptmp(1) 238 4 + Cmat3(n,D3_RB_RB_GBB) 239 5 *ptmp(2)*ptmp(2) ) 240c ------------------------------------------------------------------- 241c Third derivative contributions: drdgdg 242c ------------------------------------------------------------------- 243c Alpha spin terms 244 term_rgg1a = 2.0d0*( 2.0d0*Cmat3(n,D3_RA_GAA_GAA)*ptmp(4) 245 1 + Cmat3(n,D3_RA_GAA_GAB) 246 2 *( ptmp(6) + 2.0d0*ptmp(5) ) 247 3 + 4.0d0*Cmat3(n,D3_RA_GAA_GBB) 248 4 *ptmp(7) ) 249c 250 term_rgg2a = 2.0d0*Cmat3(n,D3_RA_GAA_GAB)*ptmp(4) 251 1 + Cmat3(n,D3_RA_GAB_GAB) 252 2 *( ptmp(6) + 2.0d0*ptmp(5) ) 253 3 + 4.0d0*Cmat3(n,D3_RA_GAB_GBB)*ptmp(7) 254c 255 term_rgg3a = 2.0d0*( 2.0d0*Cmat3(n,D3_RA_GBB_GBB)*ptmp(7) 256 1 + Cmat3(n,D3_RA_GAB_GBB)*ptmp(5) ) 257c 258 term_rgg4a = 2.0d0*Cmat3(n,D3_RA_GAB_GBB)*ptmp(7) 259 1 + Cmat3(n,D3_RA_GAB_GAB)*ptmp(5) 260c 261 term_rgg5a = 4.0d0*( 2.0d0*Cmat3(n,D3_RA_GAA_GAA)*ptmp(4) 262 1 + Cmat3(n,D3_RA_GAA_GAB) 263 2 *( ptmp(6) + ptmp(5) ) 264 3 + 2.0d0*Cmat3(n,D3_RA_GAA_GBB) 265 4 *ptmp(7) )*ptmp(1) 266 5 + 4.0d0*( 2.0d0*Cmat3(n,D3_RB_GAA_GAA)*ptmp(4) 267 6 + Cmat3(n,D3_RB_GAA_GAB) 268 7 *( ptmp(6) + ptmp(5) ) 269 8 + 2.0d0*Cmat3(n,D3_RB_GAA_GBB) 270 9 *ptmp(7) )*ptmp(2) 271c 272 term_rgg6a = 2.0d0*( 2.0d0*Cmat3(n,D3_RA_GAA_GAB)*ptmp(4) 273 1 + Cmat3(n,D3_RA_GAB_GAB) 274 2 *( ptmp(6) + ptmp(5) ) 275 3 + 2.0d0*Cmat3(n,D3_RA_GAB_GBB) 276 4 *ptmp(7) )*ptmp(1) 277 5 + 2.0d0*( 2.0d0*Cmat3(n,D3_RB_GAA_GAB)*ptmp(4) 278 6 + Cmat3(n,D3_RB_GAB_GAB) 279 7 *( ptmp(6) + ptmp(5) ) 280 8 + 2.0d0*Cmat3(n,D3_RB_GAB_GBB) 281 9 *ptmp(7) )*ptmp(2) 282c 283c Beta spin terms 284 term_rgg1b = 2.0d0*( 2.0d0*Cmat3(n,D3_RB_GAA_GAA)*ptmp(4) 285 1 + Cmat3(n,D3_RB_GAA_GAB) 286 2 *( ptmp(6) + 2.0d0*ptmp(5) ) 287 3 + 4.0d0*Cmat3(n,D3_RB_GAA_GBB) 288 4 *ptmp(7) ) 289c 290 term_rgg2b = 2.0d0*Cmat3(n,D3_RB_GAA_GAB)*ptmp(4) 291 1 + Cmat3(n,D3_RB_GAB_GAB) 292 2 *( ptmp(6) + 2.0d0*ptmp(5) ) 293 3 + 4.0d0*Cmat3(n,D3_RB_GAB_GBB)*ptmp(7) 294c 295 term_rgg3b = 2.0d0*( 2.0d0*Cmat3(n,D3_RB_GBB_GBB)*ptmp(7) 296 1 + Cmat3(n,D3_RB_GAB_GBB)*ptmp(5) ) 297c 298 term_rgg4b = 2.0d0*Cmat3(n,D3_RB_GAB_GBB)*ptmp(7) 299 1 + Cmat3(n,D3_RB_GAB_GAB)*ptmp(5) 300c 301 term_rgg5b = 4.0d0*( 2.0d0*Cmat3(n,D3_RA_GAA_GBB)*ptmp(4) 302 1 + Cmat3(n,D3_RA_GAB_GBB) 303 2 *( ptmp(6) + ptmp(5) ) 304 3 + 2.0d0*Cmat3(n,D3_RA_GBB_GBB) 305 4 *ptmp(7) )*ptmp(1) 306 5 + 4.0d0*( 2.0d0*Cmat3(n,D3_RB_GAA_GBB)*ptmp(4) 307 6 + Cmat3(n,D3_RB_GAB_GBB) 308 7 *( ptmp(6) + ptmp(5) ) 309 8 + 2.0d0*Cmat3(n,D3_RB_GBB_GBB) 310 9 *ptmp(7) )*ptmp(2) 311c 312 term_rgg6b = 2.0d0*( 2.0d0*Cmat3(n,D3_RA_GAA_GAB)*ptmp(4) 313 1 + Cmat3(n,D3_RA_GAB_GAB) 314 2 *( ptmp(6) + ptmp(5) ) 315 3 + 2.0d0*Cmat3(n,D3_RA_GAB_GBB) 316 4 *ptmp(7) )*ptmp(1) 317 5 + 2.0d0*( 2.0d0*Cmat3(n,D3_RB_GAA_GAB)*ptmp(4) 318 6 + Cmat3(n,D3_RB_GAB_GAB) 319 7 *( ptmp(6) + ptmp(5) ) 320 8 + 2.0d0*Cmat3(n,D3_RB_GAB_GBB) 321 9 *ptmp(7) )*ptmp(2) 322c term_rgg6b = term_rgg6a 323c ------------------------------------------------------------------- 324c Third derivative contributions: dgdgdg 325c ------------------------------------------------------------------- 326c Alpha spin terms 327 term_ggg1a = 4.0d0*( 2.0d0*Cmat3(n,D3_GAA_GAA_GAA) 328 1 *ptmp(4) 329 2 + Cmat3(n,D3_GAA_GAA_GAB) 330 3 *( ptmp(6) + 2.0d0*ptmp(5) ) 331 4 + 4.0d0*Cmat3(n,D3_GAA_GAA_GBB) 332 5 *ptmp(7) ) 333c 334 term_ggg2a = 2.0d0*( 2.0d0*Cmat3(n,D3_GAA_GAA_GAB) 335 1 *ptmp(4) 336 2 + Cmat3(n,D3_GAA_GAB_GAB) 337 3 *( ptmp(6) + 2.0d0*ptmp(5) ) 338 4 + 4.0d0*Cmat3(n,D3_GAA_GAB_GBB) 339 5 *ptmp(7) ) 340c 341 term_ggg3a = 4.0d0*( 2.0d0*Cmat3(n,D3_GAA_GBB_GBB) 342 1 *ptmp(7) 343 2 + Cmat3(n,D3_GAA_GAB_GBB)*ptmp(5) ) 344c 345 term_ggg4a = 2.0d0*( 2.0d0*Cmat3(n,D3_GAA_GAB_GBB) 346 1 *ptmp(7) 347 2 + Cmat3(n,D3_GAA_GAB_GAB)*ptmp(5) ) 348c 349 term_ggg5a = 2.0d0*( 2.0d0*Cmat3(n,D3_GAA_GAA_GAB) 350 1 *ptmp(4) 351 2 + Cmat3(n,D3_GAA_GAB_GAB) 352 3 *( ptmp(6) + 2.0d0*ptmp(5) ) 353 4 + 4.0d0*Cmat3(n,D3_GAA_GAB_GBB) 354 5 *ptmp(7) ) 355c 356 term_ggg6a = 2.0d0*Cmat3(n,D3_GAA_GAB_GAB)*ptmp(4) 357 1 + Cmat3(n,D3_GAB_GAB_GAB) 358 2 *( ptmp(6) + 2.0d0*ptmp(5) ) 359 3 + 4.0d0*Cmat3(n,D3_GAB_GAB_GBB)*ptmp(7) 360c 361 term_ggg7a = 2.0d0*( 2.0d0*Cmat3(n,D3_GAB_GBB_GBB) 362 1 *ptmp(7) 363 2 + Cmat3(n,D3_GAB_GAB_GBB) 364 3 *ptmp(5) ) 365c 366 term_ggg8a = 2.0d0*Cmat3(n,D3_GAB_GAB_GBB)*ptmp(7) 367 1 + Cmat3(n,D3_GAB_GAB_GAB)*ptmp(5) 368c 369c Beta spin terms 370 term_ggg1b = 2.0d0*( 2.0d0*Cmat3(n,D3_GAA_GAA_GAB) 371 1 *ptmp(4) 372 2 + Cmat3(n,D3_GAA_GAB_GAB) 373 3 *( ptmp(6) + 2.0d0*ptmp(5) ) 374 4 + 4.0d0*Cmat3(n,D3_GAA_GAB_GBB) 375 5 *ptmp(7) ) 376c 377 term_ggg2b = 2.0d0*Cmat3(n,D3_GAA_GAB_GAB)*ptmp(4) 378 1 + Cmat3(n,D3_GAB_GAB_GAB) 379 2 *( ptmp(6) + 2.0d0*ptmp(5) ) 380 3 + 4.0d0*Cmat3(n,D3_GAB_GAB_GBB)*ptmp(7) 381c 382 term_ggg3b = 2.0d0*( 2.0d0*Cmat3(n,D3_GAB_GBB_GBB) 383 1 *ptmp(7) 384 2 + Cmat3(n,D3_GAB_GAB_GBB)*ptmp(5) ) 385c 386 term_ggg4b = 2.0d0*Cmat3(n,D3_GAB_GAB_GBB)*ptmp(7) 387 1 + Cmat3(n,D3_GAB_GAB_GAB)*ptmp(5) 388c 389 term_ggg5b = 4.0d0*( 2.0d0*Cmat3(n,D3_GAA_GAA_GBB) 390 1 *ptmp(4) 391 2 + Cmat3(n,D3_GAA_GAB_GBB) 392 3 *( ptmp(6) + 2.0d0*ptmp(5) ) 393 4 + 4.0d0*Cmat3(n,D3_GAA_GBB_GBB) 394 5 *ptmp(7) ) 395c 396 term_ggg6b = 2.0d0*( 2.0d0*Cmat3(n,D3_GAA_GAB_GBB) 397 1 *ptmp(4) 398 2 + Cmat3(n,D3_GAB_GAB_GBB) 399 3 *( ptmp(6) + 2.0d0*ptmp(5) ) 400 4 + 4.0d0*Cmat3(n,D3_GAB_GBB_GBB) 401 5 *ptmp(7) ) 402c 403 term_ggg7b = 4.0d0*( 2.0d0*Cmat3(n,D3_GBB_GBB_GBB) 404 1 *ptmp(7) 405 2 + Cmat3(n,D3_GAB_GBB_GBB)*ptmp(5) ) 406c 407 term_ggg8b = 2.0d0*( 2.0d0*Cmat3(n,D3_GAB_GBB_GBB) 408 1 *ptmp(7) 409 2 + Cmat3(n,D3_GAB_GAB_GBB)*ptmp(5) ) 410c ------------------------------------------------------------------- 411c Third derivative contributions: drdg 412c ------------------------------------------------------------------- 413c Alpha spin terms 414 term_rg1a = 2.0d0*( Cmat2(n,D2_RA_GAA)*ptmp(8) 415 1 + Cmat2(n,D2_RA_GAB)*ptmp(9) 416 2 + Cmat2(n,D2_RA_GBB)*ptmp(10) ) 417c 418 term_rg2a = 4.0d0*( Cmat2(n,D2_RA_GAA)*ptmp(1) 419 1 + Cmat2(n,D2_RB_GAA)*ptmp(2) ) 420c 421 term_rg3a = 2.0d0*( Cmat2(n,D2_RA_GAB)*ptmp(1) 422 1 + Cmat2(n,D2_RB_GAB)*ptmp(2) ) 423c 424c Beta spin terms 425 term_rg1b = 2.0d0*( Cmat2(n,D2_RB_GAA)*ptmp(8) 426 1 + Cmat2(n,D2_RB_GAB)*ptmp(9) 427 2 + Cmat2(n,D2_RB_GBB)*ptmp(10) ) 428c 429 term_rg2b = 2.0d0*( Cmat2(n,D2_RA_GAB)*ptmp(1) 430 1 + Cmat2(n,D2_RB_GAB)*ptmp(2) ) 431c 432 term_rg3b = 4.0d0*( Cmat2(n,D2_RA_GBB)*ptmp(1) 433 1 + Cmat2(n,D2_RB_GBB)*ptmp(2) ) 434c ------------------------------------------------------------------- 435c Third derivative contributions: dgdg 436c ------------------------------------------------------------------- 437c Alpha spin terms 438 term_gg1a = 4.0d0*( 2.0d0*Cmat2(n,D2_GAA_GAA)*ptmp(4) 439 1 + Cmat2(n,D2_GAA_GAB) 440 2 *( ptmp(6) + ptmp(5) ) 441 3 + 2.0d0*Cmat2(n,D2_GAA_GBB)*ptmp(7) ) 442c 443 term_gg2a = 4.0d0*( Cmat2(n,D2_GAA_GAA)*ptmp(8) 444 1 + Cmat2(n,D2_GAA_GAB)*ptmp(9) 445 2 + Cmat2(n,D2_GAA_GBB)*ptmp(10) ) 446c 447 term_gg3a = 2.0d0*( 2.0d0*Cmat2(n,D2_GAA_GAB)*ptmp(4) 448 1 + Cmat2(n,D2_GAB_GAB) 449 2 *( ptmp(6) + ptmp(5) ) 450 3 + 2.0d0*Cmat2(n,D2_GAB_GBB)*ptmp(7) ) 451c 452 term_gg4a = 2.0d0*( Cmat2(n,D2_GAA_GAB)*ptmp(8) 453 1 + Cmat2(n,D2_GAB_GAB)*ptmp(9) 454 2 + Cmat2(n,D2_GAB_GBB)*ptmp(10) ) 455c 456c Beta spin terms 457 term_gg1b = 2.0d0*( 2.0d0*Cmat2(n,D2_GAA_GAB)*ptmp(4) 458 1 + Cmat2(n,D2_GAB_GAB) 459 2 *( ptmp(6) + ptmp(5) ) 460 3 + 2.0d0*Cmat2(n,D2_GAB_GBB)*ptmp(7) ) 461c 462 term_gg2b = 2.0d0*( Cmat2(n,D2_GAA_GAB)*ptmp(8) 463 1 + Cmat2(n,D2_GAB_GAB)*ptmp(9) 464 2 + Cmat2(n,D2_GAB_GBB)*ptmp(10) ) 465c 466 term_gg3b = 4.0d0*( 2.0d0*Cmat2(n,D2_GAA_GBB)*ptmp(4) 467 1 + Cmat2(n,D2_GAB_GBB) 468 2 *( ptmp(6) + ptmp(5) ) 469 3 + 2.0d0*Cmat2(n,D2_GBB_GBB)*ptmp(7) ) 470c 471 term_gg4b = 4.0d0*( Cmat2(n,D2_GAA_GBB)*ptmp(8) 472 1 + Cmat2(n,D2_GAB_GBB)*ptmp(9) 473 2 + Cmat2(n,D2_GBB_GBB)*ptmp(10) ) 474c ------------------------------------------------------------------- 475c Terms that do not involve the gradient with respect to the basis 476c functions. 477c ------------------------------------------------------------------- 478c 479c Recall that: 480c ptmp(1) = prho_alpha 481c ptmp(2) = prho_beta 482c ptmp(3) = prho_alpha + prho_beta 483c ptmp(4) = delrho_alpha * pdelrho_alpha 484c ptmp(5) = delrho_alpha * pdelrho_beta 485c ptmp(6) = delrho_beta * pdelrho_alpha 486c ptmp(7) = delrho_beta * pdelrho_beta 487c ptmp(8) = pdelrho_alpha * pdelrho_alpha 488c ptmp(9) = pdelrho_alpha * pdelrho_beta = pdelrho_beta * pdelrho_alpha 489c ptmp(10) = pdelrho_beta * pdelrho_beta 490c 491 prho(n,1,ipert) = term_rrr1a 492 1 + term_rrg1a*ptmp(1) + term_rrg2a*ptmp(2) 493 2 + term_rgg1a*ptmp(4) + term_rgg2a*ptmp(6) 494 3 + term_rgg3a*ptmp(7) + term_rgg4a*ptmp(5) 495 4 + term_rg1a 496c 497 prho(n,2,ipert) = term_rrr1b 498 1 + term_rrg1b*ptmp(1) + term_rrg2b*ptmp(2) 499 2 + term_rgg1b*ptmp(4) + term_rgg2b*ptmp(6) 500 3 + term_rgg3b*ptmp(7) + term_rgg4b*ptmp(5) 501 4 + term_rg1b 502c ------------------------------------------------------------------- 503c Terms that involve the gradient with respect to the basis 504c functions. 505c ------------------------------------------------------------------- 506c These terms are multiplied by the ground state density gradient 507c Multiplied by delrho_alpha (alpha coefficient) 508 term_pdelrho1a = term_rrg3a 509 1 + term_rgg5a 510 2 + term_ggg1a*ptmp(4) + term_ggg2a*ptmp(6) 511 3 + term_ggg3a*ptmp(7) + term_ggg4a*ptmp(5) 512 4 + term_gg2a 513c Multiplied by delrho_beta (alpha coefficient) 514 term_pdelrho1b = term_rrg4a 515 1 + term_rgg6a 516 2 + term_ggg5a*ptmp(4) + term_ggg6a*ptmp(6) 517 3 + term_ggg7a*ptmp(7) + term_ggg8a*ptmp(5) 518 4 + term_gg4a 519c Multiplied by delrho_alpha (beta coefficient) 520 term_pdelrho2a = term_rrg3b 521 1 + term_rgg5b 522 2 + term_ggg1b*ptmp(4) + term_ggg2b*ptmp(6) 523 3 + term_ggg3b*ptmp(7) + term_ggg4b*ptmp(5) 524 4 + term_gg2b 525c Multiplied by delrho_beta (beta coefficient) 526 term_pdelrho2b = term_rrg4b 527 1 + term_rgg6b 528 2 + term_ggg5b*ptmp(4) + term_ggg6b*ptmp(6) 529 3 + term_ggg7b*ptmp(7) + term_ggg8b*ptmp(5) 530 4 + term_gg4b 531c ------------------------------------------------------------------- 532c Daniel (3-19-13): These can be written more efficiently 533 pdra(1) = pdelrho(n,1,1,ipert) 534 pdra(2) = pdelrho(n,2,1,ipert) 535 pdra(3) = pdelrho(n,3,1,ipert) 536 pdrb(1) = pdelrho(n,1,2,ipert) 537 pdrb(2) = pdelrho(n,2,2,ipert) 538 pdrb(3) = pdelrho(n,3,2,ipert) 539c 540 pdelrho(n,1,1,ipert) = term_rg2a*pdra(1) 541 1 + term_rg3a*pdrb(1) 542 2 + term_gg1a*pdra(1) 543 3 + term_gg3a*pdrb(1) 544 4 + term_pdelrho1a*delrho(n,1,1) 545 5 + term_pdelrho1b*delrho(n,1,2) 546c 547 pdelrho(n,2,1,ipert) = term_rg2a*pdra(2) 548 1 + term_rg3a*pdrb(2) 549 2 + term_gg1a*pdra(2) 550 3 + term_gg3a*pdrb(2) 551 4 + term_pdelrho1a*delrho(n,2,1) 552 5 + term_pdelrho1b*delrho(n,2,2) 553c 554 pdelrho(n,3,1,ipert) = term_rg2a*pdra(3) 555 1 + term_rg3a*pdrb(3) 556 2 + term_gg1a*pdra(3) 557 3 + term_gg3a*pdrb(3) 558 4 + term_pdelrho1a*delrho(n,3,1) 559 5 + term_pdelrho1b*delrho(n,3,2) 560c 561 pdelrho(n,1,2,ipert) = term_rg2b*pdra(1) 562 1 + term_rg3b*pdrb(1) 563 2 + term_gg1b*pdra(1) 564 3 + term_gg3b*pdrb(1) 565 4 + term_pdelrho2a*delrho(n,1,1) 566 5 + term_pdelrho2b*delrho(n,1,2) 567c 568 pdelrho(n,2,2,ipert) = term_rg2b*pdra(2) 569 1 + term_rg3b*pdrb(2) 570 2 + term_gg1b*pdra(2) 571 3 + term_gg3b*pdrb(2) 572 4 + term_pdelrho2a*delrho(n,2,1) 573 5 + term_pdelrho2b*delrho(n,2,2) 574c 575 pdelrho(n,3,2,ipert) = term_rg2b*pdra(3) 576 1 + term_rg3b*pdrb(3) 577 2 + term_gg1b*pdra(3) 578 3 + term_gg3b*pdrb(3) 579 4 + term_pdelrho2a*delrho(n,3,1) 580 5 + term_pdelrho2b*delrho(n,3,2) 581c 582 enddo 583 endif 584 else if (triplet) then 585c -------------------------------------------------------------- 586c Restricted Triplet case 587c (TDDFT excitation energy gradients) 588c (More stuff later) 589c -------------------------------------------------------------- 590 if (.not. grad) then ! local functionals 591 do n = 1, nq 592c Note that there is a cancellation of terms here (see comment at 593c right). 594 term_rrr = Amat3(n,D3_RA_RA_RA) ! rarara + rararb 595 1 - Amat3(n,D3_RA_RB_RB) ! - rararb - rarbrb 596c Here, we multiply by the perturbed density twice, mimicing what was 597c done for the XC-kernel where the perturbed density is multiplied in 598c once. 599 prho(n,1,ipert) = 600 1 term_rrr*prho(n,1,ipert)*prho(n,1,ipert) 601 enddo 602 else ! gradient dependent functionals 603 do n = 1, nq 604 ptmp(1) = prho(n,1,ipert) ! perturbed density 605c 606c Here we collect the various dot products needed to construct the 607c total functional derivative for a GGA. 608 ptmp(2) = delrho(n,1,1)*pdelrho(n,1,1,ipert) ! delrho*del-perturbed density 609 1 + delrho(n,2,1)*pdelrho(n,2,1,ipert) 610 2 + delrho(n,3,1)*pdelrho(n,3,1,ipert) 611c 612 ptmp(3) = pdelrho(n,1,1,ipert)*pdelrho(n,1,1,ipert) ! del-perturbed density*del-perturbed density 613 1 + pdelrho(n,2,1,ipert)*pdelrho(n,2,1,ipert) 614 2 + pdelrho(n,3,1,ipert)*pdelrho(n,3,1,ipert) 615c ------------------------------------------------------------------- 616c Third derivative contributions: drdrdr 617c ------------------------------------------------------------------- 618 term_rrr = Amat3(n,D3_RA_RA_RA) - Amat3(n,D3_RA_RB_RB) 619c ------------------------------------------------------------------- 620c Third derivative contributions: drdrdg 621c ------------------------------------------------------------------- 622 term_rrg1 = 2.0d0*( 2.0d0*Cmat3(n,D3_RA_RA_GAA) 623 1 + Cmat3(n,D3_RA_RA_GAB) 624 2 - 2.0d0*Cmat3(n,D3_RA_RB_GBB) 625 3 - Cmat3(n,D3_RA_RB_GAB) ) 626c 627 term_rrg2 = 2.0d0*Cmat3(n,D3_RA_RA_GAA) 628 1 + Cmat3(n,D3_RA_RA_GAB) 629 2 - 2.0d0*Cmat3(n,D3_RA_RB_GBB) 630 3 - Cmat3(n,D3_RA_RB_GAB) 631c ------------------------------------------------------------------- 632c Third derivative contributions: drdgdg 633c ------------------------------------------------------------------- 634 term_rgg1 = 4.0d0*( Cmat3(n,D3_RA_GAA_GAA) 635 1 + Cmat3(n,D3_RA_GAA_GAB) 636 2 - Cmat3(n,D3_RA_GBB_GBB) 637 3 - Cmat3(n,D3_RA_GAB_GBB) ) 638c 639 term_rgg2 = 8.0d0*( Cmat3(n,D3_RA_GAA_GAA) 640 1 + Cmat3(n,D3_RA_GAA_GAB) 641 2 - Cmat3(n,D3_RA_GAA_GBB) 642 3 - Cmat3(n,D3_RA_GAB_GBB) ) 643c ------------------------------------------------------------------- 644c Third derivative contributions: dgdgdg 645c ------------------------------------------------------------------- 646 term_ggg = 8.0d0*( Cmat3(n,D3_GAA_GAA_GAA) 647 1 + Cmat3(n,D3_GAA_GAA_GAB) 648 3 - Cmat3(n,D3_GAA_GAA_GBB) 649 4 - Cmat3(n,D3_GAA_GAB_GBB) ) 650c ------------------------------------------------------------------- 651c Second derivative contributions: drdg 652c ------------------------------------------------------------------- 653 term_rg1 = 2.0d0*( Cmat2(n,D2_RA_GAA) 654 1 - Cmat2(n,D2_RA_GBB) ) 655c 656 term_rg2 = 4.0d0*( Cmat2(n,D2_RA_GAA) 657 1 - Cmat2(n,D2_RA_GBB) ) 658c ------------------------------------------------------------------- 659c Second derivative contributions: dgdg 660c ------------------------------------------------------------------- 661 term_gg1 = 8.0d0*( Cmat2(n,D2_GAA_GAA) 662 1 - Cmat2(n,D2_GAA_GBB) ) 663c 664 term_gg2 = 4.0d0*( Cmat2(n,D2_GAA_GAA) 665 1 - Cmat2(n,D2_GAA_GBB) ) 666c ------------------------------------------------------------------- 667c Terms that do not involve the gradient with respect to the basis 668c functions. 669c ------------------------------------------------------------------- 670c 671c Recall that: 672c ptmp(1) = prho 673c ptmp(2) = pdelrho*delrho 674c ptmp(3) = pdelrho*pdelrho 675c 676 term_prho = term_rrr*ptmp(1)*ptmp(1) 677 1 + term_rrg1*ptmp(2)*ptmp(1) 678 2 + term_rgg1*ptmp(2)*ptmp(2) 679 3 + term_rg1*ptmp(3) 680c ------------------------------------------------------------------- 681c Terms that involve the gradient with respect to the basis 682c functions. 683c ------------------------------------------------------------------- 684 term_pdelrho = term_rrg2*ptmp(1)*ptmp(1) 685 1 + term_rgg2*ptmp(2)*ptmp(1) 686 2 + term_ggg*ptmp(2)*ptmp(2) 687 3 + term_gg2*ptmp(3) 688c ------------------------------------------------------------------- 689 prho(n,1,ipert) = term_prho 690c 691 pdelrho(n,1,1,ipert) = 692 1 term_rg2*pdelrho(n,1,1,ipert)*ptmp(1) 693 2 + term_gg1*pdelrho(n,1,1,ipert)*ptmp(2) 694 3 + term_pdelrho*delrho(n,1,1) 695c 696 pdelrho(n,2,1,ipert) = 697 1 term_rg2*pdelrho(n,2,1,ipert)*ptmp(1) 698 2 + term_gg1*pdelrho(n,2,1,ipert)*ptmp(2) 699 2 + term_pdelrho*delrho(n,2,1) 700c 701 pdelrho(n,3,1,ipert) = 702 1 term_rg2*pdelrho(n,3,1,ipert)*ptmp(1) 703 2 + term_gg1*pdelrho(n,3,1,ipert)*ptmp(2) 704 3 + term_pdelrho*delrho(n,3,1) 705 enddo 706 endif 707 else ! singlet case 708c -------------------------------------------------------------- 709c Restricted case 710c (TDDFT excitation energy gradients) 711c (More stuff later) 712c -------------------------------------------------------------- 713 if (.not. grad) then ! local functionals 714 do n = 1, nq 715 term_rrr = Amat3(n,D3_RA_RA_RA) ! rarara + 2*rararb + rarbrb 716 1 + 2.0d0*Amat3(n,D3_RA_RA_RB) 717 2 + Amat3(n,D3_RA_RB_RB) 718c Here, we multiply by the perturbed density twice, mimicing what was 719c done for the XC-kernel where the perturbed density is multiplied in 720c once. 721 prho(n,1,ipert) = 722 1 term_rrr*prho(n,1,ipert)*prho(n,1,ipert) 723 enddo 724 else ! gradient dependent functionals 725 do n = 1, nq 726 ptmp(1) = prho(n,1,ipert) ! perturbed density 727c 728c Here we collect the various dot products needed to construct the 729c total functional derivative for a GGA. 730 ptmp(2) = delrho(n,1,1)*pdelrho(n,1,1,ipert) ! delrho*del-perturbed density 731 1 + delrho(n,2,1)*pdelrho(n,2,1,ipert) 732 2 + delrho(n,3,1)*pdelrho(n,3,1,ipert) 733c 734 ptmp(3) = pdelrho(n,1,1,ipert)*pdelrho(n,1,1,ipert) ! del-perturbed density*del-perturbed density 735 1 + pdelrho(n,2,1,ipert)*pdelrho(n,2,1,ipert) 736 2 + pdelrho(n,3,1,ipert)*pdelrho(n,3,1,ipert) 737c ------------------------------------------------------------------- 738c Third derivative contributions: drdrdr 739c ------------------------------------------------------------------- 740c Correct terms: term_rrr 741 term_rrr = Amat3(n,D3_RA_RA_RA) 742 1 + 2.0d0*Amat3(n,D3_RA_RA_RB) 743 2 + Amat3(n,D3_RA_RB_RB) 744c ------------------------------------------------------------------- 745c Third derivative contributions: drdrdg 746c ------------------------------------------------------------------- 747c Correct terms: term_rrg1, term_rrg2, term_rrg3 748 term_rrg1 = 2.0d0*( Cmat3(n,D3_RA_RA_GAA) 749 1 + 1.5d0*Cmat3(n,D3_RA_RA_GAB) 750 2 + 2.0d0*Cmat3(n,D3_RA_RA_GBB) 751 3 + Cmat3(n,D3_RA_RB_GBB) 752 4 + 0.5d0*Cmat3(n,D3_RA_RB_GAB)) 753c 754 term_rrg2 = 2.0d0*( Cmat3(n,D3_RA_RA_GAA) 755 1 + 0.5d0*Cmat3(n,D3_RA_RA_GAB) 756 2 + 2.0d0*Cmat3(n,D3_RA_RB_GAA) 757 3 + 1.5d0*Cmat3(n,D3_RA_RB_GAB) 758 4 + Cmat3(n,D3_RA_RB_GBB)) 759c 760 term_rrg3 = 2.0d0*( Cmat3(n,D3_RA_RA_GAA) 761 1 + Cmat3(n,D3_RA_RA_GAB) 762 2 + 2.0d0*Cmat3(n,D3_RA_RB_GAA) 763 3 + Cmat3(n,D3_RA_RB_GAB) 764 4 + Cmat3(n,D3_RA_RA_GBB)) 765c ------------------------------------------------------------------- 766c Third derivative contributions: drdgdg 767c ------------------------------------------------------------------- 768c Correct terms: term_rgg1, term_rgg2, term_rgg3 769 term_rgg1 = 4.0d0*( Cmat3(n,D3_RA_GAA_GAA) 770 1 + 2.0d0*Cmat3(n,D3_RA_GAA_GAB) 771 2 + Cmat3(n,D3_RA_GAB_GAB) 772 3 + 2.0d0*Cmat3(n,D3_RA_GAA_GBB) 773 4 + 2.0d0*Cmat3(n,D3_RA_GAB_GBB) 774 5 + Cmat3(n,D3_RA_GBB_GBB)) 775c 776 term_rgg2 = 4.0d0*( Cmat3(n,D3_RA_GAA_GAA) 777 1 + 2.0d0*Cmat3(n,D3_RA_GAA_GAB) 778 2 + Cmat3(n,D3_RA_GAB_GAB) 779 3 + 3.0d0*Cmat3(n,D3_RA_GAA_GBB) 780 4 + 2.0d0*Cmat3(n,D3_RA_GAB_GBB)) 781c 782 term_rgg3 = 4.0d0*( Cmat3(n,D3_RA_GAA_GAA) 783 1 + 2.0d0*Cmat3(n,D3_RA_GAA_GAB) 784 2 + Cmat3(n,D3_RA_GAB_GAB) 785 3 + 2.0d0*Cmat3(n,D3_RA_GBB_GBB) 786 4 + 2.0d0*Cmat3(n,D3_RA_GAB_GBB) 787 5 + Cmat3(n,D3_RA_GAA_GBB)) 788c ------------------------------------------------------------------- 789c Third derivative contributions: dgdgdg 790c ------------------------------------------------------------------- 791c Correct terms: term_ggg 792 term_ggg = 4.0d0*( 2.0d0*Cmat3(n,D3_GAA_GAA_GAA) 793 1 + 6.0d0*Cmat3(n,D3_GAA_GAA_GAB) 794 2 + 6.0d0*Cmat3(n,D3_GAA_GAB_GAB) 795 3 + Cmat3(n,D3_GAB_GAB_GAB) 796 4 + 4.0d0*Cmat3(n,D3_GAA_GAA_GBB) 797 5 + 6.0d0*Cmat3(n,D3_GAA_GAB_GBB) 798 6 + 2.0d0*Cmat3(n,D3_GAA_GBB_GBB)) 799c ------------------------------------------------------------------- 800c Second derivative contributions: drdg 801c ------------------------------------------------------------------- 802c Correct terms: term_rg1, term_rg2, term_rg3 803 term_rg1 = 2.0d0*( Cmat2(n,D2_RA_GAA) 804 1 + Cmat2(n,D2_RA_GAB) 805 2 + Cmat2(n,D2_RA_GBB)) 806c 807 term_rg2 = 2.0d0*( Cmat2(n,D2_RA_GAA) 808 1 + 1.5d0*Cmat2(n,D2_RA_GAB)) 809c 810 term_rg3 = 2.0d0*( Cmat2(n,D2_RA_GAA) 811 1 + 0.5d0*Cmat2(n,D2_RA_GAB) 812 2 + 2.0d0*Cmat2(n,D2_RA_GBB)) 813c ------------------------------------------------------------------- 814c Second derivative contributions: dgdg 815c ------------------------------------------------------------------- 816c Correct terms: term_gg1, term_gg2, term_gg3 817 term_gg1 = 4.0d0*Cmat2(n,D2_GAA_GAA) 818 1 + 4.0d0*Cmat2(n,D2_GAA_GBB) 819 2 + 8.0d0*Cmat2(n,D2_GAA_GAB) 820 3 + 2.0d0*Cmat2(n,D2_GAB_GAB) 821c 822 term_gg2 = 4.0d0*Cmat2(n,D2_GAA_GAA) 823 1 + 8.0d0*Cmat2(n,D2_GAA_GAB) 824 2 + 3.0d0*Cmat2(n,D2_GAB_GAB) 825c 826 term_gg3 = 4.0d0*Cmat2(n,D2_GAA_GAA) 827 1 + 8.0d0*Cmat2(n,D2_GAA_GBB) 828 2 + 8.0d0*Cmat2(n,D2_GAA_GAB) 829 3 + Cmat2(n,D2_GAB_GAB) 830c ------------------------------------------------------------------- 831c Terms that do not involve the gradient with respect to the basis 832c functions. 833c ------------------------------------------------------------------- 834c 835c Recall that: 836c ptmp(1) = prho 837c ptmp(2) = pdelrho*delrho 838c ptmp(3) = pdelrho*pdelrho 839c 840 term_prho = term_rrr*ptmp(1)*ptmp(1) 841 1 + (term_rrg1 + term_rrg2)*ptmp(2)*ptmp(1) 842 2 + term_rgg1*ptmp(2)*ptmp(2) 843 3 + term_rg1*ptmp(3) 844c ------------------------------------------------------------------- 845c Terms that involve the gradient with respect to the basis 846c functions. 847c ------------------------------------------------------------------- 848 term_pdelrho = term_rrg3*ptmp(1)*ptmp(1) 849 1 + (term_rgg2 + term_rgg3)*ptmp(2)*ptmp(1) 850 2 + term_ggg*ptmp(2)*ptmp(2) 851 3 + term_gg1*ptmp(3) 852c ------------------------------------------------------------------- 853 prho(n,1,ipert) = term_prho 854c 855 pdelrho(n,1,1,ipert) = 856 1 (term_rg2 + term_rg3)*pdelrho(n,1,1,ipert)*ptmp(1) 857 2 + (term_gg2 + term_gg3)*pdelrho(n,1,1,ipert)*ptmp(2) 858 3 + term_pdelrho*delrho(n,1,1) 859 pdelrho(n,2,1,ipert) = 860 1 (term_rg2 + term_rg3)*pdelrho(n,2,1,ipert)*ptmp(1) 861 2 + (term_gg2 + term_gg3)*pdelrho(n,2,1,ipert)*ptmp(2) 862 2 + term_pdelrho*delrho(n,2,1) 863 pdelrho(n,3,1,ipert) = 864 1 (term_rg2 + term_rg3)*pdelrho(n,3,1,ipert)*ptmp(1) 865 2 + (term_gg2 + term_gg3)*pdelrho(n,3,1,ipert)*ptmp(2) 866 3 + term_pdelrho*delrho(n,3,1) 867 enddo 868 endif 869 endif 870 enddo 871c 872c Put delrho back the way it was since it may be used later on 873c Daniel (3-6-13): It seems strange that we divide the perturbed density 874c and perturbed density gradient by 2 above, but shouldn't multiply them 875c by 2 here. I guess we are somehow building perturbed densities that 876c are 2x what they should be. The code as written below gives the 877c correct answer. Note that this is identical to what is done in 878c xc_cpks_coeff. 879 if (ipol.eq.1 .and. grad) then 880 call dscal(nq*3*ipol,2d0,delrho,1) 881 endif 882c TEST 883c write(6,*) 'Amat3(n,D3_RA_RA_RA)', Amat3(n,D3_RA_RA_RA) 884c write(6,*) 'Amat3(n,D3_RA_RA_RB)', Amat3(n,D3_RA_RA_RB) 885c write(6,*) 'Amat3(n,D3_RA_RB_RB)', Amat3(n,D3_RA_RB_RB) 886c write(6,*) 'Amat3(n,D3_RB_RB_RB)', Amat3(n,D3_RB_RB_RB) 887c TEST 888c ------ 889c Return 890c ------ 891 return 892 end 893