1c 2C$Id$ 3c 4c Modified to handle CPKS while reusing existing code 5c 6c BGJ - 8/98 7c 8C> \ingroup nwdft_xc 9C> @{ 10C> 11C> \brief Evaluate various matrix elements over the functionals 12C> 13C> This routine evaluates a variety of matrix elements that involve the 14C> density functional. 15C> 16C> Note that there are some inconsistencies with respect to the kinetic 17C> energy density. NWChem defines the kinetic energy density as 18C> \f{eqnarray*}{ 19C> \tau_1 &=& \sum_i \frac{1}{2}\langle\phi_i|\nabla^2|\phi_i\rangle 20C> \f} 21C> In most density functionals the kinetic energy density is defined as 22C> \f{eqnarray*}{ 23C> \tau_2 &=& \sum_i \langle\phi_i|\nabla^2|\phi_i\rangle 24C> \f} 25C> the reason being that the factor 1/2 can easily be absorpted into the 26C> density functional parameters. As a result the density functional routines 27C> evaluate 28C> \f{eqnarray*}{ 29C> e &=& f(\rho,\gamma,\tau_2) \\\\ 30C> Amat = &=& \frac{\partial f(\rho,\gamma,\tau_2)}{\partial\rho} \\\\ 31C> Cmat = &=& \frac{\partial f(\rho,\gamma,\tau_2)}{\partial\gamma} \\\\ 32C> Mmat = &=& \frac{\partial f(\rho,\gamma,\tau_2)}{\partial\tau_2} \\\\ 33C> \f} 34C> The NWDFT module in NWChem handles this situation by conceptually wrapping 35C> the functional routines used (the actual implementation just scales 36C> \f$\tau\f$ by a factor 2) 37C> \f{eqnarray*}{ 38C> e &=& g(\rho,\gamma,\tau_1) \\\\ 39C> g(\rho,\gamma,\tau_1) &=& f(\rho,\gamma,2\tau_1) \\\\ 40C> \f} 41C> As a result, however, the density functional routines now return 42C> \f{eqnarray*}{ 43C> e &=& g(\rho,\gamma,\tau_1) \\\\ 44C> Amat = &=& \frac{\partial g(\rho,\gamma,\tau_1)}{\partial\rho} \\\\ 45C> Cmat = &=& \frac{\partial g(\rho,\gamma,\tau_1)}{\partial\gamma} \\\\ 46C> Mmat = &=& \frac{\partial g(\rho,\gamma,\tau_1)}{\partial\tau_2} \\\\ 47C> \f} 48C> where there is an inconsistency in the derivative wrt. to \f$\tau\f$ as 49C> \f$g\f$ is a function of \f$\tau_1\f$ but differentiated wrt. \f$\tau_2\f$. 50C> 51C> In the NWXC module the differentation techniques deployed do not allow 52C> such inconsistencies and \f$g\f$ is strictly differentiated wrt. 53C> \f$\tau_1\f$. This results in a factor \f$2\f$ difference that we need to 54C> compensate for. 55C> 56 Subroutine xc_tabcd(what,l3d_dum, 57 , tol_rho, Tmat, TTmat, Amat, Bmat, Cmat, Dmat, 58 N Emat,Fmat,qxyz,xyz, 59 & chi, delchi, heschi, 60 N curatoms,ncuratoms,nmat, 61 I ipol, nq, nbf, max_at_bf, max_at_bf2, 62 G GRAD, basis, natoms, iniz, ifin, 63 & g_vxc, ibf, rchi_atom, rdelchi_atom, 64 & rdens_atom, cetobfr,kske,Mmat,kslap,Gmat,Lmat,scr) 65c 66c We are using xc_tabcd for CPKS purposes 67c 68 implicit none 69#include "mafdecls.fh" 70#include "global.fh" 71#include "dftpara.fh" 72#include "dft2drv.fh" 73#include "dist.fh" 74#include "dft_fdist.fh" 75#include "util.fh" 76c 77 Logical GRAD !< [Input] .True. if functional depends on density 78 !< gradient 79 integer what !< [Input] What should be calculated: 80 !< - what=0: Do everything (?) 81 !< - what=1: CPKS_LHS 82 !< - what=2: CPKS_RHS 83 !< - what=3: NMR_RHS 84 integer basis !< [Input] The basis set handle 85 integer max_at_bf !< [Input] The maximum number of basis functions 86 !< on an atom 87 integer max_at_bf2 !< [Input] The maximum number of basis 88 !< functions on an atom 89 integer nmat !< [Input] Number of XC matrices (alpha + beta sets) 90 !< to compute, e.g. the number of perturbations for 91 !< CPSCF. 92 integer ipol !< [Input] The number of spin channels 93 integer nq !< [Input] The number of grid points 94 integer nbf !< [Input] The number of basis functions 95 integer natoms !< [Input] The number of atoms 96 integer ncuratoms !< [Input] The number of current "active" atoms 97 integer curatoms(*) !< [Input] Mapping array for current atoms 98 integer jatcur, nu, nu1, indT 99 double precision tol_rho !< [Input] The electron density threshold 100 integer g_vxc(*) !< [Input] The GA handle for Fock matrix 101 !< derivatives 102 double precision Tmat(*) !< [Scratch] 103 double precision TTmat(*) !< [Scratch] 104 double precision rchi_atom(natoms) !< [Input] The maximum basis 105 !< function radius for each 106 !< atom in this call 107 double precision rdelchi_atom(natoms) !< [Input] The maximum 108 !< radius of the basis function gradient for each atom in 109 !< this call 110 double precision rdens_atom(natoms,natoms,ipol) !< Not used 111 integer cetobfr(2,natoms) !< [Input] Mapping from center to 112 !< basis functions (how different from `iniz` and `ifin`?): 113 !< - cetobfr(1,*): first basis function 114 !< - cetobfr(2,*): last basis function 115 logical l3d_dum !< [Junk] 116c 117c Sampling Matrices for the XC Potential & Energy 118c 119 double precision Amat(nq,ipol,*) !< [Input] The derivative wrt rho 120 double precision Cmat(nq,3,ipol,*) !< [Input] The derivative wrt 121 !< rgamma 122 double precision Mmat(nq,ipol,*) !< [Input] The derivative wrt tau 123 double precision Lmat(nq,ipol,*) !< [Input] The derivative wrt lap 124 double precision Gmat(nq,max_at_bf) !< [Scratch] 125c 126c nmr 127 double precision Emat(nq,max_at_bf) !< [Scratch] NMR only 128 double precision Fmat(nq,3,max_at_bf) !< [Scratch] NMR only 129 double precision qxyz(3,*) !< [Input] Grid point coordinates 130 double precision xyz(3,*) !< [Input] Nuclear coordinates 131 logical kske !< [Input] .True. if functional depends on kinetic 132 !< energy density 133 logical kslap !< [Input] .True. if functional depends on laplacian 134c#elif defined(TABCD_CPKS_LHS) 135c 136c Note: Meaning of dimensioning of Amat and Cmat changes for 137c second derivatives, simulating "overloading" of 138c Amat and Cmat 139c 140c Sampling Matrices for the XC part of integrand when making 141c multiple matrices, e.g. XC part of perturbations for CPSCF 142c 143c double precision Amat(nq,ipol,nmat), Cmat(nq,3,ipol,nmat) 144c#elif defined(TABCD_CPKS_RHS) 145c 146c For explicit nuclear derivatives of XC matrix, the same functional 147c derivative values are combined with different basis fn derivatives 148c 149c double precision Amat(nq,ipol), Cmat(nq,3,ipol) 150 151c 152c Sampling Matrices for [Products of] Basis Functions & Gradients 153c 154 integer imat ! XC matrix loop index 155 double precision Bmat(nq,max_at_bf) !< [Scratch] 156 double precision Dmat(nq,3,max_at_bf) !< [Scratch] 157 integer iniz(natoms) !< [Input] The first basis function for each 158 !< atom 159 integer ifin(natoms) !< [Input] The last basis function for each 160 !< atom 161c 162c Basis Functions & Gradients 163c 164 double precision chi(nq,nbf) !< [Input] The value of the basis 165 !< functions at the grid points 166 double precision delchi(nq,3,nbf) !< [Input] The value of the 167 !< gradient of the basis 168 !< functions at the grid points 169 double precision heschi(nq,6,*) !< [Input] The value of the 170 !< Hessian of the basis 171 !< functions at the grid points 172 integer ibf(nbf) !< [Input] The rank of the basis function for 173 !< every basis function (why do we need this?) 174 double precision A_MAX, C_MAX, AC_MAX, FUNC_MAXI, 175 & B_MAX, D_MAX, BD_MAX, FUNC_MAXJ 176 integer iat, inizia, ifinia, nbfia, nnia, ifirst, ilast 177 integer jat, inizja, ifinja, nbfja, nnja, jfirst, jlast 178 integer ii, mu, mu1 179 integer n,lastjat 180 double precision chi1 181 double precision scr(nq) !< [Scratch] Temporary stuff in 182 !< quadrature 183 double precision dabsmax 184 double precision tolrho15 185 external dabsmax 186 logical l3d !< [Input] .True. if XC-matrices stored in a 3D GA 187 integer jrsh,jrsh2,n3d,idir,jdir 188 integer g_update(2) 189cnmr 190 double precision Rij(3) ! vector R(jat) - R(iat) 191 integer inia, iq, ix, ix1, ix2 192c 193c Indexing array for basis function hessian columns as if 194c it were a 3x3 matrix 195c 196 integer indh(3,3) 197 logical w01,w02,w013,dofull 198 double precision raa,rbb 199 integer nbhandl1,nbhandl2 200 logical nbfirst1,nbfirst2,doitt 201 integer sizeblk 202#include "nwc_const.fh" 203 integer nonzero(nw_max_atom),natleft, 204 A iat0,jat0 205 data indh / 1, 2, 3, 206 & 2, 4, 5, 207 & 3, 5, 6 / 208 data nbhandl1 /0./ 209 data nbhandl2 /0./ 210 save nbhandl1 211 save nbhandl2 212c 213c 0: l3d=.f. & n3d=1 214ccc rhs: l3d=.true. & n3d=3 215ccc lhs: l3d=.true. & n3d=1 216c 217 call starttimer(monitor_tabcd) 218 call starttimer(monitor_screen0) 219c lingering nbacc from previous calls 220 if(nbhandl1.ne.0) call ga_nbwait(nbhandl1) 221 if(nbhandl2.ne.0) call ga_nbwait(nbhandl2) 222 l3d = ga_ndim(g_vxc).eq.3 223 natleft=0 224 do iat = 1, natoms 225 if (iniz(iat).ne.0) then 226 natleft=natleft+1 227 nonzero(natleft)=iat 228 endif 229 enddo 230 tolrho15=tol_rho**1.25d0 231 if (what.eq.0) then 232 n3d=1 233 elseif (what.eq.1) then 234 n3d=1 235 elseif (what.eq.2) then 236 n3d=3 237 elseif (what.eq.3) then 238 n3d=3 239 else 240 call errquit(' wrong what value for xctabcd ',0,0) 241 endif 242 w01=what.eq.0.or.what.eq.1 243 w013=w01.or.what.eq.3 244 w02=what.eq.0.or.what.eq.2 245 dofull=what.ne.0.or.l3d 246 nbfirst1=.true. 247 nbfirst2=.true. 248 call endtimer(monitor_screen0) 249c 250c Beginning of loop over multiple XC matrices 251c 252 do 500 imat = 1,nmat 253c 254c Compute the matrix product for the XC potential and energy: 255c 256c T = transpose(A*B) + transpose(C*D) 257c 258 259 call starttimer(monitor_screen1) 260 261 A_MAX = dabsmax(nq*ipol,Amat(1,1,imat)) 262 if (GRAD) then 263 C_MAX = dabsmax(nq*3*ipol,Cmat(1,1,1,imat)) 264 else 265 C_MAX = 0d0 266 endif 267 AC_MAX = max(A_MAX,C_MAX) 268 269 call endtimer(monitor_screen1) 270 271c 272c repl stuff 273c 274 if (xcreplicated.and.dorepxc) then 275 g_update(1)=k_repxc(1) 276 g_update(2)=k_repxc(2) 277 else 278 g_update(1)=g_vxc(1) 279 g_update(2)=g_vxc(2) 280 endif 281 do 430 iat0=1,natleft 282 call starttimer(monitor_screen2) 283 iat=nonzero(iat0) 284 inizia = iniz(iat) 285 FUNC_MAXI = rchi_atom(iat) 286 if(GRAD) FUNC_MAXI = max(FUNC_MAXI,rdelchi_atom(iat)) 287 doitt=(AC_MAX*FUNC_MAXI).ge.tol_rho 288 call endtimer(monitor_screen2) 289 if (what.eq.2.or.doitt) then 290 ifinia = ifin(iat) 291 ifirst = cetobfr(1,iat) 292 ilast = cetobfr(2,iat) 293 nnia = ifinia - inizia + 1 294 nbfia = ilast - ifirst + 1 295 do ii = 1, ipol 296 do mu = 1, nnia 297 298 call starttimer(monitor_mult1) 299 300 mu1 = mu+inizia-1 301 if (GRAD) then 302 do n = 1, nq 303 chi1 = chi(n,mu1) 304 305 306 Bmat(n,mu) = Amat(n,ii,imat)*chi1 + 307 & delchi(n,1,mu1)*Cmat(n,1,ii,imat) + 308 & delchi(n,2,mu1)*Cmat(n,2,ii,imat) + 309 & delchi(n,3,mu1)*Cmat(n,3,ii,imat) 310 Dmat(n,1,mu) = Cmat(n,1,ii,imat)*chi1 311 Dmat(n,2,mu) = Cmat(n,2,ii,imat)*chi1 312 Dmat(n,3,mu) = Cmat(n,3,ii,imat)*chi1 313 314 if(kske) then 315 Dmat(n,1,mu) = Dmat(n,1,mu) + 316 & Mmat(n,ii,imat)*delchi(n,1,mu1) 317 Dmat(n,2,mu) = Dmat(n,2,mu) + 318 & Mmat(n,ii,imat)*delchi(n,2,mu1) 319 Dmat(n,3,mu) = Dmat(n,3,mu) + 320 & Mmat(n,ii,imat)*delchi(n,3,mu1) 321 endif 322 323 if(kslap) then 324 Dmat(n,1,mu) = Dmat(n,1,mu) + 325 & 2d0*Lmat(n,ii,imat)*delchi(n,1,mu1) 326 Dmat(n,2,mu) = Dmat(n,2,mu) + 327 & 2d0*Lmat(n,ii,imat)*delchi(n,2,mu1) 328 Dmat(n,3,mu) = Dmat(n,3,mu) + 329 & 2d0*Lmat(n,ii,imat)*delchi(n,3,mu1) 330 Bmat(n,mu) = Bmat(n,mu) + Lmat(n,ii,imat)* 331 & (heschi(n,1,mu1)+heschi(n,4,mu1)+heschi(n,6,mu1)) 332 Gmat(n,mu) = Lmat(n,ii,imat)*chi1 333 endif 334 335 336 enddo 337 else 338 do n = 1, nq 339 Bmat(n,mu) = chi(n,mu1)*Amat(n,ii,imat) 340 enddo 341 endif ! GRAD 342 enddo ! mu 343c Monitoring 344 345 call endtimer(monitor_mult1) 346 347c 348 call starttimer(monitor_screen3) 349 B_MAX = dabsmax(nnia*nq,Bmat) 350 if (GRAD) then 351 D_MAX = dabsmax(nnia*nq*3,Dmat) 352 else 353 D_MAX = 0d0 354 endif 355 BD_MAX = max(B_MAX,D_MAX) 356c 357 lastjat=iat0 358 if (what.eq.2) lastjat=natleft 359 if (what.eq.3) lastjat=iat0-1 360 call endtimer(monitor_screen3) 361 do 168 jat0=1,lastjat 362 jat=nonzero(jat0) 363 if(what.eq.2) then 364c 365c To fit better into existing structure, loop over full square 366c of atom pairs and only compute nuclear derivative contribution 367c from jat. Also, this way we only need check jatcur once and 368c for all, and do not have to check iatcur at all. 369c 370 jatcur = curatoms(jat) 371 if (jatcur.eq.0) goto 168 372 endif 373 call starttimer(monitor_screen4) 374 inizja = iniz(jat) 375 FUNC_MAXJ = rchi_atom(jat) 376 if(grad) FUNC_MAXJ = max(rdelchi_atom(jat),FUNC_MAXJ) 377 doitt=(BD_MAX*FUNC_MAXJ).ge.tol_rho 378 call endtimer(monitor_screen4) 379 if (what.eq.2.or.doitt) then 380 381c Monitoring 382 383 call starttimer(monitor_mult2) 384 385 if (what.eq.3) then 386 Rij(1) = 0.5d0*(xyz(1,jat)-xyz(1,iat)) 387 Rij(2) = 0.5d0*(xyz(2,jat)-xyz(2,iat)) 388 Rij(3) = 0.5d0*(xyz(3,jat)-xyz(3,iat)) 389 endif 390 ifinja = ifin(jat) 391 jfirst = cetobfr(1,jat) 392 jlast = cetobfr(2,jat) 393 nbfja = jlast - jfirst + 1 394 nnja = ifinja - inizja + 1 395 sizeblk=n3d*nbfia*nbfja 396 if (what.eq.2.or.what.eq.3) 397 Y call dcopy(sizeblk, 0d0,0, TTmat,1) 398c 399c Loop over x, y, z directions for derivative XC mats 400c 401 do jdir = 1,n3d 402c 403 if (what.eq.3) then 404 ix1 = mod(jdir,3)+1 405 ix2 = mod(jdir+1,3)+1 406 raa=rij(ix1) 407 rbb=rij(ix2) 408 do iq = 1, nq 409 scr(iq) = raa*qxyz(ix2,iq) - rbb*qxyz(ix1,iq) 410 enddo 411 do inia = 1, nnia 412 do iq = 1, nq 413 Emat(iq,inia) = scr(iq)*Bmat(iq,inia) 414 enddo 415 if (GRAD) then 416 do iq = 1, nq 417 Emat(iq,inia) = Emat(iq,inia)+ 418 & (raa*Dmat(iq,ix2,inia) 419 & - rbb*Dmat(iq,ix1,inia)) 420 enddo 421 endif 422 enddo 423 if (GRAD) then 424 do inia = 1, nnia 425 do iq = 1, nq 426 Fmat(iq,1,inia) = scr(iq) 427 & * Dmat(iq,1,inia) 428 Fmat(iq,2,inia) = scr(iq) 429 & * Dmat(iq,2,inia) 430 Fmat(iq,3,inia) = scr(iq) 431 & * Dmat(iq,3,inia) 432 enddo 433 enddo 434 endif 435 endif ! what.eq.3 436c 437c Daniel (2-7-13): Here, we build matrix elements of the XC-kernel when 438c what = 1 for CPKS RHS with fxc*drhonuc*chi_mu in Bmat. 439 if (w01) then 440 call dgemm('T', 'N', nnia, nnja, nq, 1.d0, Bmat, 441 & nq, chi(1,inizja), nq, 0.d0, Tmat, nnia) 442 elseif (what.eq.3) then 443 call dgemm('T', 'N', nnia, nnja, nq, 1.0d0, Emat, 444 & nq, chi(1,inizja), nq, 0.0d0, Tmat, nnia) 445 else 446c Note the sign change for a nuclear derivative, and also that the 447c leading dimension of delchi must be set correctly 448 call dgemm('T', 'N', nnia, nnja, nq, -1.d0, Bmat, 449 & nq, delchi(1,jdir,inizja), nq*3, 0.d0, Tmat, 450 & nnia) 451 endif 452 if (GRAD) then 453 if (w01) then 454 call dgemm('T', 'N', nnia, nnja, 3*nq, 455 & 1.d0, Dmat, 3*nq, delchi(1,1,inizja), 456 & 3*nq, 1.d0, Tmat, nnia) 457 elseif (what.eq.3) then 458 call dgemm('T', 'N', nnia, nnja, 3*nq, 459 & 1.0d0, Fmat, 3*nq, delchi(1,1,inizja), 460 & 3*nq, 1.0d0, Tmat, nnia) 461 else 462 indT = 0 463 do nu = 1, nnja 464 nu1 = nu+inizja-1 465 do mu = 1, nnia 466 indT = indT + 1 467 Tmat(indT) = Tmat(indT)- 468 * ddot(nq,Dmat(1,1,mu),1, 469 & heschi(1,indh(1,jdir),nu1),1) - 470 * ddot(nq,Dmat(1,2,mu),1, 471 & heschi(1,indh(2,jdir),nu1),1) - 472 * ddot(nq,Dmat(1,3,mu),1, 473 & heschi(1,indh(3,jdir),nu1),1) 474 enddo 475 enddo 476 endif 477 endif ! GRAD 478 if(kslap) then 479 if (w01) then 480 indT=0 481 do nu = 1, nnja 482 nu1 = nu+inizja-1 483 do mu = 1, nnia 484 indT = indT + 1 485 Tmat(indT) = Tmat(indT) + 486 & ddot(nq,Gmat(1,mu),1,heschi(1,1,nu1),1) + 487 & ddot(nq,Gmat(1,mu),1,heschi(1,4,nu1),1) + 488 & ddot(nq,Gmat(1,mu),1,heschi(1,6,nu1),1) 489 enddo 490 enddo 491 endif 492 endif 493c Daniel (2-7-13): For the CPKS RHS stuff, the first call has what=1, 494c so n3d=1. Also, dofull is true in that case. 495 if (n3d.eq.1) then 496 call dfill(max_at_bf2, 0.d0, TTmat, 1) 497 if (dofull) then 498 call scat_mat(TTmat, Tmat, nbfia, nbfja, nnia, 499 & nnja,ifirst,jfirst,ibf(inizia),ibf(inizja)) 500 else 501 call scat_matup(TTmat, Tmat, nbfia, nbfja, nnia, 502 & nnja,ifirst,jfirst,ibf(inizia),ibf(inizja)) 503 endif 504c Daniel (2-7-13): For the CPKS RHS stuff, the second call has what=2, 505c so n3d=3. Also, dofull is true in that case. 506 else 507 call scat_mat3(n3d,jdir, 508 & TTmat, Tmat, nbfia, nbfja, nnia, 509 & nnja,ifirst,jfirst,ibf(inizia),ibf(inizja)) 510 endif 511 512 enddo ! jdir (loop over x, y, z directions for nmr) 513c Monitoring 514 515 call endtimer(monitor_mult2) 516 517 doitt=.true. 518 call starttimer(monitor_screen5) 519 if (what.eq.0) then 520 doitt=dabsmax(sizeblk,ttmat).ge.tolrho15 521 jrsh=ii 522 elseif (what.eq.1) then 523 doitt=.true. 524 jrsh=imat+(ii-1)*nmat 525 elseif (what.eq.3) then 526 jrsh=(ii-1)*n3d+1 527 else 528 doitt=dabsmax(sizeblk,ttmat).ge.tol_rho 529 jrsh=1+(jat-1)*3+(ii-1)*3*natoms 530 endif 531 call endtimer(monitor_screen5) 532 if (doitt) then 533 534 jrsh2=jrsh+n3d-1 535 536c Monitoring 537 538 call updist(monitor_size_ga_acc1, sizeblk) 539 call starttimer( monitor_comm_ga_acc1) 540 541c Daniel (2-7-13): l3d is always true for the CPKS RHS stuff. 542 if (l3d) then 543 call dft_3dacc(g_vxc, ttmat, 544 & jrsh,jrsh2, 545 % ifirst, ilast, jfirst, jlast, nbfia) 546 else 547 if (dftnbacc) then 548 if (.not.nbfirst1) then 549 call starttimer( monitor_wait1) 550 call ga_nbwait(nbhandl1) 551 call endtimer( monitor_wait1) 552 endif 553 nbfirst1=.false. 554 call upd_atombl_nb(g_update(ii), 555 . basis,iat,jat,ttmat,nbhandl1) 556 else 557 if (truerepxc) then 558 call xc_atom_blockd(dbl_mb(k_repxc(ii)), 559 N nbf_ld,basis,iat,jat,ttmat) 560 else 561 if (truerepxc) then 562 call xc_atom_block(dbl_mb(k_repxc(ii)), 563 N nbf_ld,basis,jat,iat,tmat) 564 else 565 call upd_atom_block(g_update(ii), 566 . basis,iat,jat,ttmat) 567 endif 568 endif 569 endif 570 endif ! l3d 571c Monitoring 572 573 call endtimer( monitor_comm_ga_acc1) 574 575c Daniel (2-7-13): l3d is always true for the CPKS RHS stuff. 576 if (what.ne.0.or.l3d) then 577c 578c check to see if can skip and use ga_symmetrize 579c 580c Daniel (2-7-13): This part is always done for the second call to 581c xc_tabcd in CPKS RHS stuff. For the first call, it only happens 582c for off-diagonal terms. 583 if ((w013.and.iat.ne.jat).or.what.eq.2) then 584c For CPKS RHS, we update with transpose even for iat = jat, 585c since that is necessary to get both contributions 586c mu * del(nu) and del(mu) * nu 587 588 589 call starttimer(monitor_comp_transp) 590 591c Daniel (2-7-13): This happens for the first call for CPKS RHS stuff 592c (n3d=1). 593 if (n3d.eq.1) then 594 call transp_mat(TTmat, Tmat, 595 , nbfia, nbfja) 596c Daniel (2-7-13): This happens for the second call for CPKS RHS stuff 597c (n3d=3). 598 else 599 if (what.eq.3) then 600 call dscal(n3d*nbfia*nbfja,-1.0d0,TTmat,1) 601 endif 602 call transp_mat3(n3d,TTmat, Tmat, 603 , nbfia, nbfja) 604 endif 605 606 call endtimer(monitor_comp_transp) 607 608 609c Monitoring 610 611 call starttimer(monitor_comm_ga_acc2) 612 613 614c Daniel (2-7-13): For CPKS RHS stuff, l3d=true always. 615 if (l3d) then 616 call dft_3dacc(g_vxc, tmat, 617 & jrsh,jrsh2, 618 % jfirst, jlast, ifirst, ilast, nbfja) 619 else 620 if (dftnbacc) then 621 if (.not.nbfirst2) then 622 call ga_nbwait(nbhandl2) 623 endif 624 nbfirst2 = .false. 625 call upd_atombl_nb(g_update(ii), 626 . basis,jat,iat,tmat,nbhandl2) 627 else 628 call upd_atom_block(g_update(ii),basis, 629 J jat,iat,tmat) 630 endif 631 endif 632c Monitoring 633 634 call endtimer(monitor_comm_ga_acc2) 635 636 637 endif ! (w013.and.iat.ne.jat).or.what.eq.2 638 endif ! what.ne.3.or.l3d 639 endif ! doitt 640 endif ! what.eq.2.or.doitt 641 168 continue ! jat0 loop 642 enddo ! ipol loop 643 endif ! what.eq.2.or.doitt 644 430 continue ! iat0 loop 645 500 continue ! imat loop 646 call endtimer(monitor_tabcd) 647 648 return 649 end 650 subroutine upd_atombl_nb(g_array, basis, iat, jat, buf, 651 $ nbhandle) 652 implicit none 653#include "errquit.fh" 654#include "global.fh" 655#include "bas.fh" 656c 657 integer g_array, basis, iat, jat 658 integer nbhandle 659 double precision buf(*) 660 logical status 661c 662 integer ilo, ihi, jlo, jhi, idim, jdim 663c 664c add atom block buf info of the matrix g_array (over basis functions) 665c 666 status= bas_ce2bfr(basis, iat, ilo, ihi) 667 status=status.and. bas_ce2bfr(basis, jat, jlo, jhi) 668 if (.not. status) 669 $ call errquit('upd_atom_block: ce2bfr failed', 0, BASIS_ERR) 670c 671 idim = ihi - ilo + 1 672 jdim = jhi - jlo + 1 673c 674c clearing lingering nbacc 675c 676c if(nbhandle.ne.0) call ga_nbwait(nbhandle) 677 nbhandle=0 678 if (idim.gt.0 .and. jdim.gt.0) 679 $ call ga_nbacc(g_array, ilo, ihi, jlo, jhi, buf, idim, 680 $ 1.0d0,nbhandle) 681c 682 end 683C> 684C> @} 685