1c 2c Computes explicit nuclear 2nd derivatives of the XC energy 3c 4c BGJ - 8/98 5c 6c $Id$ 7c 8C> \ingroup hess 9C> @{ 10C> 11C> \file xc_d2expl.F 12C> Explicit 2nd derivatives wrt nuclear coordinates of the DFT 13C> functional 14C> 15C> \brief Computes explicit nuclear 2nd derivatives of the XC energy 16C> 17C> The explicit 2nd derivatives of the XC energy with respect to the 18C> nuclear coordinates are defined as 19C> \f{eqnarray*}{ 20C> \frac{\mathrm{d}^2E_{xc}}{\mathrm{d}x\mathrm{d}y} &=& \sum_{ij}\int 21C> \frac{\mathrm{d}^2f_{xc}}{\mathrm{d}p_i\mathrm{d}p_j} 22C> \frac{\mathrm{d}p_i}{\mathrm{d}x}\frac{\mathrm{d}p_j}{\mathrm{d}y} 23C> \mathrm{d}r 24C> + \sum_i\int\frac{\mathrm{d}f_{xc}}{\mathrm{d}p_i} 25C> \frac{\mathrm{d}^2p_i}{\mathrm{d}x\mathrm{d}y} 26C> \mathrm{d}r 27C> \f} 28C> where \f$p \in \{\rho_\alpha,\rho_\beta,\gamma_{\alpha\alpha}, 29C> \gamma_{\alpha\beta},\gamma_{\beta\beta}\}\f$ and the indeces 30C> \f$i\f$ and \f$j\f$ run over this set. 31c 32c d2Exc / d2fxc dp(i) dp(j) / dfxc d2p(i) 33c ----- = sum sum | ----------- ----- ----- + sum | ----- ------ 34c dx dy i j / dp(i) dp(j) dx dy i / dp(i) dx dy 35c 36c where "p" represents a density parameter in the set 37c { ra, rb, gaa, gab, gbb }, and the i, j indices run over this set 38c 39c First step: separable part of 2nd derivative contribution 40c 41 Subroutine xc_d2expl(tol_rho, scr, 42 & Amat, Amat2, Acof2, Cmat, Cmat2, Ccof2, Mmat, Mmat2, Mcof2, 43 & F, Pmat, ff, s, 44 & chi, delchi, heschi, d3chi, 45 & curatoms, ncuratoms, ipol, nq, nbf, max_at_bf, GRAD, basis, 46 & natoms, iniz, ifin, drho, ddelrho, dttau, 47 & delrho, g_dens, hess, ibf, 48 & rchi_atom, rdelchi_atom, rdens_atom, cetobfr, kske) 49c 50 implicit none 51#include "dft2drv.fh" 52#include "stdio.fh" 53#include "util.fh" 54#include "global.fh" 55c 56 logical kske !< [Input] .True. if the functional is a meta-GGA 57 integer basis !< [Input] The basis set handle 58 integer max_at_bf !< [Input] The maximum number of basis functions 59 !< on any atom 60 integer ipol !< [Input] The number of spin channels 61 integer nq !< [Input] The number of grid points 62 integer nbf !< [Input] The number of basis functions 63 integer natoms !< [Input] The number of atoms 64 integer ncuratoms !< [Input] The number of currently "active" 65 !< atoms 66 integer curatoms(natoms) !< [Input] The mapping array for 67 !< currently active atoms 68 double precision tol_rho !< [Input] The tolerance on the density 69 logical GRAD !< [Input] .True. if the functional is a GGA 70c 71c Explicit first derivatives of density wrt current nuclei [input] 72c 73 double precision drho(nq,ipol,3,ncuratoms) !< [Input] Derivative 74 !< of rho wrt nuclear coordinates 75 double precision ddelrho(nq,3,ipol,3,ncuratoms) !< [Input] 76 !< Derivative of the density gradient wrt nuclear coordinates 77 double precision dttau(nq,ipol,3,ncuratoms) !< [Input] Derivative 78 !< of the kinetic energy density wrt nuclear coordinates 79c 80c Space for separable coefficients of first derivatives of density 81c 82 double precision Acof2(nq,ipol,3) !< [Scratch] Derivative rho 83 double precision Ccof2(nq,3,ipol,3) !< [Scratch] Derivative grad 84 double precision Mcof2(nq,ipol,3) !< [Scratch] Derivative tau 85c 86c Spin density gradients 87c 88 double precision delrho(nq,3,ipol) !< [Input] Density gradient 89c 90 integer g_dens(ipol) !< [Input] GA handle for density matrices 91c 92c Hessian matrix (updated) 93c 94 double precision hess(3,natoms,3,natoms) !< [In/Output] Hessian 95c 96 double precision scr(nq,15) !< [Scratch] matrix 97 double precision rchi_atom(natoms) !< [Input] Screening parameters 98 double precision rdelchi_atom(natoms) !< [Input] Screening 99 !< parameters 100 double precision rdens_atom(natoms,natoms,ipol) !< [Input] 101 !< Screening parameters 102 integer cetobfr(2,natoms) !< [Input] First and last basis 103 !< function of the atoms 104c 105 double precision Pmat(max_at_bf*max_at_bf) !< [Scratch] vector 106 double precision F(max_at_bf*max_at_bf) !< [Scratch] vector 107 double precision ff(nq,3,*) !< [Scratch] arrays 108 double precision s(nq,max_at_bf) !< [Scratch] arrays 109c 110c Sampling Matrices for the XC Functional Derivatives 111c 112 double precision Amat(nq,ipol) !< [Input] Derivative of functional 113 !< wrt rho 114 double precision Cmat(nq,3,ipol) !< [Input] Derivative of the 115 !< functional wrt the norm of the electron density gradient. 116 !< The call to `transform_Cmat` changes the contents to 117 !< \f$\gamma_{\alpha\alpha}\cdot\frac{\partial \rho_\alpha}{\partial x} 118 !< \ldots \gamma_{\beta\beta}\cdot\frac{\partial \rho_\beta}{\partial z}\f$ 119 double precision Mmat(nq,ipol) !< [Input] Derivative of functional 120 !< wrt kinetic energy density 121c 122 double precision Amat2(nq,NCOL_AMAT2) !< [Input] 2nd derivative 123 !< of functional wrt rho 124 double precision Cmat2(nq,NCOL_CMAT2) !< [Input] 2nd derivative 125 !< of functional wrt gamma 126 double precision Mmat2(nq,NCOL_MMAT2) !< [Input] 2nd derivative 127 !< of functional wrt kinetic energy density 128c 129c Sampling Matrices for [Products of] Basis Functions & Gradients 130c 131 integer iniz(natoms) !< [Input] Start something 132 integer ifin(natoms) !< [Input] End something 133c 134c Basis Functions & Derivatives 135c 136 double precision chi(nq,nbf) !< [Input] Basis function values 137 double precision delchi(nq,3,nbf) !< [Input] Basis function 138 !< derivative values 139 double precision heschi(nq,6,nbf) !< [Input] Basis function 140 !< 2nd derivative values 141 double precision d3chi(nq,10,nbf) !< [Input] Basis function 142 !< 3rd derivative values 143c 144 integer ibf(nbf) !< [Input] Some mapping table 145c 146c local declarations 147c 148 logical oprint 149 double precision A_MAX, C_MAX, AC_MAX, FUNC_MAXI, FUNC_MAXJ 150 double precision FUNC_MAX, DELFUNC_MAX, tol_rho_tmp 151 integer iatcur, jatcur 152 integer iat, inizia, ifinia, nbfia, nnia, ifirst, ilast, idim 153 integer jat, inizja, ifinja, nbfja, nnja, jfirst, jlast, jdim 154 integer ii, mu, nu, icount 155 integer n 156 double precision aaa, fdchix, fdchiy, fdchiz, 157 & ccc1, ccc2, ccc3 158 double precision T(3,3),t6(6) 159 integer idir, jdir 160c 161c The following parameter definitions must be consistent with 162c those in routine xc_eval_basis, or this routine will not work 163c 164 integer iixx,iixy,iixz, 165 & iiyy,iiyz, 166 & iizz 167c 168 parameter ( iixx=1,iixy=2,iixz=3, 169 & iiyy=4,iiyz=5, 170 & iizz=6 ) 171c 172 integer iixxx,iixxy,iixxz, 173 & iixyy,iixyz, 174 & iixzz, 175 & iiyyy,iiyyz, 176 & iiyzz, 177 & iizzz 178c 179 parameter ( iixxx=1,iixxy=2,iixxz=3, 180 & iixyy=4,iixyz=5, 181 & iixzz=6, 182 & iiyyy=7,iiyyz=8, 183 & iiyzz=9, 184 & iizzz=10 ) 185c 186 double precision duefac 187 double precision dabsmax 188 external dabsmax 189c 190c d2Exc / d2fxc dp(i) dp(j) / dfxc d2p(i) 191c ----- = sum sum | ----------- ----- ----- + sum | ----- ------ 192c dx dy i j / dp(i) dp(j) dx dy i / dp(i) dx dy 193c 194c where "p" represents a density parameter in the set 195c { ra, rb, gaa, gab, gbb }, and the i, j indices run over this set 196c 197c First step: separable part of 2nd derivative contribution 198c 199 oprint= util_print('xc_hessian',print_debug) 200 do 10 iat = 1, natoms 201 iatcur = curatoms(iat) 202 if (iatcur.eq.0) goto 10 203c 204c Form Acof2, Ccof2 for iatcur 205c 206 call dcopy(nq*ipol*3,drho(1,1,1,iatcur),1,Acof2,1) 207 if (grad) then 208 call dcopy(nq*ipol*9,ddelrho(1,1,1,1,iatcur),1,Ccof2,1) 209 endif 210 if (kske) then 211 call dcopy(nq*ipol*3,dttau(1,1,1,iatcur),1,Mcof2,1) 212 end if 213c 214 call xc_cpks_coeff(Acof2, Ccof2, Mcof2, 215 & Amat2, Cmat2, Cmat, Mmat2, 216 & delrho,3, ipol, nq, grad, kske, .false.) 217c 218 do 20 jat = 1, iat 219 jatcur = curatoms(jat) 220 if (jatcur.eq.0) goto 20 221c 222 call dfill(9,0.d0,T,1) 223 if (ipol.eq.1) then 224 if (.not.GRAD) then 225 do jdir = 1, 3 226 do idir = 1, 3 227 T(idir,jdir) = T(idir,jdir) + 228 & ddot(nq,Acof2(1,1,idir),1, 229 . drho(1,1,jdir,jatcur),1) 230 enddo 231 enddo 232 else 233 do jdir = 1, 3 234 do idir = 1, 3 235 do n = 1, nq 236 T(idir,jdir) = T(idir,jdir) 237 & + Acof2(n,1,idir)*drho(n,1,jdir,jatcur) 238 & + Ccof2(n,1,1,idir)*ddelrho(n,1,1,jdir,jatcur) 239 & + Ccof2(n,2,1,idir)*ddelrho(n,2,1,jdir,jatcur) 240 & + Ccof2(n,3,1,idir)*ddelrho(n,3,1,jdir,jatcur) 241 enddo 242 enddo 243 enddo 244 endif 245 else 246 do jdir = 1, 3 247 do idir = 1, 3 248 T(idir,jdir) = T(idir,jdir) + 249 & ddot(nq,Acof2(1,1,idir),1, 250 . drho(1,1,jdir,jatcur),1) + 251 & ddot(nq,Acof2(1,2,idir),1, 252 . drho(1,2,jdir,jatcur),1) 253 enddo 254 enddo 255 if (GRAD) then 256 do jdir = 1, 3 257 do idir = 1, 3 258 do n = 1, nq 259 T(idir,jdir) = T(idir,jdir) 260 & + Ccof2(n,1,1,idir)*ddelrho(n,1,1,jdir,jatcur) 261 & + Ccof2(n,2,1,idir)*ddelrho(n,2,1,jdir,jatcur) 262 & + Ccof2(n,3,1,idir)*ddelrho(n,3,1,jdir,jatcur) 263 & + Ccof2(n,1,2,idir)*ddelrho(n,1,2,jdir,jatcur) 264 & + Ccof2(n,2,2,idir)*ddelrho(n,2,2,jdir,jatcur) 265 & + Ccof2(n,3,2,idir)*ddelrho(n,3,2,jdir,jatcur) 266 enddo 267 enddo 268 enddo 269 endif 270 endif 271c 272c Update Hessian block(s) 273c 274 do jdir = 1,3 275 do idir = 1,3 276 hess(idir,iat,jdir,jat) = hess(idir,iat,jdir,jat) 277 & + T(idir,jdir) 278 if (iat.ne.jat) then 279 hess(jdir,jat,idir,iat) = hess(jdir,jat,idir,iat) 280 & + T(idir,jdir) 281 endif 282 enddo 283 enddo 284c 285 20 continue 286 10 continue 287c 288c Second step: remaining terms involving functional first derivatives 289c and density parameter second derivatives 290c 291c We now need Cmat in the delrho form 292c 293 if (GRAD) call transform_Cmat(delrho, Cmat, ipol, nq) 294c 295 A_MAX = dabsmax(nq*ipol,Amat) 296 C_MAX = dabsmax(nq*3*ipol,Cmat) 297 AC_MAX = max(A_MAX,C_MAX) 298c 299#if 0 300 write(6,*) ' xc_d2expl: AMAT ' 301 call output(amat, 1, nq, 1, ipol, nq, ipol, 1) 302 if (GRAD) then 303 write(6,*) ' xc_d2expl: CMAT ' 304 call output(cmat, 1, nq, 1, 3*ipol, nq, 3*ipol, 1) 305 endif 306 write(6,*) ' xc_d2expl: chi ' 307 call output(chi, 1, nq, 1, nbf, nq, nbf, 1) 308 if (GRAD) then 309 write(6,*) ' xc_d2expl: delchi ' 310 call output(delchi, 1, nq, 1, 3*nbf, nq, 3*nbf, 1) 311 endif 312#endif 313c 314c Screening is accomplished by: p(r) <= |Xi(r)|*|Xj(r)|*|Dij| 315c Xi(r) is screened on desired accuracy/max(|Xj(r)|)*max(|Dij|) 316c Dij is screened on desired accuracy/max(|Xi(r)|)*max(|Xj(r)|) 317c 318 FUNC_MAX = dabsmax(natoms,rchi_atom) 319 DELFUNC_MAX = dabsmax(natoms,rdelchi_atom) 320c 321 do 230 iat = 1, natoms 322 inizia = iniz(iat) 323 if (inizia.eq.0)goto 230 324 iatcur = curatoms(iat) 325 ifinia = ifin(iat) 326 ifirst = cetobfr(1,iat) 327 ilast = cetobfr(2,iat) 328 nbfia = ilast-ifirst+1 329 nnia = ifinia-inizia+1 330c 331c screening parameters 332c 333 FUNC_MAXI = max(rchi_atom(iat),rdelchi_atom(iat)) 334 FUNC_MAXJ = max(FUNC_MAX,DELFUNC_MAX) 335#if 0 336 if (ipol.gt.1)then 337 P_MAXJ_A = dabsmax(natoms,rdens_atom(1,iat,1)) 338 P_MAXJ_B = dabsmax(natoms,rdens_atom(1,iat,2)) 339 P_MAXJ = MAX(P_MAXJ_A, P_MAXJ_B) 340 else 341 P_MAXJ = dabsmax(natoms,rdens_atom(1,iat,1)) 342 endif 343 if (FUNC_MAXI*FUNC_MAXJ*P_MAXJ.lt.tol_rho) goto 225 344c !!! Cutoff temporarily commented out !!! 345#endif 346 do 220 jat = 1, iat 347 inizja = iniz(jat) 348 if (inizja.eq.0)goto 220 349 jatcur = curatoms(jat) 350 if (iatcur .eq. 0 .and. jatcur .eq. 0) goto 220 351 ifinja = ifin(jat) 352 jfirst = cetobfr(1,jat) 353 jlast = cetobfr(2,jat) 354 nbfja = jlast-jfirst+1 355 nnja = ifinja-inizja+1 356c 357c screening parameters 358c 359 FUNC_MAXJ = max(rchi_atom(jat),rdelchi_atom(jat)) 360#if 0 361 if (ipol.eq.1)then 362 P_MAXIJ = rdens_atom(iat,jat,1) 363 else 364 P_MAXIJ = max(rdens_atom(iat,jat,1), 365 & rdens_atom(iat,jat,2)) 366 endif 367 if (FUNC_MAXI*FUNC_MAXJ*P_MAXIJ.lt.tol_rho) goto 215 368c !!! Cutoff temporarily commented out !!! 369#endif 370 tol_rho_tmp = tol_rho/(FUNC_MAXI*FUNC_MAXJ) 371c 372 do 210 ii = 1, ipol 373c 374c screening parameters 375c 376#if 0 377 P_MAXIJ = rdens_atom(iat,jat,ii) 378 if (FUNC_MAXI*FUNC_MAXJ*P_MAXIJ.lt.tol_rho)goto 210 379c !!! Cutoff temporarily commented out !!! 380#endif 381c 382 call get_atom_block(g_dens(ii), basis, 383 & iat, jat, Pmat, idim, jdim) 384c 385 call gat_mat(F, Pmat, nbfia, nbfja, nnia, nnja,ifirst, 386 & jfirst, ibf(inizia), ibf(inizja)) 387c 388c Three terms to compute 389c 390c First term: Xiat(r)*hessXjat(r)*Diat,jat -> hess(jat,jat) 391c GC: delXiat(r)*hessXjat(r)*Diat,jat -> hess(jat,jat) 392c GC: Xiat(r)* d3Xjat(r)*Diat,jat -> hess(jat,jat) 393c 394 if (jatcur .ne. 0) then 395 call dfill(9,0.d0,T,1) 396 call dgemm('n','n',nq,nnja,nnia,1d0, 397 A chi(1,inizia),nq,F,nnia,0d0,s,nq) 398 if (grad) 399 G call dgemm('n','n',nq*3,nnja,nnia,1d0, 400 A delchi(1,1,inizia),nq*3,F,nnia,0d0,ff,nq*3) 401 t6(1)=0d0 402 t6(2)=0d0 403 t6(3)=0d0 404 t6(4)=0d0 405 t6(5)=0d0 406 t6(6)=0d0 407 do mu=inizja,ifinja 408 if (GRAD) then 409 do n = 1, nq 410 ff(n,1,1) = 411 S Amat(n,ii)*s(n,mu-inizja+1) 412 & + Cmat(n,1,ii)*ff(n,1,mu-inizja+1) 413 & + Cmat(n,2,ii)*ff(n,2,mu-inizja+1) 414 & + Cmat(n,3,ii)*ff(n,3,mu-inizja+1) 415 enddo 416 else 417 do n = 1, nq 418 ff(n,1,1) = Amat(n,ii)* 419 * s(n,mu-inizja+1) 420 enddo 421 endif 422 call dgemv('t',nq, 6, 1d0,heschi(1,1,mu),nq, 423 Y ff(1,1,1),1,1d0,t6,1) 424 enddo 425 T(1,1) = T(1,1)+t6(1) 426 T(1,2) = T(1,2)+t6(2) 427 T(1,3) = T(1,3)+t6(3) 428 T(2,2) = T(2,2)+t6(4) 429 T(2,3) = T(2,3)+t6(5) 430 T(3,3) = T(3,3)+t6(6) 431 if (GRAD) then 432 do mu = inizja, ifinja 433 do n = 1, nq 434 ccc1 = Cmat(n,1,ii)*s(n,mu-inizja+1) 435 ccc2 = Cmat(n,2,ii)*s(n,mu-inizja+1) 436 ccc3 = Cmat(n,3,ii)*s(n,mu-inizja+1) 437 T(1,1) = T(1,1) + d3chi(n,iixxx,mu)*ccc1 438 & + d3chi(n,iixxy,mu)*ccc2 439 & + d3chi(n,iixxz,mu)*ccc3 440 T(1,2) = T(1,2) + d3chi(n,iixxy,mu)*ccc1 441 & + d3chi(n,iixyy,mu)*ccc2 442 & + d3chi(n,iixyz,mu)*ccc3 443 T(1,3) = T(1,3) + d3chi(n,iixxz,mu)*ccc1 444 & + d3chi(n,iixyz,mu)*ccc2 445 & + d3chi(n,iixzz,mu)*ccc3 446 T(2,2) = T(2,2) + d3chi(n,iixyy,mu)*ccc1 447 & + d3chi(n,iiyyy,mu)*ccc2 448 & + d3chi(n,iiyyz,mu)*ccc3 449 T(2,3) = T(2,3) + d3chi(n,iixyz,mu)*ccc1 450 & + d3chi(n,iiyyz,mu)*ccc2 451 & + d3chi(n,iiyzz,mu)*ccc3 452 T(3,3) = T(3,3) + d3chi(n,iixzz,mu)*ccc1 453 & + d3chi(n,iiyzz,mu)*ccc2 454 & + d3chi(n,iizzz,mu)*ccc3 455 enddo 456 enddo 457 endif 458c 459 duefac=1d0 460 if (iat.ne.jat) duefac=2d0 461c 462 T(2,1) = T(1,2) 463 T(3,1) = T(1,3) 464 T(3,2) = T(2,3) 465 do jdir = 1,3 466 do idir = 1,3 467 hess(idir,jat,jdir,jat) = 468 & hess(idir,jat,jdir,jat) + T(idir,jdir)*duefac 469 enddo 470 enddo 471 endif 472c 473c Second term: hessXiat(r)* Xjat(r)*Diat,jat -> hess(iat,iat) 474c GC: hessXiat(r)*delXjat(r)*Diat,jat -> hess(iat,iat) 475c GC: d3Xiat(r)* Xjat(r)*Diat,jat -> hess(iat,iat) 476c 477 if (iatcur .ne. 0) then 478 call dfill(9,0.d0,T,1) 479 call dgemm('n','t',nq,nnia,nnja,1d0, 480 A chi(1,inizja),nq,F,nnia,0d0,s,nq) 481 if (grad) 482 G call dgemm('n','t',nq*3,nnia,nnja,1d0, 483 A delchi(1,1,inizja),nq*3,F,nnia,0d0,ff,nq*3) 484 t6(1)=0d0 485 t6(2)=0d0 486 t6(3)=0d0 487 t6(4)=0d0 488 t6(5)=0d0 489 t6(6)=0d0 490 do nu=inizia,ifinia 491 if (GRAD) then 492 do n = 1, nq 493 ff(n,1,1) = 494 S Amat(n,ii)*s(n,nu-inizia+1) 495 & + Cmat(n,1,ii)*ff(n,1,nu-inizia+1) 496 & + Cmat(n,2,ii)*ff(n,2,nu-inizia+1) 497 & + Cmat(n,3,ii)*ff(n,3,nu-inizia+1) 498 enddo 499 else 500 do n = 1, nq 501 ff(n,1,1) = 502 A Amat(n,ii)*s(n,nu-inizia+1) 503 enddo 504 endif 505 call dgemv('t',nq, 6, 1d0,heschi(1,1,nu),nq, 506 Y ff(1,1,1),1,1d0,t6,1) 507 enddo 508 T(1,1) = T(1,1)+t6(1) 509 T(1,2) = T(1,2)+t6(2) 510 T(1,3) = T(1,3)+t6(3) 511 T(2,2) = T(2,2)+t6(4) 512 T(2,3) = T(2,3)+t6(5) 513 T(3,3) = T(3,3)+t6(6) 514 if (GRAD) then 515 do nu = inizia, ifinia 516c 517 do n = 1, nq 518 ccc1 = Cmat(n,1,ii)*s(n,nu-inizia+1) 519 ccc2 = Cmat(n,2,ii)*s(n,nu-inizia+1) 520 ccc3 = Cmat(n,3,ii)*s(n,nu-inizia+1) 521 T(1,1) = T(1,1) + d3chi(n,iixxx,nu)*ccc1 522 & + d3chi(n,iixxy,nu)*ccc2 523 & + d3chi(n,iixxz,nu)*ccc3 524 T(1,2) = T(1,2) + d3chi(n,iixxy,nu)*ccc1 525 & + d3chi(n,iixyy,nu)*ccc2 526 & + d3chi(n,iixyz,nu)*ccc3 527 T(1,3) = T(1,3) + d3chi(n,iixxz,nu)*ccc1 528 & + d3chi(n,iixyz,nu)*ccc2 529 & + d3chi(n,iixzz,nu)*ccc3 530 T(2,2) = T(2,2) + d3chi(n,iixyy,nu)*ccc1 531 & + d3chi(n,iiyyy,nu)*ccc2 532 & + d3chi(n,iiyyz,nu)*ccc3 533 T(2,3) = T(2,3) + d3chi(n,iixyz,nu)*ccc1 534 & + d3chi(n,iiyyz,nu)*ccc2 535 & + d3chi(n,iiyzz,nu)*ccc3 536 T(3,3) = T(3,3) + d3chi(n,iixzz,nu)*ccc1 537 & + d3chi(n,iiyzz,nu)*ccc2 538 & + d3chi(n,iizzz,nu)*ccc3 539 enddo 540 enddo 541 endif 542c 543 duefac=1d0 544 if (iat.ne.jat) duefac=2d0 545c 546 T(2,1) = T(1,2) 547 T(3,1) = T(1,3) 548 T(3,2) = T(2,3) 549 do jdir = 1,3 550 do idir = 1,3 551 hess(idir,iat,jdir,iat) = 552 & hess(idir,iat,jdir,iat) + T(idir,jdir)*duefac 553 enddo 554 enddo 555 endif 556c 557c Third term: delXiat(r)*del(T)Xjat(r)*Diat,jat -> hess(iat,jat) 558c 559 if (jatcur .ne. 0 .and. iatcur .ne. 0) then 560 call dfill(9,0.d0,T,1) 561 icount = 0 562 call dgemm('n','n',nq*3,nnja,nnia,1d0, 563 A delchi(1,1,inizia),nq*3,F,nnia,0d0,ff,nq*3) 564 do mu = inizja, ifinja 565 do n = 1, nq 566 fdchix = Amat(n,ii)*ff(n,1,mu-inizja+1) 567 fdchiy = Amat(n,ii)*ff(n,2,mu-inizja+1) 568 fdchiz = Amat(n,ii)*ff(n,3,mu-inizja+1) 569 T(1,1) = T(1,1) + fdchix*delchi(n,1,mu) 570 T(1,2) = T(1,2) + fdchix*delchi(n,2,mu) 571 T(1,3) = T(1,3) + fdchix*delchi(n,3,mu) 572 T(2,1) = T(2,1) + fdchiy*delchi(n,1,mu) 573 T(2,2) = T(2,2) + fdchiy*delchi(n,2,mu) 574 T(2,3) = T(2,3) + fdchiy*delchi(n,3,mu) 575 T(3,1) = T(3,1) + fdchiz*delchi(n,1,mu) 576 T(3,2) = T(3,2) + fdchiz*delchi(n,2,mu) 577 T(3,3) = T(3,3) + fdchiz*delchi(n,3,mu) 578 enddo 579 enddo 580 if (GRAD) then 581 do mu = inizja, ifinja 582 call dfill(nq*6,0.d0,scr,1) 583 do nu = inizia, ifinia 584 icount = icount+1 585 aaa = F(icount) 586 if (abs(aaa).gt.tol_rho_tmp)then 587 call daxpy(nq*6,aaa,heschi(1,1,nu),1, 588 & scr,1) 589 endif 590 enddo 591 do n = 1,nq 592 s(n,1) = Cmat(n,1,ii)*scr(n,iixx) 593 & + Cmat(n,2,ii)*scr(n,iixy) 594 & + Cmat(n,3,ii)*scr(n,iixz) 595 s(n,2) = Cmat(n,1,ii)*scr(n,iixy) 596 & + Cmat(n,2,ii)*scr(n,iiyy) 597 & + Cmat(n,3,ii)*scr(n,iiyz) 598 s(n,3) = Cmat(n,1,ii)*scr(n,iixz) 599 & + Cmat(n,2,ii)*scr(n,iiyz) 600 & + Cmat(n,3,ii)*scr(n,iizz) 601 enddo 602 call dgemm('t','n',3,3,nq,1d0, 603 A s,nq,delchi(1,1,mu),nq,1d0,T,3) 604 do n = 1, nq 605 s(n,1) = Cmat(n,1,ii)*heschi(n,iixx,mu) 606 & + Cmat(n,2,ii)*heschi(n,iixy,mu) 607 & + Cmat(n,3,ii)*heschi(n,iixz,mu) 608 s(n,2) = Cmat(n,1,ii)*heschi(n,iixy,mu) 609 & + Cmat(n,2,ii)*heschi(n,iiyy,mu) 610 & + Cmat(n,3,ii)*heschi(n,iiyz,mu) 611 s(n,3) = Cmat(n,1,ii)*heschi(n,iixz,mu) 612 & + Cmat(n,2,ii)*heschi(n,iiyz,mu) 613 & + Cmat(n,3,ii)*heschi(n,iizz,mu) 614 enddo 615 call dgemm('t','n',3,3,nq,1d0, 616 A ff(1,1,mu-inizja+1),nq,s,nq,1d0,T,3) 617 enddo 618 endif 619c 620c This term always comes with a factor of 2 in front 621c 622 duefac=2d0 623c 624 do jdir = 1,3 625 do idir = 1,3 626 hess(idir,iat,jdir,jat) = 627 & hess(idir,iat,jdir,jat) + T(idir,jdir)*duefac 628 if (iat.ne.jat) then 629 hess(jdir,jat,idir,iat) = 630 & hess(jdir,jat,idir,iat) + T(idir,jdir)*duefac 631 endif 632 enddo 633 enddo 634 endif 635 210 continue 636 220 continue 637 230 continue 638 if(oprint.and.ga_nodeid().eq.0) then 639 write(luout,*) ' xc_d2expl: hess ' 640 call output(hess,1,3*natoms,1,3*natoms,3*natoms,3*natoms,1) 641 endif 642c 643 return 644 end 645C> @} 646