1czora...Build the zora contribution to the fock matrix 2czora...This is routine is analogous to the xc_tabcd routine that 3czora...contructs the exchange-correlation contribution 4czora...TODO: Clean up and more comments 5c 6 subroutine zora_tabcd(what,l3d, 7 , tol_rho, Tmat, TTmat, Amat, Bmat, Cmat, Dmat, 8 N Emat,Fmat,qxyz,xyz, 9 & chi, delchi, heschi, 10 N curatoms,ncuratoms,nmat, 11 I ipol, nq, nbf, max_at_bf, max_at_bf2, 12 G GRAD, basis, natoms, iniz, ifin, 13 & g_zora, ibf, rchi_atom, rdelchi_atom, 14 & rdens_atom, cetobfr,kske,Mmat,scr) 15c 16 implicit none 17#include "mafdecls.fh" 18#include "dftpara.fh" 19#include "dft2drv.fh" 20#include "dist.fh" 21#include "dft_fdist.fh" 22c 23czora...start 24#include "zora.fh" 25czora...end 26c 27c 28 Logical GRAD 29 integer what ! [in] 30cwhat=0 31cwhat=1 CPKS_LHS 32cwhat=2 CPKS_RHS 33cwhat=2 NMR_RHS 34 integer basis 35 integer max_at_bf ! [input] 36 integer max_at_bf2 ! [input] 37 integer nmat ! Number of XC matrices (alpha + beta sets) to make [input] 38c ! e.g. number of perturbations for CPSCF 39 integer imat ! XC matrix loop index 40 integer ipol ! [input] 41 integer nq ! [input] 42 integer nbf ! [input] 43 integer natoms ! [input] 44 integer ncuratoms ! Number of current "active" atoms [input] 45 integer curatoms(*) ! Mapping array for current atoms [input] 46 integer jatcur, nu, nu1, indT 47 double precision tol_rho 48 double precision Tmat(*), TTmat(*) 49 double precision rchi_atom(natoms) 50 double precision rdelchi_atom(natoms) 51 double precision rdens_atom(natoms,natoms,ipol) 52 integer cetobfr(2,natoms) 53c 54czora...start 55 integer g_zora(*) 56czora...end 57c 58c Sampling Matrices for the XC Potential & Energy 59c 60 double precision Amat(nq,ipol,*), Cmat(nq,3,ipol,*) 61 double precision Mmat(nq,ipol) 62c nmr 63 double precision Emat(nq,max_at_bf) 64 double precision Fmat(nq,3,max_at_bf) 65 double precision qxyz(3,*) ! [input] grid point coordinates 66 double precision xyz(3,*) ! [input] nuclear coordinates 67 logical kske 68c#elif defined(TABCD_CPKS_LHS) 69c 70c Note: Meaning of dimensioning of Amat and Cmat changes for 71c second derivatives, simulating "overloading" of 72c Amat and Cmat 73c 74c Sampling Matrices for the XC part of integrand when making 75c multiple matrices, e.g. XC part of perturbations for CPSCF 76c 77c double precision Amat(nq,ipol,nmat), Cmat(nq,3,ipol,nmat) 78c#elif defined(TABCD_CPKS_RHS) 79c 80c For explicit nuclear derivatives of XC matrix, the same functional 81c derivative values are combined with different basis fn derivatives 82c 83c double precision Amat(nq,ipol), Cmat(nq,3,ipol) 84 85c 86c Sampling Matrices for [Products of] Basis Functions & Gradients 87c 88 double precision Bmat(nq,max_at_bf) 89 double precision Dmat(nq,3,max_at_bf) 90 integer iniz(natoms), ifin(natoms) 91c 92c Basis Functions & Gradients 93c 94 double precision chi(nq,nbf), delchi(nq,3,nbf) 95 double precision heschi(nq,6,*) 96 integer ibf(nbf) 97 double precision A_MAX, C_MAX, AC_MAX, FUNC_MAXI, 98 & B_MAX, D_MAX, BD_MAX, FUNC_MAXJ 99 integer iat, inizia, ifinia, nbfia, nnia, ifirst, ilast 100 integer jat, inizja, ifinja, nbfja, nnja, jfirst, jlast 101 integer ii, mu, mu1 102 integer n,lastjat 103 double precision chi1 104 double precision scr(nq) 105 double precision dabsmax 106 double precision tolrho15 107 external dabsmax 108 logical l3d 109 integer jrsh,jrsh2,n3d,idir,jdir 110 integer g_update(2) 111cnmr 112 double precision Rij(3) ! vector R(jat) - R(iat) 113 integer inia, iq, ix, ix1, ix2 114c 115c Indexing array for basis function hessian columns as if 116c it were a 3x3 matrix 117c 118 integer indh(3,3) 119 logical w01,w02,w013,dofull 120 double precision ddot 121 double precision raa,rbb 122 integer nbhandl1,nbhandl2 123 logical nbfirst1,nbfirst2,doitt 124 integer sizeblk 125#include "nwc_const.fh" 126 integer nonzero(nw_max_atom),natleft, 127 A iat0,jat0 128 external ddot 129 data indh / 1, 2, 3, 130 & 2, 4, 5, 131 & 3, 5, 6 / 132 133c 0: l3d=.f. & n3d=1 134ccc rhs: l3d=.true. & n3d=3 135ccc lhs: l3d=.true. & n3d=1 136c 137 call starttimer(monitor_tabcd) 138 call starttimer(monitor_screen0) 139 natleft=0 140 do iat = 1, natoms 141 if(iniz(iat).ne.0) then 142 natleft=natleft+1 143 nonzero(natleft)=iat 144 endif 145 enddo 146 tolrho15=tol_rho**1.25d0 147 if(what.eq.0) then 148 n3d=1 149 elseif(what.eq.1) then 150 n3d=1 151 elseif(what.eq.2) then 152 n3d=3 153 elseif(what.eq.3) then 154 n3d=3 155 else 156 call errquit(' wrong what value for xctabcd ',0,0) 157 endif 158 w01=what.eq.0.or.what.eq.1 159 w013=w01.or.what.eq.3 160 w02=what.eq.0.or.what.eq.2 161 dofull=what.ne.0.or.l3d 162 nbfirst1=.true. 163 nbfirst2=.true. 164 call endtimer(monitor_screen0) 165c 166c 167c Beginning of loop over multiple XC matrices 168c 169 do 500 imat = 1,nmat 170c 171c Compute the matrix product for the XC potential and energy: 172c 173c T = transpose(A*B) + transpose(C*D) 174c 175 176 call starttimer(monitor_screen1) 177 178 A_MAX = dabsmax(nq*ipol,Amat(1,1,imat)) 179 if (GRAD) then 180 C_MAX = dabsmax(nq*3*ipol,Cmat(1,1,1,imat)) 181 else 182 C_MAX = 0d0 183 endif 184 AC_MAX = max(A_MAX,C_MAX) 185 186 call endtimer(monitor_screen1) 187 188c 189c repl stuff 190c 191 if(xcreplicated.and.dorepxc) then 192 g_update(1)=k_repxc(1) 193 g_update(2)=k_repxc(2) 194 else 195 g_update(1)=g_zora(1) 196 g_update(2)=g_zora(2) 197 endif 198 do 430 iat0=1,natleft 199 call starttimer(monitor_screen2) 200 iat=nonzero(iat0) 201 inizia = iniz(iat) 202 FUNC_MAXI = rchi_atom(iat) 203 if(GRAD) FUNC_MAXI = max(FUNC_MAXI,rdelchi_atom(iat)) 204 doitt=(AC_MAX*FUNC_MAXI).ge.tol_rho 205 call endtimer(monitor_screen2) 206 if(what.eq.2.or.doitt) then 207 ifinia = ifin(iat) 208 ifirst = cetobfr(1,iat) 209 ilast = cetobfr(2,iat) 210 nnia = ifinia - inizia + 1 211 nbfia = ilast - ifirst + 1 212 do ii = 1, ipol 213 do mu = 1, nnia 214 215 call starttimer(monitor_mult1) 216 217 mu1 = mu+inizia-1 218 if (GRAD) then 219 do n = 1, nq 220 chi1 = chi(n,mu1) 221 222 Bmat(n,mu) = 0.d0 223 224 if (so_term .eq. 1) then ! z-component 225 Dmat(n,1,mu) = delchi(n,1,mu1)*Amat(n,ii,imat) 226 Dmat(n,2,mu) = delchi(n,2,mu1)*Amat(n,ii,imat) 227 Dmat(n,3,mu) = 0.d0 228 else if (so_term .eq. 2) then ! y-component 229 Dmat(n,3,mu) = delchi(n,3,mu1)*Amat(n,ii,imat) 230 Dmat(n,1,mu) = delchi(n,1,mu1)*Amat(n,ii,imat) 231 Dmat(n,2,mu) = 0.d0 232 else if (so_term .eq. 3) then ! x-component 233 Dmat(n,2,mu) = delchi(n,2,mu1)*Amat(n,ii,imat) 234 Dmat(n,3,mu) = delchi(n,3,mu1)*Amat(n,ii,imat) 235 Dmat(n,1,mu) = 0.d0 236 else if (so_term .eq. 0) then ! spin-free 237 Dmat(n,1,mu) = delchi(n,1,mu1)*Amat(n,ii,imat) 238 Dmat(n,2,mu) = delchi(n,2,mu1)*Amat(n,ii,imat) 239 Dmat(n,3,mu) = delchi(n,3,mu1)*Amat(n,ii,imat) 240 end if 241 242 enddo ! nq: over grid points 243 else 244 do n = 1, nq 245 Bmat(n,mu) = 0.d0 246 enddo ! nq: over grid points 247 endif 248 249 enddo ! mu: over the basis set 250 251c Monitoring 252 253 call endtimer(monitor_mult1) 254c 255 call starttimer(monitor_screen3) 256 B_MAX = dabsmax(nnia*nq,Bmat) 257 if (GRAD) then 258 D_MAX = dabsmax(nnia*nq*3,Dmat) 259 else 260 D_MAX = 0d0 261 endif 262 BD_MAX = max(B_MAX,D_MAX) 263c 264 lastjat=iat0 265 if(what.eq.2) lastjat=natleft 266 if(what.eq.3) lastjat=iat0-1 267 call endtimer(monitor_screen3) 268 do 168 jat0=1,lastjat 269 jat=nonzero(jat0) 270 if(what.eq.2) then 271c 272c To fit better into existing structure, loop over full square 273c of atom pairs and only compute nuclear derivative contribution 274c from jat. Also, this way we only need check jatcur once and 275c for all, and don't have to check iatcur at all. 276c 277 jatcur = curatoms(jat) 278 if (jatcur.eq.0) goto 168 279 endif 280 call starttimer(monitor_screen4) 281 inizja = iniz(jat) 282 FUNC_MAXJ = rchi_atom(jat) 283 if(grad) FUNC_MAXJ = max(rchi_atom(jat),FUNC_MAXJ) 284 doitt=(BD_MAX*FUNC_MAXJ).ge.tol_rho 285 call endtimer(monitor_screen4) 286 if(what.eq.2.or.doitt) then 287 288c Monitoring 289 290 call starttimer(monitor_mult2) 291 292 if(what.eq.3) then 293 Rij(1) = 0.5d0*(xyz(1,jat)-xyz(1,iat)) 294 Rij(2) = 0.5d0*(xyz(2,jat)-xyz(2,iat)) 295 Rij(3) = 0.5d0*(xyz(3,jat)-xyz(3,iat)) 296 endif 297 ifinja = ifin(jat) 298 jfirst = cetobfr(1,jat) 299 jlast = cetobfr(2,jat) 300 nbfja = jlast - jfirst + 1 301 nnja = ifinja - inizja + 1 302 sizeblk=n3d*nbfia*nbfja 303 if(what.eq.2.or.what.eq.3) 304 & call dcopy(sizeblk, 0d0,0, TTmat,1) 305c 306c Loop over x, y, z directions for derivative XC mats 307c 308 do jdir = 1,n3d 309c 310 if(what.eq.3) then 311 ix1 = mod(jdir,3)+1 312 ix2 = mod(jdir+1,3)+1 313 raa=rij(ix1) 314 rbb=rij(ix2) 315 do iq = 1, nq 316 scr(iq) = raa*qxyz(ix2,iq) - rbb*qxyz(ix1,iq) 317 enddo 318 do inia = 1, nnia 319 do iq = 1, nq 320 Emat(iq,inia) = scr(iq)*Bmat(iq,inia) 321 enddo 322 if (GRAD) then 323 do iq = 1, nq 324 Emat(iq,inia) = 0.d0 325 enddo 326 endif 327 enddo 328 if (GRAD) then 329 do inia = 1, nnia 330 do iq = 1, nq 331 Fmat(iq,1,inia) = 0.d0 332 Fmat(iq,2,inia) = 0.d0 333 Fmat(iq,3,inia) = 0.d0 334 enddo 335 enddo 336 endif 337 endif 338c 339 if (GRAD) then 340 indT = 0 341 do nu = 1, nnja 342 nu1 = nu+inizja-1 343 do mu = 1, nnia 344 indT = indT + 1 345 346czora... build Tmat 347 348 if (so_term .eq. 1) then ! z-component 349 Tmat(indT) = 350 & ddot(nq,Dmat(1,1,mu),1,delchi(1,2,nu1),1) 351 & - ddot(nq,Dmat(1,2,mu),1,delchi(1,1,nu1),1) 352 else if (so_term .eq. 2) then ! y-component 353 Tmat(indT) = 354 & ddot(nq,Dmat(1,3,mu),1,delchi(1,1,nu1),1) 355 & - ddot(nq,Dmat(1,1,mu),1,delchi(1,3,nu1),1) 356 else if (so_term .eq. 3) then ! x-component 357 Tmat(indT) = 358 & ddot(nq,Dmat(1,2,mu),1,delchi(1,3,nu1),1) 359 & - ddot(nq,Dmat(1,3,mu),1,delchi(1,2,nu1),1) 360 else if (so_term .eq. 0) then ! spin-free 361 Tmat(indT) = 362 & ddot(nq,Dmat(1,1,mu),1,delchi(1,1,nu1),1) 363 & + ddot(nq,Dmat(1,2,mu),1,delchi(1,2,nu1),1) 364 & + ddot(nq,Dmat(1,3,mu),1,delchi(1,3,nu1),1) 365 end if 366 367 enddo ! mu: over the basis set 368 enddo ! nu: over the basis set 369 370 endif 371 372 if(n3d.eq.1) then 373 call dfill(max_at_bf2, 0.d0, TTmat, 1) 374 if(dofull) then 375 call scat_mat(TTmat, Tmat, nbfia, nbfja, nnia, 376 & nnja,ifirst,jfirst,ibf(inizia),ibf(inizja)) 377 else 378 if (so_term.eq.0) then 379 call scat_matup(TTmat, Tmat, nbfia, nbfja, nnia, 380 & nnja,ifirst,jfirst,ibf(inizia),ibf(inizja)) 381 else 382 call scat_mat(TTmat, Tmat, nbfia, nbfja, nnia, 383 & nnja,ifirst,jfirst,ibf(inizia),ibf(inizja)) 384 end if 385 endif 386 else 387 call scat_mat3(n3d,jdir, 388 & TTmat, Tmat, nbfia, nbfja, nnia, 389 & nnja,ifirst,jfirst,ibf(inizia),ibf(inizja)) 390 endif 391 392 enddo ! jdir (loop over x, y, z directions for nmr) 393 394c Monitoring 395 396 call endtimer(monitor_mult2) 397 398 doitt=.true. 399 call starttimer(monitor_screen5) 400 if(what.eq.0) then 401 doitt=dabsmax(sizeblk,ttmat).ge.tolrho15 402 jrsh=ii 403 elseif(what.eq.1) then 404 doitt=.true. 405 jrsh=imat+(ii-1)*nmat 406 elseif(what.eq.3) then 407 jrsh=(ii-1)*n3d+1 408 else 409 doitt=dabsmax(sizeblk,ttmat).ge.tol_rho 410 jrsh=1+(jat-1)*3+(ii-1)*3*natoms 411 endif 412 call endtimer(monitor_screen5) 413 if(doitt) then 414 415 jrsh2=jrsh+n3d-1 416 417c Monitoring 418 419 call updist(monitor_size_ga_acc1, sizeblk) 420 call starttimer( monitor_comm_ga_acc1) 421 422 if(l3d) then 423 call dft_3dacc(g_zora, ttmat, 424 & jrsh,jrsh2, 425 & ifirst, ilast, jfirst, jlast, nbfia) 426 else 427 if(dftnbacc) then 428 429 if(.not.nbfirst1) then 430 call starttimer( monitor_wait1) 431 call ga_nbwait(nbhandl1) 432 call endtimer( monitor_wait1) 433 endif 434 nbfirst1=.false. 435 call upd_atombl_nb(g_update(ii), 436 & basis,iat,jat,ttmat,nbhandl1) 437 else 438 call upd_atom_block(g_update(ii), 439 & basis,iat,jat,ttmat) 440 441 442 endif 443 endif 444 445 446 447c Monitoring 448 449 call endtimer( monitor_comm_ga_acc1) 450 451 452 if(what.ne.0.or.l3d) then 453c 454c check to see if can skip and use ga_symmetrize 455c 456 if ((w013.and.iat.ne.jat).or.what.eq.2) then 457c For CPKS RHS, we update with transpose even for iat = jat, 458c since that is necessary to get both contributions 459c mu * del(nu) and del(mu) * nu 460 461 462 call starttimer(monitor_comp_transp) 463 464 if(n3d.eq.1) then 465 call transp_mat(TTmat, Tmat, 466 , nbfia, nbfja) 467 else 468 if(what.eq.3) then 469 call dscal(n3d*nbfia*nbfja,-1.0d0,TTmat,1) 470 endif 471 call transp_mat3(n3d,TTmat, Tmat, 472 , nbfia, nbfja) 473 endif 474 475 call endtimer(monitor_comp_transp) 476 477 478c Monitoring 479 480 call starttimer(monitor_comm_ga_acc2) 481 482 483 if(l3d) then 484 485 call dft_3dacc(g_zora, tmat, 486 & jrsh,jrsh2, 487 % jfirst, jlast, ifirst, ilast, nbfja) 488 else 489 if(dftnbacc) then 490 491 if(.not.nbfirst2) then 492 call ga_nbwait(nbhandl2) 493 endif 494 nbfirst2=.false. 495 call upd_atombl_nb(g_update(ii), 496 . basis,jat,iat,tmat,nbhandl2) 497 else 498 499 call upd_atom_block(g_update(ii),basis, 500 J jat,iat,tmat) 501 endif 502 endif 503c Monitoring 504 505 call endtimer(monitor_comm_ga_acc2) 506 507 508 endif 509 endif 510 511 endif 512 endif 513 514 168 continue 515 enddo ! ipol loop 516 endif 517 430 continue 518 500 continue 519 520 return 521 end 522c $Id$ 523