1#define STAG 1 2* $Id$ 3c Modified for construction of multiple sets of perturbed densities 4c and explicit derivative densities wrt nuclei 5c 6c BGJ - 8/98 7c 8C> \ingroup nwdft_xc 9C> @{ 10C> 11C> \file xc_rhogen.F 12C> Evaluate density and derivatives 13C> 14C> @} 15C> 16C> \ingroup nwdft_xc 17C> @{ 18C> 19C> \brief Evaluate the density and their derivatives on a grid 20C> 21C> 22 Subroutine xc_rhogen(what_in, 23 T tol_rho, basis, g_dens, max_at_bf, 24 N natoms, curatoms, ncuratoms, npert, 25 I ipol, nq, nbf, mbf, GRAD, ipol2, 26 & F, Pmat, ff, ffd, 27 C chi, delchi, heschi, 28 I ibf, iniz, ifin, 29 & rho, delrho, lap, rchi_atom, 30 & rdelchi_atom, rdens_atom, cetobfr, wmax, 31 & ttau, kske, dolap) 32 implicit none 33#include "global.fh" 34#include "errquit.fh" 35#include "mafdecls.fh" 36#include "dftpara.fh" 37#include "dist.fh" 38#include "dft_fdist.fh" 39 40 Logical GRAD !< [Input] .true. when using gradient corrected 41 !< functional 42 Logical kske !< [Input] .true. when using kinetic energy density 43 !< functional 44 integer what_in !< [Input] What to calculate 45 !< - what=0: calculate density 46 !< - what=1: calculate perturbed density (derivative 47 !< with respect to electron position) 48 !< - what=2: calculate derivative with respect to 49 !< nuclear coordinates 50 logical dolap !< [Input] .true. if the Laplacian of the 51 !< density is required 52c what=-1 dplot aka density contours 53c what=0 dens 54c what=1 pert 55c what=2 nucder 56 integer basis !< [Input] basis set handle 57 integer ipol !< [Input] no. of spin channels 58 integer ipol2 !< [Input] no. of spin channels in density 59 !< - 1 if closed shell 60 !< - 3 if open shell 61 integer nbf !< [Input] no. of basis functions 62 integer mbf !< [Input] "restricted" no. of basis functions 63 integer max_at_bf !< [Input] max no. bf per atom 64 integer nq !< [Input] no. of quadrature points 65 integer natoms !< [Input] no. of atoms 66 double precision wmax !< [Input] max weight 67 integer npert !< [Input] number of perturbed densities 68 integer curatoms(*) !< [Input] indexing array for "active" atoms 69 integer ncuratoms !< [Input] number of currently active atoms 70 integer g_dens(*) !< [Input] GA handle for DM 71 integer ibf(mbf) !< [Input] mapping of nbf_ao -> mbf 72 integer iniz(natoms) !< [Input] mapping of nbf_ao -> mbf 73 integer ifin(natoms) !< [Input] mapping of nbf_ao -> mbf 74 double precision tol_rho !< [Input] accuracy for rho evaluation 75 double precision chi(nq,mbf) !< [Input] function values 76 double precision delchi(nq,3,mbf)!< [Input] function gradients 77 double precision heschi(nq,6,mbf)!< [Input] function hessians 78c 79 double precision ttau(nq,ipol,*) !< [Output] Total Kohn-Sham K.E.density 80c 81 double precision lap(nq,ipol2,*) 82c 83 double precision delrho(nq,3,ipol,*) !< [Output] Derivative of density 84 double precision Pmat(*) !< [Scratch] scratch vector 85 double precision F(max_at_bf*max_at_bf) !< [Scratch] scratch vector 86 double precision ff(nq,*) !< [Scratch] scratch array 87 double precision ffd(nq,3,*) !< [Scratch] scratch array 88 double precision rho(nq,ipol2,*) !< [Output] The density 89 double precision rchi_atom(natoms) !< [Input] Screening parameters 90 double precision rdelchi_atom(natoms) !< [Input] Screening parameters 91 double precision rdens_atom(natoms,natoms,ipol) !< [Input] Screening parameters 92 integer cetobfr(2,natoms) !< [Input] Centers to basis functions 93c 94c local declarations 95c 96 integer i0, ii, mu, n, npol 97 integer ipert ! perturbation loop index 98 integer iat, inizia, ifirst, ilast, nbfia, nnia, iat_in 99 integer ifinia, ifinja 100 integer jat, inizja, jfirst, jlast, nbfja, nnja 101 integer iatcur, jatcur 102 double precision FUNC_MAX, DELFUNC_MAX, FUNC_MAXI, FUNC_MAXJ 103 double precision P_MAXJ, P_MAXJ_A, P_MAXJ_B, P_MAXIJ 104 double precision dabsmax 105 integer g_keepd(2) 106 integer nbhandl 107 integer jj 108 logical doffd,doitt 109 external dabsmax 110 external xc_rhoscreen 111 integer xc_rhoscreen 112 integer nonzero,nonz0 113 integer i_nz,l_nz 114 integer sizeblk, gindx 115 integer what 116 logical zapnegatives 117#ifdef DEBUG 118 integer ga_nodeid 119 external ga_nodeid 120#endif 121 g_keepd(1)=0 122 g_keepd(2)=0 123 call starttimer(monitor_xcrho) 124 if(what_in.ge.0) then 125 what=what_in 126 zapnegatives=.true. 127 else 128 what=-what_in 129 zapnegatives=.false. 130 endif 131c 132c Evaluate the charge density and its gradient at each of the 133c sampling points 134c 135 doffd=(what.eq.2.and.grad).or.(what.eq.0.and.kske).or. 136 & (what.eq.1.and.kske).or.(what.eq.0.and.dolap) 137 138 npol = (ipol*(ipol+1))/2 139c to keep compilers quiet 140 iatcur=1 141 jatcur=1 142c 143 if(what.eq.0) then 144 call dcopy(nq*npol,0d0,0,rho,1) 145 if (grad) call dcopy(3*nq*ipol,0d0,0,delrho,1) 146 if (kske) call dcopy(nq*ipol,0d0,0,ttau,1) ! total 147 if (dolap) call dcopy(nq*ipol2,0d0,0,lap,1) ! total 148c 149c repl stuff 150c 151 g_keepd(1)=g_dens(1) 152 if(ipol.eq.2) g_keepd(2)=g_dens(2) 153 if(xcreplicated.and.dorepdm) then 154 g_dens(1)=g_repdm(1) 155 if(ipol.eq.2) g_dens(2)=g_repdm(2) 156 endif 157 elseif(what.eq.1) then 158 call dcopy(nq*ipol*npert,0d0,0,rho,1) 159 if (grad) call dcopy(3*nq*ipol*npert,0d0,0,delrho,1) 160 if (kske) call dcopy(nq*ipol*npert,0d0,0,ttau,1) ! total 161 elseif(what.eq.2) then 162 call dcopy(nq*ipol*3*ncuratoms,0d0,0,rho,1) 163 if (grad)call dcopy(3*nq*ipol*3*ncuratoms,0d0,0,delrho,1) 164 else 165 call errquit('wrong what value',0,0) 166 endif 167c 168c Screening is accomplished by: p(r) <= |Xi(r)|*|Xj(r)|*|Dij| 169c Xi(r) is screened on desired accuracy/max(|Xj(r)|)*max(|Dij|) 170c Dij is screened on desired accuracy/max(|Xi(r)|)*max(|Xj(r)|) 171c 172 i0=ipol-1 173c 174 FUNC_MAX = dabsmax(natoms,rchi_atom) 175 DELFUNC_MAX=0d0 176 if(grad) DELFUNC_MAX = dabsmax(natoms,rdelchi_atom) 177c 178 nonzero=0 179 if(dftnbget) then 180 if (.not.ma_push_get 181 & (mt_int,(ipol*npert*natoms*(natoms+1))/2, 182 N 'nzmap map',l_nz,i_nz)) 183 & call errquit('xcrho:push_get failed', 13, MA_ERR) 184 nonzero=xc_rhoscreen(grad,ipol,natoms,npert, 185 I iniz, 186 W tol_rho,wmax, 187 O int_mb(i_nz), 188 R rchi_atom,rdelchi_atom,rdens_atom) 189 190 if(nonzero.eq.0) goto 1688 191c 192c prefetch first DM block 193c 194 nonz0=1 195 call xc_getdmblock(int_mb(i_nz),nonz0,natoms,cetobfr, 196 G g_dens(1), 197 A Pmat,nbhandl) 198 endif 199 200#ifdef STAG 201 do 230 iat_in = 1+ga_nodeid(), natoms+ga_nodeid() 202 iat=mod(iat_in,natoms) 203 if(iat.eq.0) iat=natoms 204#else 205 do 230 iat = 1, natoms 206#endif 207 inizia = iniz(iat) 208 if (inizia.eq.0)goto 230 209 if(what.eq.2) then 210 iatcur = curatoms(iat) 211 endif 212 ifinia = ifin(iat) 213 ifirst = cetobfr(1,iat) 214 ilast = cetobfr(2,iat) 215 nbfia = ilast-ifirst+1 216 nnia = ifinia-inizia+1 217c 218c screening parameters 219c 220 FUNC_MAXI = rchi_atom(iat) 221 if(grad) 222 . FUNC_MAXI = max(FUNC_MAXI,rdelchi_atom(iat)) 223 FUNC_MAXJ = max(FUNC_MAX,DELFUNC_MAX) 224 if(what.lt.2) then 225 if (ipol.gt.1)then 226 P_MAXJ_A = dabsmax(natoms,rdens_atom(1,iat,1)) 227 P_MAXJ_B = dabsmax(natoms,rdens_atom(1,iat,2)) 228 P_MAXJ = MAX(P_MAXJ_A, P_MAXJ_B) 229 else 230 P_MAXJ=0d0 231 do jat=1,iat 232 if(iniz(jat).ne.0) 233 . P_MAXJ=max(P_MAXJ,rdens_atom(jat,iat,1)) 234 enddo 235 endif 236 if (wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXJ.lt.tol_rho) goto 225 237 endif 238 do 220 jat = 1, iat 239 inizja = iniz(jat) 240 if (inizja.eq.0)goto 220 241 if(what.eq.2) then 242 jatcur = curatoms(jat) 243 if (jatcur.eq.0.and.iatcur.eq.0) goto 220 244 endif 245 call starttimer(monitor_rscreen0) 246 ifinja = ifin(jat) 247 jfirst = cetobfr(1,jat) 248 jlast = cetobfr(2,jat) 249 nbfja = jlast-jfirst+1 250 nnja = ifinja-inizja+1 251c 252c screening parameters 253c 254 FUNC_MAXJ=rchi_atom(jat) 255 if(grad) 256 . FUNC_MAXJ = max(FUNC_MAXJ,rdelchi_atom(jat)) 257 doitt=.true. 258 if(what.lt.2) then 259 P_MAXIJ = rdens_atom(iat,jat,1) 260 if(ipol.eq.2) P_MAXIJ = max(P_MAXIJ, 261 & rdens_atom(iat,jat,2)) 262 doitt=(wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXIJ.ge.tol_rho) 263 endif 264 call endtimer(monitor_rscreen0) 265 if (doitt) then 266 do ii = 1, ipol 267c 268c screening parameters 269c 270 if((what.gt.1).or.(rdens_atom(iat,jat,ii)* 271 R wmax*FUNC_MAXI*FUNC_MAXJ.ge.tol_rho)) then 272c 273c Loop over perturbations 274c 275 do 215 ipert = 1,npert 276 sizeblk=nbfia*nbfja 277 call updist(monitor_size_ga_get, sizeblk) 278c 279 if(dftnbget) then 280 call starttimer(monitor_wait3) 281 call ga_nbwait(nbhandl) 282 call endtimer(monitor_wait3) 283 call gat_mat(F, Pmat, nbfia, nbfja, nnia, nnja, 284 I ifirst, jfirst, ibf(inizia), ibf(inizja)) 285 nonz0=nonz0+1 286 if((npert*ipol).gt.1) then 287 gindx=ipert+(ii-1)*npert+1 288 289 if(gindx.gt.npert*ipol) 290 + gindx=mod(gindx,npert*ipol) 291 else 292 gindx=1 293 endif 294 if(nonz0.le.nonzero) 295 C call xc_getdmblock(int_mb(i_nz),nonz0,natoms, 296 C cetobfr,g_dens(gindx), 297 A Pmat,nbhandl) 298 if (wmax*FUNC_MAXI*FUNC_MAXJ* 299 . dabsmax(sizeblk,Pmat).lt.tol_rho)goto 215 300 else 301 call starttimer(monitor_gaget) 302 if(truerepdm) then 303 call xc_dmget( 304 + dbl_mb(k_repdm(ipert+(ii-1)*npert)),nbf_ld, 305 % ifirst, ilast, jfirst, jlast, Pmat,nbfia) 306 else 307 call ga_get(g_dens(ipert+(ii-1)*npert), 308 % ifirst, ilast, jfirst, jlast, Pmat,nbfia) 309 endif 310 call endtimer(monitor_gaget) 311 if (wmax*FUNC_MAXI*FUNC_MAXJ* 312 . dabsmax(sizeblk,Pmat).lt.tol_rho)goto 215 313 call gat_mat(F, Pmat, nbfia, nbfja, nnia, nnja, 314 + ifirst, jfirst, ibf(inizia), ibf(inizja)) 315 endif 316c 317 call starttimer(monitor_rhocomp) 318 if(iat.ne.jat) call dscal(nnia*nnja,2d0,F,1) 319c 320c Compute Xiat(r)*Xjat(r)*Diat,jat 321c 322 call dgemm('n','n',nq,nnja,nnia,1d0, 323 A chi(1,inizia),nq,F,nnia,0d0,ff,nq) 324 if(what.lt.2) then 325 jj=i0+ii 326 if(what.eq.1) jj=ii 327 do mu=inizja,ifinja 328 do n=1,nq 329 rho(n,jj,ipert) = rho(n,jj,ipert) + 330 P chi(n,mu)*ff(n,mu-inizja+1) 331 enddo 332 enddo 333 endif 334 if(doffd) then 335 call dgemm('n','n',nq*3,nnja,nnia,1d0, 336 A delchi(1,1,inizia),nq*3,F,nnia,0d0,ffd,nq*3) 337c 338c build tau for meta GGA 339c 340 if(kske) then 341 do mu=inizja,ifinja 342 do n=1,nq 343 ttau(n,ii,ipert) = ttau(n,ii,ipert)+0.5d0*( 344 & delchi(n,1,mu)*ffd(n,1,mu-inizja+1)+ 345 & delchi(n,2,mu)*ffd(n,2,mu-inizja+1)+ 346 & delchi(n,3,mu)*ffd(n,3,mu-inizja+1)) 347 enddo 348 enddo 349 endif 350 351 if (dolap) then 352 do mu=inizja,ifinja 353 do n=1,nq 354 ! total 355 lap(n,ii,ipert) = lap(n,ii,ipert) + 2d0* 356 & (delchi(n,1,mu)*ffd(n,1,mu-inizja+1) + 357 & delchi(n,2,mu)*ffd(n,2,mu-inizja+1) + 358 & delchi(n,3,mu)*ffd(n,3,mu-inizja+1)) 359cold & delchi(n,3,mu)*ffd(n,3,mu-inizja+1) + 360cold & heschi(n,1,mu)*ff(n,mu-inizja+1) + 361cold & heschi(n,4,mu)*ff(n,mu-inizja+1) + 362cold & heschi(n,6,mu)*ff(n,mu-inizja+1)) 363 enddo 364 enddo 365 endif 366 endif 367c 368 if((what.eq.2.and.jatcur.ne.0).or. 369c 370c We need the "grad" code at zero order in the 371c nuclear deriv case, but we can skip this part if 372c iat is not active 373c 374 O (what.ne.2.and.grad)) then 375 if(what.eq.0.or.what.eq.1) then 376 call xc_dchiff(nq,inizja,ifinja, 377 P delrho(1,1,ii,ipert),delchi, 378 F ff) 379cc DMR/Begin 380 if (dolap) then 381 do mu=inizja,ifinja 382 do n=1,nq 383 lap(n,ii,ipert) = lap(n,ii,ipert) + 384 & ff(n,mu-inizja+1)*(heschi(n,1,mu) + 385 & heschi(n,4,mu) + heschi(n,6,mu)) 386 enddo 387 enddo 388 endif 389cc DMR/End 390 elseif(what.eq.2) then 391 call xc_dchiffp(nq,ipol2,inizja,ifinja, 392 P rho,delchi,ff, 393 M ii,jat) 394 if (grad) then 395c 396c Compute nuclear gradient of delrho 397c 398 call xc_drhonuc(nq,ipol,inizja,ifinja, 399 D delrho,heschi,delchi,ff,ffd, 400 I ii,jat) 401 endif 402 endif 403 endif 404c 405 if((what.eq.2.and.iatcur.ne.0).or. 406 O (what.ne.2.and.grad)) then 407c 408c Compute delXiat(r)*Xjat(r)*Diat,jat 409c 410 call dgemm('n','t',nq,nnia,nnja,1d0, 411 A chi(1,inizja),nq,F,nnia,0d0,ff,nq) 412 if(what.lt.2) then 413 call xc_dchiff(nq,inizia,ifinia, 414 P delrho(1,1,ii,ipert),delchi, 415 F ff) 416cc DMR/Begin 417 if (dolap) then 418 do mu=inizia,ifinia 419 do n=1,nq 420 lap(n,ii,ipert) = lap(n,ii,ipert) + 421 & ff(n,mu-inizia+1)*(heschi(n,1,mu) + 422 & heschi(n,4,mu) + heschi(n,6,mu)) 423 enddo 424 enddo 425 endif 426cc DMR/End 427 elseif(what.eq.2) then 428 call xc_dchiffp(nq,ipol2,inizia,ifinia, 429 P rho,delchi,ff, 430 M ii,iat) 431 if(grad) then 432 call dgemm('n','t',nq*3,nnia,nnja,1d0, 433 A delchi(1,1,inizja),nq*3,F,nnia,0d0,ffd, 434 + nq*3) 435c 436c Compute nuclear gradient of delrho 437c 438 call xc_drhonuc(nq,ipol,inizia,ifinia, 439 D delrho,heschi,delchi,ff,ffd, 440 I ii,iat) 441 442 endif 443 endif 444 endif 445 call endtimer(monitor_rhocomp) 446 215 continue 447 endif 448 enddo 449 endif 450 220 continue 451 225 continue 452 230 continue 453 if(zapnegatives) then 454c 455c Enforce results that are compatible with the laws of physics. 456c I.e. Rho >= 0, Tau >= 0, Rho=0 ==> Grad(Rho)=0, and 457c Tau=0 ==> Grad(Rho)=0. The reason for doing this is that density 458c functionals have been designed to be valid for physically 459c sensible data points. They may fail spectacularly for data points 460c outside the physically sensible range. Alternatively to screening 461c the data here it might be done in the functionals themselves. 462c However, this requires a great deal of care to make sure that the 463c right limits are taken. 464c 465c Rho should be non-negative. Numerical errors in the construction 466c may leave small negative values that are formally incorrect. 467c Hence filter negative values out. 468c 469c Also note that rho.eq.0 implies delrho=0, so this is enforced as 470c well. Proof: Rho is a non-negative quantity. Hence if rho.eq.0 471c then rho is at a minimum. A minimum being an extremum it follows 472c that its gradient must be zero. QED. 473c 474 if (what.eq.0) then 475 do ii = 1, ipol 476 jj = i0 + ii 477 do n = 1, nq 478 if (rho(n,jj,1).le.0.0d0) then 479 rho(n,jj,1) = 0.0d0 480 if (grad) then 481 delrho(n,1,ii,1) = 0.0d0 482 delrho(n,2,ii,1) = 0.0d0 483 delrho(n,3,ii,1) = 0.0d0 484 endif 485 endif 486 enddo 487 enddo 488 else if (what.eq.1) then 489 do ii = 1, ipol 490 jj = ii 491 do n = 1, nq 492 if (rho(n,jj,npert).le.0.0d0) then 493 rho(n,jj,npert) = 0.0d0 494 if (grad) then 495 delrho(n,1,ii,npert) = 0.0d0 496 delrho(n,2,ii,npert) = 0.0d0 497 delrho(n,3,ii,npert) = 0.0d0 498 endif 499 endif 500 enddo 501 enddo 502 endif 503c 504c Tau should be non-negative. Numerical errors in the construction 505c may leave small negative values that can cause trouble in the 506c density functionals. Hence filter negative values out. 507c 508c Also note that tau.eq.0 implies delrho=0, so this is enforced as 509c well. Proof: tau =1/2 \sum_i (\nabla\phi_i)\cdot(\nabla\phi_i) 510c hence tau.eq.0 iff (\nabla\phi_i)=0 for all i. Delrho is 511c defined as delrho = \sum_i\nabla\phi_i^2 hence if the gradient 512c of the orbital is zero for all orbitals i then delrho must be 513c zero as well. QED. 514c 515 if (kske) then 516 if (what.eq.0) then 517 do ii = 1, ipol 518 do n = 1, nq 519 if (ttau(n,ii,1).le.0.0d0) then 520 ttau(n,ii,1) = 0.0d0 521 if (grad) then 522 delrho(n,1,ii,1) = 0.0d0 523 delrho(n,2,ii,1) = 0.0d0 524 delrho(n,3,ii,1) = 0.0d0 525 endif 526 endif 527 enddo 528 enddo 529 else if (what.eq.1) then 530 do ii = 1, ipol 531 do n = 1, nq 532 if (ttau(n,ii,npert+1).le.0.0d0) then 533 ttau(n,ii,npert+1) = 0.0d0 534 if (grad) then 535 delrho(n,1,ii,npert+1) = 0.0d0 536 delrho(n,2,ii,npert+1) = 0.0d0 537 delrho(n,3,ii,npert+1) = 0.0d0 538 endif 539 endif 540 enddo 541 enddo 542 endif 543 endif 544 endif 545c 546 call starttimer(monitor_rhocomp2) 547 if(what.eq.0) then 548c 549c Only construct total densities for regular case 550c 551 if (ipol.eq.2)then 552 call dcopy(nq, rho(1,2,1), 1, rho(1,1,1), 1) 553 call daxpy(nq, 1.d0, rho(1,3,1), 1, rho(1,1,1), 1) 554 endif 555 endif 556 if(what.eq.0) then 557 if(xcreplicated.and.dorepdm) then 558 g_dens(1)=g_keepd(1) 559 if(ipol.eq.2) g_dens(2)=g_keepd(2) 560 endif 561 endif 562 call endtimer(monitor_rhocomp2) 563c 564 1688 continue 565 if(dftnbget) then 566 if (.not.ma_pop_stack(l_nz)) 567 & call errquit('xcrho:pop_stack failed', 13, MA_ERR) 568 endif 569 570 call endtimer(monitor_xcrho) 571 return 572 end 573c 574 subroutine xc_drhonuc(nq,ipol,n0,n1, 575 D delrho,heschi,delchi,ff,ffd, 576 I ii,iat) 577 implicit none 578 integer nq,ipol,ii,n0,n1,iat 579 double precision delrho(nq,3,ipol,3,*) 580 double precision heschi(nq,6,*) 581 double precision delchi(nq,3,*) 582 double precision ff(nq,*) 583 double precision ffd(nq,3,*) 584c 585 integer n,mu,mu1 586c 587 do mu=n0,n1 588 mu1=mu-n0+1 589 do n = 1, nq 590 delrho(n,1,ii,1,iat) = delrho(n,1,ii,1,iat) - 591 & heschi(n,1,mu)*ff(n,mu1) - 592 - delchi(n,1,mu)*ffd(n,1,mu1) 593 delrho(n,2,ii,1,iat) = delrho(n,2,ii,1,iat) - 594 & heschi(n,2,mu)*ff(n,mu1) - 595 - delchi(n,1,mu)*ffd(n,2,mu1) 596 delrho(n,3,ii,1,iat) = delrho(n,3,ii,1,iat) - 597 & heschi(n,3,mu)*ff(n,mu1) - 598 - delchi(n,1,mu)*ffd(n,3,mu1) 599 delrho(n,1,ii,2,iat) = delrho(n,1,ii,2,iat) - 600 & heschi(n,2,mu)*ff(n,mu1) - 601 - delchi(n,2,mu)*ffd(n,1,mu1) 602 delrho(n,2,ii,2,iat) = delrho(n,2,ii,2,iat) - 603 & heschi(n,4,mu)*ff(n,mu1) - 604 - delchi(n,2,mu)*ffd(n,2,mu1) 605 delrho(n,3,ii,2,iat) = delrho(n,3,ii,2,iat) - 606 & heschi(n,5,mu)*ff(n,mu1) - 607 - delchi(n,2,mu)*ffd(n,3,mu1) 608 delrho(n,1,ii,3,iat) = delrho(n,1,ii,3,iat) - 609 & heschi(n,3,mu)*ff(n,mu1) - 610 - delchi(n,3,mu)*ffd(n,1,mu1) 611 delrho(n,2,ii,3,iat) = delrho(n,2,ii,3,iat) - 612 & heschi(n,5,mu)*ff(n,mu1) - 613 - delchi(n,3,mu)*ffd(n,2,mu1) 614 delrho(n,3,ii,3,iat) = delrho(n,3,ii,3,iat) - 615 & heschi(n,6,mu)*ff(n,mu1) - 616 - delchi(n,3,mu)*ffd(n,3,mu1) 617 enddo 618 enddo 619 return 620 end 621 subroutine xc_dchiff(nq,n0,n1, 622 P delrho,delchi,ff) 623 implicit none 624 integer nq,n0,n1 625 double precision delrho(nq,3) 626 double precision delchi(nq,3,*) 627 double precision ff(nq,*) 628c 629 integer n,mu,mu1 630c 631 do mu=n0,n1 632 mu1=mu-n0+1 633 do n = 1, nq 634 delrho(n,1) = delrho(n,1) + delchi(n,1,mu)*ff(n,mu1) 635 delrho(n,2) = delrho(n,2) + delchi(n,2,mu)*ff(n,mu1) 636 delrho(n,3) = delrho(n,3) + delchi(n,3,mu)*ff(n,mu1) 637 enddo 638 enddo 639 return 640 end 641 subroutine xc_dchiffp(nq,ipol2,n0,n1, 642 P rho,delchi,ff, 643 M ii,iat) 644 implicit none 645 integer nq,ipol2,n0,n1 646 double precision rho(nq,ipol2,3,*) 647 double precision delchi(nq,3,*) 648 double precision ff(nq,*) 649 integer ii,iat 650c 651 integer n,mu,mu1 652c 653 do mu=n0,n1 654 mu1=mu-n0+1 655 do n = 1, nq 656 rho(n,ii,1,iat) = rho(n,ii,1,iat)-delchi(n,1,mu)*ff(n,mu1) 657 rho(n,ii,2,iat) = rho(n,ii,2,iat)-delchi(n,2,mu)*ff(n,mu1) 658 rho(n,ii,3,iat) = rho(n,ii,3,iat)-delchi(n,3,mu)*ff(n,mu1) 659 enddo 660 enddo 661 return 662 end 663 664c@@@@@@@@@@@@@@@@@ 665 integer function xc_rhoscreen(grad,ipol,natoms,npert, 666 I iniz, 667 W tol_rho,wmax, 668 O nz, 669 R rchi_atom,rdelchi_atom,rdens_atom) 670 implicit none 671#include "global.fh" 672 logical grad 673 integer ipol,natoms,npert 674 integer iniz(natoms) ! mapping of nbf_ao -> mbf 675 double precision tol_rho,wmax 676 double precision rchi_atom(*) 677 double precision rdelchi_atom(*) 678 double precision rdens_atom(natoms,natoms,ipol) 679c 680c output 681c 682 integer nz(*) 683c 684 integer iat,jat,inizia,inizja,iat_in 685 integer ipert 686 double precision FUNC_MAXI,FUNC_MAXJ,DELFUNC_MAX,FUNC_MAX, 687 , P_MAXJ,P_MAXJ_A,P_MAXJ_B,P_MAXIJ 688 integer nonzero,ii 689 double precision dabsmax 690 external dabsmax 691c 692 FUNC_MAX = dabsmax(natoms,rchi_atom) 693 DELFUNC_MAX=0d0 694 if(grad) DELFUNC_MAX = dabsmax(natoms,rdelchi_atom) 695 696 nonzero=0 697cscat 698#ifdef STAG 699 do iat_in = 1+ga_nodeid(), natoms+ga_nodeid() 700 iat=mod(iat_in,natoms) 701 if(iat.eq.0) iat=natoms 702#else 703 do iat = 1, natoms 704#endif 705 inizia = iniz(iat) 706 if (inizia.ne.0) then 707c 708c screening parameters 709c 710 FUNC_MAXI = rchi_atom(iat) 711 if(grad) 712 . FUNC_MAXI = max(FUNC_MAXI,rdelchi_atom(iat)) 713 FUNC_MAXJ = max(FUNC_MAX,DELFUNC_MAX) 714 if (ipol.gt.1)then 715 P_MAXJ_A = dabsmax(natoms,rdens_atom(1,iat,1)) 716 P_MAXJ_B = dabsmax(natoms,rdens_atom(1,iat,2)) 717 P_MAXJ = MAX(P_MAXJ_A, P_MAXJ_B) 718 else 719 P_MAXJ=0d0 720 do jat=1,iat 721 if(iniz(jat).ne.0) 722 . P_MAXJ=max(P_MAXJ,rdens_atom(jat,iat,1)) 723 enddo 724 endif 725 if (wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXJ.ge.tol_rho) then 726 do jat = 1, iat 727 inizja = iniz(jat) 728 if (inizja.ne.0) then 729c 730c screening parameters 731c 732 FUNC_MAXJ=rchi_atom(jat) 733 if(grad) 734 . FUNC_MAXJ = max(FUNC_MAXJ,rdelchi_atom(jat)) 735 P_MAXIJ = rdens_atom(iat,jat,1) 736 if(ipol.eq.2) P_MAXIJ = max(P_MAXIJ, 737 & rdens_atom(iat,jat,2)) 738 if(wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXIJ.ge.tol_rho)then 739c 740 do ipert = 1,npert 741 do ii = 1, ipol 742c 743c screening parameters 744c 745 P_MAXIJ = rdens_atom(iat,jat,ii) 746 if (wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXIJ. 747 P ge.tol_rho) then 748 nonzero=nonzero+1 749 nz(nonzero)=(iat-1)*natoms+jat 750 endif 751 enddo 752 enddo 753 endif 754 endif 755 enddo ! loop jat 756 endif 757 endif 758 enddo ! loop iat 759 xc_rhoscreen=nonzero 760 return 761 end 762 subroutine xc_getdmblock(nz,nonz0,natoms,cetobfr,g_dens, 763 A Pmat,nbhandl) 764 implicit none 765#include "dist.fh" 766#include "dft_fdist.fh" 767 integer nz(*) 768 integer natoms 769 integer nbhandl 770 integer g_dens 771 integer cetobfr(2,*) 772 integer nonz0 773 double precision pmat(*) 774c 775 integer ij0 776 integer iat0,jat0,ifirst,ilast,jfirst,jlast,nbfia 777c 778 ij0=nz(nonz0) 779 iat0=(ij0-1)/natoms+1 780 jat0=ij0-(iat0-1)*natoms 781 ifirst = cetobfr(1,iat0) 782 ilast = cetobfr(2,iat0) 783 jfirst = cetobfr(1,jat0) 784 jlast = cetobfr(2,jat0) 785 nbfia= ilast-ifirst+1 786 call starttimer(monitor_ganbget) 787 call ga_nbget(g_dens, 788 % ifirst, ilast, jfirst, jlast, Pmat,nbfia,nbhandl) 789 call endtimer(monitor_ganbget) 790 return 791 end 792 subroutine xc_dmget(repdm, nbf_ld, 793 % ilo, ihi, jlo, jhi, Pmat,nbfia) 794 implicit none 795 integer ilo, ihi, jlo, jhi, nbfia,nbf_ld 796 double precision pmat(nbfia,*) 797 double precision repdm(*) 798 integer ij,nnn 799 integer i,j 800c 801 nnn=nbfia 802 if(ilo.ne.jlo) then 803 do j=jlo,jhi 804 ij=((j-1)*(2*(nbf_ld+1)-j)+1)/2+ilo-j+1 805 call dcopy(nnn,repdm(ij),1,pmat(1,j-jlo+1),1) 806 enddo 807 else 808c diag block: copy only lower tr 809 do j=jlo,jhi 810 ij=((j-1)*(2*(nbf_ld+1)-j)+1)/2+1 811 call dcopy(nnn,repdm(ij),1,pmat(j-jlo+1,j-jlo+1),1) 812 nnn=nnn-1 813 enddo 814c copy offdiag terms (aka transp) 815 do j=1,nbfia 816 do i=j+1,nbfia 817 pmat(j,i)=pmat(i,j) 818 enddo 819 enddo 820 endif 821 return 822 end 823c 824c Essentially a clone of the old xc_rhogen.F without the new negative density screening 825c We need this for transition densities which can be -ve 826c 827 subroutine td_rhogen(what, 828 T tol_rho, basis, g_dens, max_at_bf, 829 N natoms, curatoms, ncuratoms, npert, 830 I ipol, nq, nbf, mbf, GRAD, ipol2, 831 & F, Pmat, ff, ffd, 832 C chi, delchi, heschi, 833 I ibf, iniz, ifin, 834 & rho, delrho, lap, rchi_atom, 835 & rdelchi_atom, rdens_atom, cetobfr, wmax, 836 & ttau, kske, dolap) 837c 838 implicit none 839c 840#include "errquit.fh" 841#include "mafdecls.fh" 842#include "dftpara.fh" 843#include "dist.fh" 844#include "dft_fdist.fh" 845 846 Logical GRAD !< [Input] .true. when using gradient corrected 847 !< functional 848 Logical kske !< [Input] .true. when using kinetic energy density 849 !< functional 850 integer what !< [Input] What to calculate 851 !< - what=0: calculate density 852 !< - what=1: calculate perturbed density (derivative 853 !< with respect to electron position) 854 !< - what=2: calculate derivative with respect to 855 !< nuclear coordinates 856 logical dolap ! true if lap is required [in] 857 858c what=0 dens 859c what=1 pert 860c what=2 nucder 861 integer basis !< [Input] basis set handle 862 integer ipol !< [Input] no. of spin channels 863 integer ipol2 !< [Input] no. of spin channels in density 864 !< - 1 if closed shell 865 !< - 3 if open shell 866 integer nbf !< [Input] no. of basis functions 867 integer mbf !< [Input] "restricted" no. of basis functions 868 integer max_at_bf !< [Input] max no. bf per atom 869 integer nq !< [Input] no. of quadrature points 870 integer natoms !< [Input] no. of atoms 871 double precision wmax !< [Input] max weight 872 integer npert !< [Input] number of perturbed densities 873 integer curatoms(*) !< [Input] indexing array for "active" atoms 874 integer ncuratoms !< [Input] number of currently active atoms 875 integer g_dens(*) !< [Input] GA handle for DM 876 integer ibf(mbf) !< [Input] mapping of nbf_ao -> mbf 877 integer iniz(natoms) !< [Input] mapping of nbf_ao -> mbf 878 integer ifin(natoms) !< [Input] mapping of nbf_ao -> mbf 879 double precision tol_rho !< [Input] accuracy for rho evaluation 880 double precision chi(nq,mbf) !< [Input] function values 881 double precision delchi(nq,3,mbf)!< [Input] function gradients 882 double precision heschi(nq,6,mbf)!< [Input] function hessians 883c 884 double precision ttau(nq,ipol,*) !< [Output] Total Kohn-Sham K.E.density 885c 886 double precision lap(nq,ipol2,*) 887c 888 double precision delrho(nq,3,ipol,*) !< [Output] Derivative of density 889 double precision Pmat(*) !< [Scratch] scratch vector 890 double precision F(max_at_bf*max_at_bf) !< [Scratch] scratch vector 891 double precision ff(nq,*) !< [Scratch] scratch array 892 double precision ffd(nq,3,*) !< [Scratch] scratch array 893 double precision rho(nq,ipol2,*) !< [Output] The density 894 double precision rchi_atom(natoms) !< [Input] Screening parameters 895 double precision rdelchi_atom(natoms) !< [Input] Screening parameters 896 double precision rdens_atom(natoms,natoms,ipol) !< [Input] Screening parameters 897 integer cetobfr(2,natoms) !< [Input] Centers to basis functions 898c 899c local declarations 900c 901 integer i0, ii, mu, n, npol 902 integer ipert ! perturbation loop index 903 integer iat, inizia, ifirst, ilast, nbfia, nnia 904 integer ifinia, ifinja 905 integer jat, inizja, jfirst, jlast, nbfja, nnja 906 integer iatcur, jatcur 907 double precision FUNC_MAX, DELFUNC_MAX, FUNC_MAXI, FUNC_MAXJ 908 double precision P_MAXJ, P_MAXJ_A, P_MAXJ_B, P_MAXIJ 909 double precision dabsmax 910 integer g_keepd(2) 911 integer nbhandl 912 integer jj 913 logical doffd,doitt 914 external dabsmax 915 external xc_rhoscreen 916 integer xc_rhoscreen 917 integer nonzero,nonz0 918 integer i_nz,l_nz 919 integer sizeblk, gindx 920#ifdef DEBUG 921 integer ga_nodeid 922 external ga_nodeid 923#endif 924 call starttimer(monitor_xcrho) 925 g_keepd(1)=0 926 g_keepd(2)=0 927c 928c Evaluate the charge density and its gradient at each of the 929c sampling points 930c 931 doffd=(what.eq.2.and.grad).or.(what.eq.0.and.kske).or. 932 & (what.eq.1.and.kske).or.(what.eq.0.and.dolap) 933 934 npol = (ipol*(ipol+1))/2 935c to keep compilers quiet 936 iatcur=1 937 jatcur=1 938c 939 if(what.eq.0) then 940 call dcopy(nq*npol,0.D0,0,rho,1) 941 if (grad) call dcopy(3*nq*ipol,0.D0,0,delrho,1) 942 if (kske) call dcopy(nq*ipol,0.D0,0,ttau,1) ! total 943 if (dolap) call dfill(nq*ipol2,0.D0,lap,1) ! total 944c 945c repl stuff 946c 947 g_keepd(1)=g_dens(1) 948 if(ipol.eq.2) g_keepd(2)=g_dens(2) 949 if(xcreplicated.and.dorepdm) then 950 g_dens(1)=g_repdm(1) 951 if(ipol.eq.2) g_dens(2)=g_repdm(2) 952 endif 953 elseif(what.eq.1) then 954 call dfill(nq*ipol*npert,0.D0,rho,1) 955 if (grad) call dfill(3*nq*ipol*npert,0.D0,delrho,1) 956 if (kske) call dfill(nq*ipol*npert,0.D0,ttau,1) ! total 957 elseif(what.eq.2) then 958 call dfill(nq*ipol*3*ncuratoms,0.D0,rho,1) 959 if (grad)call dfill(3*nq*ipol*3*ncuratoms,0.D0,delrho,1) 960 else 961 call errquit('wrong what value',0,0) 962 endif 963c 964c Screening is accomplished by: p(r) <= |Xi(r)|*|Xj(r)|*|Dij| 965c Xi(r) is screened on desired accuracy/max(|Xj(r)|)*max(|Dij|) 966c Dij is screened on desired accuracy/max(|Xi(r)|)*max(|Xj(r)|) 967c 968 i0=ipol-1 969c 970 FUNC_MAX = dabsmax(natoms,rchi_atom) 971 DELFUNC_MAX=0d0 972 if(grad) DELFUNC_MAX = dabsmax(natoms,rdelchi_atom) 973c 974 nonzero=0 975 if(dftnbget) then 976 if (.not.ma_push_get 977 & (mt_int,(ipol*npert*natoms*(natoms+1))/2, 978 N 'nzmap map',l_nz,i_nz)) 979 & call errquit('xcrho:push_get failed', 13, MA_ERR) 980 nonzero=xc_rhoscreen(grad,ipol,natoms,npert, 981 I iniz, 982 W tol_rho,wmax, 983 O int_mb(i_nz), 984 R rchi_atom,rdelchi_atom,rdens_atom) 985 986 if(nonzero.eq.0) goto 1688 987c 988c prefetch first DM block 989c 990 nonz0=1 991 call xc_getdmblock(int_mb(i_nz),nonz0,natoms,cetobfr, 992 G g_dens(1), 993 A Pmat,nbhandl) 994 endif 995 996 do 230 iat = 1, natoms 997 inizia = iniz(iat) 998 if (inizia.eq.0)goto 230 999 if(what.eq.2) then 1000 iatcur = curatoms(iat) 1001 endif 1002 ifinia = ifin(iat) 1003 ifirst = cetobfr(1,iat) 1004 ilast = cetobfr(2,iat) 1005 nbfia = ilast-ifirst+1 1006 nnia = ifinia-inizia+1 1007c 1008c screening parameters 1009c 1010 FUNC_MAXI = rchi_atom(iat) 1011 if(grad) 1012 . FUNC_MAXI = max(FUNC_MAXI,rdelchi_atom(iat)) 1013 FUNC_MAXJ = max(FUNC_MAX,DELFUNC_MAX) 1014 if(what.lt.2) then 1015 if (ipol.gt.1)then 1016 P_MAXJ_A = dabsmax(natoms,rdens_atom(1,iat,1)) 1017 P_MAXJ_B = dabsmax(natoms,rdens_atom(1,iat,2)) 1018 P_MAXJ = MAX(P_MAXJ_A, P_MAXJ_B) 1019 else 1020 P_MAXJ=0d0 1021 do jat=1,iat 1022 if(iniz(jat).ne.0) 1023 . P_MAXJ=max(P_MAXJ,rdens_atom(jat,iat,1)) 1024 enddo 1025 endif 1026 if (wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXJ.lt.tol_rho) goto 225 1027 endif 1028 do 220 jat = 1, iat 1029 inizja = iniz(jat) 1030 if (inizja.eq.0)goto 220 1031 if(what.eq.2) then 1032 jatcur = curatoms(jat) 1033 if (jatcur.eq.0.and.iatcur.eq.0) goto 220 1034 endif 1035 call starttimer(monitor_rscreen0) 1036 ifinja = ifin(jat) 1037 jfirst = cetobfr(1,jat) 1038 jlast = cetobfr(2,jat) 1039 nbfja = jlast-jfirst+1 1040 nnja = ifinja-inizja+1 1041c 1042c screening parameters 1043c 1044 FUNC_MAXJ=rchi_atom(jat) 1045 if(grad) 1046 . FUNC_MAXJ = max(FUNC_MAXJ,rdelchi_atom(jat)) 1047 doitt=.true. 1048 if(what.lt.2) then 1049 P_MAXIJ = rdens_atom(iat,jat,1) 1050 if(ipol.eq.2) P_MAXIJ = max(P_MAXIJ, 1051 & rdens_atom(iat,jat,2)) 1052 doitt=(wmax*FUNC_MAXI*FUNC_MAXJ*P_MAXIJ.ge.tol_rho) 1053 endif 1054 call endtimer(monitor_rscreen0) 1055 if (doitt) then 1056 do ii = 1, ipol 1057c 1058c screening parameters 1059c 1060 if((what.gt.1).or.(rdens_atom(iat,jat,ii)* 1061 R wmax*FUNC_MAXI*FUNC_MAXJ.ge.tol_rho)) then 1062c 1063c Loop over perturbations 1064c 1065 do 215 ipert = 1,npert 1066 sizeblk=nbfia*nbfja 1067 call updist(monitor_size_ga_get, sizeblk) 1068c 1069 if(dftnbget) then 1070 call starttimer(monitor_wait3) 1071 call ga_nbwait(nbhandl) 1072 call endtimer(monitor_wait3) 1073 call gat_mat(F, Pmat, nbfia, nbfja, nnia, nnja, 1074 I ifirst, jfirst, ibf(inizia), ibf(inizja)) 1075 nonz0=nonz0+1 1076 if((npert*ipol).gt.1) then 1077 gindx=ipert+(ii-1)*npert+1 1078 1079 if(gindx.gt.npert*ipol)gindx=mod(gindx,npert*ipol) 1080 else 1081 gindx=1 1082 endif 1083 if(nonz0.le.nonzero) 1084 C call xc_getdmblock(int_mb(i_nz),nonz0,natoms, 1085 C cetobfr,g_dens(gindx), 1086 A Pmat,nbhandl) 1087 if (wmax*FUNC_MAXI*FUNC_MAXJ* 1088 . dabsmax(sizeblk,Pmat).lt.tol_rho)goto 215 1089 else 1090 call starttimer(monitor_gaget) 1091 if(truerepdm) then 1092 call xc_dmget(dbl_mb(k_repdm(ipert+(ii-1)*npert)), 1093 & nbf_ld, 1094 % ifirst, ilast, jfirst, jlast, Pmat,nbfia) 1095 else 1096 call ga_get(g_dens(ipert+(ii-1)*npert), 1097 % ifirst, ilast, jfirst, jlast, Pmat,nbfia) 1098 endif 1099 call endtimer(monitor_gaget) 1100 if (wmax*FUNC_MAXI*FUNC_MAXJ* 1101 . dabsmax(sizeblk,Pmat).lt.tol_rho)goto 215 1102 call gat_mat(F, Pmat, nbfia, nbfja, nnia, nnja, ifirst, 1103 & jfirst, ibf(inizia), ibf(inizja)) 1104 endif 1105 1106c 1107 call starttimer(monitor_rhocomp) 1108 if(iat.ne.jat) call dscal(nnia*nnja,2d0,F,1) 1109c 1110c Compute Xiat(r)*Xjat(r)*Diat,jat 1111c 1112 call dgemm('n','n',nq,nnja,nnia,1d0, 1113 A chi(1,inizia),nq,F,nnia,0d0,ff,nq) 1114 if(what.lt.2) then 1115 jj=i0+ii 1116 if(what.eq.1) jj=ii 1117 do mu=inizja,ifinja 1118 do n=1,nq 1119 rho(n,jj,ipert) = rho(n,jj,ipert) + 1120 P chi(n,mu)*ff(n,mu-inizja+1) 1121 enddo 1122 enddo 1123 endif 1124 if(doffd) then 1125 call dgemm('n','n',nq*3,nnja,nnia,1d0, 1126 A delchi(1,1,inizia),nq*3,F,nnia,0d0,ffd,nq*3) 1127c 1128c build tau for meta GGA 1129c 1130 if(kske) then 1131 do mu=inizja,ifinja 1132 do n=1,nq 1133 ttau(n,ii,ipert) = ttau(n,ii,ipert)+0.5d0*( 1134 & delchi(n,1,mu)*ffd(n,1,mu-inizja+1)+ 1135 & delchi(n,2,mu)*ffd(n,2,mu-inizja+1)+ 1136 & delchi(n,3,mu)*ffd(n,3,mu-inizja+1)) 1137 enddo 1138 enddo 1139 endif 1140 1141 if (dolap) then 1142 do mu=inizja,ifinja 1143 do n=1,nq 1144 ! total 1145 lap(n,ii,ipert) = lap(n,ii,ipert) + 2d0* 1146 & (delchi(n,1,mu)*ffd(n,1,mu-inizja+1) + 1147 & delchi(n,2,mu)*ffd(n,2,mu-inizja+1) + 1148 & delchi(n,3,mu)*ffd(n,3,mu-inizja+1) + 1149 & heschi(n,1,mu)*ff(n,mu-inizja+1) + 1150 & heschi(n,4,mu)*ff(n,mu-inizja+1) + 1151 & heschi(n,6,mu)*ff(n,mu-inizja+1)) 1152 enddo 1153 enddo 1154 endif 1155 endif 1156c 1157 if((what.eq.2.and.jatcur.ne.0).or. 1158c 1159c We need the "grad" code at zero order in the nuclear deriv case, 1160c but we can skip this part if iat is not active 1161c 1162 O (what.ne.2.and.grad)) then 1163 if(what.eq.0.or.what.eq.1) then 1164 call xc_dchiff(nq,inizja,ifinja, 1165 P delrho(1,1,ii,ipert),delchi, 1166 F ff) 1167 elseif(what.eq.2) then 1168 call xc_dchiffp(nq,ipol2,inizja,ifinja, 1169 P rho,delchi,ff, 1170 M ii,jat) 1171 if (grad) then 1172c 1173c Compute nuclear gradient of delrho 1174c 1175 call xc_drhonuc(nq,ipol,inizja,ifinja, 1176 D delrho,heschi,delchi,ff,ffd, 1177 I ii,jat) 1178 endif 1179 endif 1180 endif 1181c 1182 if((what.eq.2.and.iatcur.ne.0).or. 1183 O (what.ne.2.and.grad)) then 1184c 1185c Compute delXiat(r)*Xjat(r)*Diat,jat 1186c 1187 call dgemm('n','t',nq,nnia,nnja,1d0, 1188 A chi(1,inizja),nq,F,nnia,0d0,ff,nq) 1189 if(what.lt.2) then 1190 call xc_dchiff(nq,inizia,ifinia, 1191 P delrho(1,1,ii,ipert),delchi, 1192 F ff) 1193 elseif(what.eq.2) then 1194 call xc_dchiffp(nq,ipol2,inizia,ifinia, 1195 P rho,delchi,ff, 1196 M ii,iat) 1197 if(grad) then 1198 call dgemm('n','t',nq*3,nnia,nnja,1d0, 1199 A delchi(1,1,inizja),nq*3,F,nnia,0d0,ffd,nq*3) 1200c 1201c Compute nuclear gradient of delrho 1202c 1203 call xc_drhonuc(nq,ipol,inizia,ifinia, 1204 D delrho,heschi,delchi,ff,ffd, 1205 I ii,iat) 1206 1207 endif 1208 endif 1209 endif 1210 call endtimer(monitor_rhocomp) 1211 215 continue 1212 endif 1213 enddo 1214 endif 1215 220 continue 1216 225 continue 1217 230 continue 1218 1219 call starttimer(monitor_rhocomp2) 1220 if(what.eq.0) then 1221c 1222c Only construct total densities for regular case 1223c 1224 if (ipol.eq.2)then 1225 call dcopy(nq, rho(1,2,1), 1, rho(1,1,1), 1) 1226 call daxpy(nq, 1.d0, rho(1,3,1), 1, rho(1,1,1), 1) 1227 endif 1228 endif 1229 if(what.eq.0) then 1230 if(xcreplicated.and.dorepdm) then 1231 g_dens(1)=g_keepd(1) 1232 if(ipol.eq.2) g_dens(2)=g_keepd(2) 1233 endif 1234 endif 1235 call endtimer(monitor_rhocomp2) 1236c 1237 1688 continue 1238 if(dftnbget) then 1239 if (.not.ma_pop_stack(l_nz)) 1240 & call errquit('xcrho:pop_stack failed', 13, MA_ERR) 1241 endif 1242 1243 call endtimer(monitor_xcrho) 1244 1245 return 1246 end 1247 1248 subroutine xc_delrhodotr(nq,ipol,qxyz,rho,delrho,delrhoq) 1249 implicit none 1250 integer nq,ipol ![in] 1251 double precision rho(nq,*) ! [in] 1252 double precision delrho(nq,3,*) ! [in] 1253 double precision qxyz(3,*) ! [in] 1254 double precision delrhoq(nq,*) ! [out] 1255c 1256 integer ii,n 1257 1258 if(ipol.eq.1) then 1259 do n = 1, nq 1260 delrhoq(n,1) = 3d0*rho(n,1)+ 1261 D delrho(n,1,1)*qxyz(1,n) + 1262 D delrho(n,2,1)*qxyz(2,n) + 1263 D delrho(n,3,1)*qxyz(3,n) 1264 enddo 1265 else 1266 do ii=2,3 1267 do n = 1, nq 1268 delrhoq(n,ii-1) = 3d0*rho(n,ii)+ 1269 D delrho(n,1,ii-1)*qxyz(1,n) + 1270 D delrho(n,2,ii-1)*qxyz(2,n) + 1271 D delrho(n,3,ii-1)*qxyz(3,n) 1272 enddo 1273 enddo 1274 endif 1275 return 1276 end 1277C> 1278C> @} 1279 1280 1281 1282 1283