1* 2* $Id$ 3* 4 5 6 7* ************************************ 8* * * 9* * v_bwexc_euv * 10* * * 11* ************************************ 12* 13 14 subroutine v_bwexc_euv(gga,n2ft3d,ispin,dn, 15 > x_parameter,c_parameter, 16 > stress) 17 implicit none 18 integer gga 19 integer n2ft3d 20 integer ispin 21 real*8 dn(n2ft3d,2) 22 real*8 x_parameter,c_parameter 23 real*8 stress(3,3) 24 25 26#include "bafdecls.fh" 27#include "nwpwxc.fh" 28#include "util.fh" 29 30 31c **** local variables **** 32 logical value, use_nwpwxc 33 integer nx,ny,nz,u,v,s 34 real*8 scal1,pi,scal,sum 35 integer rho(2),grx(2),gry(2),grz(2) 36 integer agr(3),fn(2),fdn(3),tmp(2),rhog(2),xce(2) 37 38 integer rhoup(2),grupx(2),grupy(2),grupz(2) 39 integer rhodn(2),grdnx(2),grdny(2),grdnz(2) 40 integer grallx(2),grally(2),grallz(2) 41 integer grad(2),gtmp(2) 42 integer xagr(2),xfn(2),xfdn(2) 43 real*8 hm(3,3),W(3,3),omega 44 real*8 dncut,dumtau 45 parameter(dncut = 1.0d-30) 46 47* **** external functions **** 48 integer G_indx 49 real*8 lattice_unitg,lattice_omega 50 external G_indx 51 external lattice_unitg,lattice_omega 52 53 call nwpw_timing_start(4) 54 call D3dB_nx(1,nx) 55 call D3dB_ny(1,ny) 56 call D3dB_nz(1,nz) 57 scal1 = 1.0d0/dble(nx*ny*nz) 58 omega = lattice_omega() 59 60 61* *** define hm **** 62 pi = 4.0d0*datan(1.0d0) 63 scal = 1.0d0/(2.0d0*pi) 64 do v=1,3 65 do u=1,3 66 hm(u,v) = scal*lattice_unitg(u,v) 67 end do 68 end do 69 70 use_nwpwxc = .false. 71 use_nwpwxc = nwpwxc_is_on() 72 73 74* ********************************** 75* ***** restricted calculation ***** 76* ********************************** 77 if (ispin.eq.1) then 78 79c ***** tempory variables needed rho,grx,gry,grz ***** 80c ***** agr,fn,fdn ***** 81 value = BA_push_get(mt_dbl,n2ft3d,'rho', rho(2), rho(1)) 82 value = value.and. 83 > BA_push_get(mt_dbl, n2ft3d,'grx',grx(2),grx(1)) 84 value = value.and. 85 > BA_push_get(mt_dbl, n2ft3d,'gry',gry(2),gry(1)) 86 value = value.and. 87 > BA_push_get(mt_dbl, n2ft3d,'grz',grz(2),grz(1)) 88 value = value.and. 89 > BA_push_get(mt_dbl, n2ft3d,'agr',agr(2),agr(1)) 90 value = value.and. 91 > BA_push_get(mt_dbl, n2ft3d,'fn',fn(2),fn(1)) 92 value = value.and. 93 > BA_push_get(mt_dbl, 2*n2ft3d,'fdn',fdn(2),fdn(1)) 94 tmp(1) = fdn(1) 95 value = value.and. 96 > BA_push_get(mt_dbl, n2ft3d,'rhog',rhog(2),rhog(1)) 97 value = value.and. 98 > BA_push_get(mt_dbl, n2ft3d,'xce',xce(2),xce(1)) 99 if (.not. value) call errquit('out of stack memory',0,0) 100 call dcopy(n2ft3d,0.0d0,0,dbl_mb(rho(1)),1) 101 call dcopy(n2ft3d,0.0d0,0,dbl_mb(agr(1)),1) 102 103 104 105c ***** calculate rho tmp1=rho(g) **** 106 call D3dB_rr_Sum(1,dn(1,1),dn(1,1),dbl_mb(rho(1))) 107 call D3dB_r_Zero_Ends(1,dbl_mb(rho(1))) 108 call D3dB_r_SMul(1,scal1,dbl_mb(rho(1)),dbl_mb(rhog(1))) 109 call D3dB_rc_fft3f(1,dbl_mb(rhog(1))) 110 call mask_C(0,dbl_mb(rhog(1))) 111 112c ***** calculated grup= grad n **** 113 call D3dB_ic_Mul(1,dbl_mb(G_indx(1)), 114 > dbl_mb(rhog(1)), 115 > dbl_mb(grx(1))) 116 call D3dB_ic_Mul(1,dbl_mb(G_indx(2)), 117 > dbl_mb(rhog(1)), 118 > dbl_mb(gry(1))) 119 call D3dB_ic_Mul(1,dbl_mb(G_indx(3)), 120 > dbl_mb(rhog(1)), 121 > dbl_mb(grz(1))) 122 call D3dB_cr_fft3b(1,dbl_mb(grx(1))) 123 call D3dB_cr_fft3b(1,dbl_mb(gry(1))) 124 call D3dB_cr_fft3b(1,dbl_mb(grz(1))) 125 126c ***** calculate agr = |grad n| **** 127 call D3dB_rr_Sqr(1,dbl_mb(grx(1)), 128 > dbl_mb(agr(1))) 129 call D3dB_rr_Sqr(1,dbl_mb(gry(1)), 130 > dbl_mb(tmp(1))) 131 132c call D3dB_rr_Sum(1,dbl_mb(tmp(1)), 133c > dbl_mb(agr(1)), 134c > dbl_mb(agr(1))) 135 call D3dB_rr_Sum2(1,dbl_mb(tmp(1)),dbl_mb(agr(1))) 136 137 call D3dB_rr_Sqr(1,dbl_mb(grz(1)), 138 > dbl_mb(tmp(1))) 139 140c call D3dB_rr_Sum(1,dbl_mb(tmp(1)), 141c > dbl_mb(agr(1)), 142c > dbl_mb(agr(1))) 143 call D3dB_rr_Sum2(1,dbl_mb(tmp(1)),dbl_mb(agr(1))) 144 145 if (use_nwpwxc) then 146c 147c Evaluate the functional 148c 149 call nwpwxc_eval_df(1,n2ft3d,dbl_mb(rho(1)),dbl_mb(agr(1)), 150 > dumtau,xce, 151 > dbl_mb(fn(1)),dbl_mb(fdn(1)),dumtau) 152cc 153cc Combine (df/d|grad a|) with (df/d(grad a|grad b)) 154cc 155 call D3dB_rr_daxpy(1,0.5d0,dbl_mb(fdn(1)+n2ft3d), 156 > dbl_mb(fdn(1))) 157cc 158cc Calculate the energy density from the functional 159c 160 do u = 1, n2ft3d 161 xce(u) = xce(u)/(dbl_mb(rho(1)-1+u)+dncut) 162 enddo 163c 164c Tranform the gradient terms from (df/d|grad n|^2) to 165c (df/d|grad n|) 166c 167 call D3dB_rr_Sqrt1(1,dbl_mb(agr(1))) 168 call D3dB_rr_Mul2(1,dbl_mb(agr(1)),dbl_mb(fdn(1))) 169c 170 else 171c call D3dB_rr_Sqrt(1,dbl_mb(agr(1)), 172c > dbl_mb(agr(1))) 173 call D3dB_rr_Sqrt1(1,dbl_mb(agr(1))) 174 175 176 177c ***** calculate fdn=df/d|grad n| **** 178 if (gga.eq.10) then 179 call gen_PBE96_BW_restricted(n2ft3d, 180 > dbl_mb(rho(1)), 181 > dbl_mb(agr(1)), 182 > x_parameter,c_parameter, 183 > dbl_mb(xce(1)), !*** not used ??*** 184 > dbl_mb(fn(1)), !*** not used ??*** 185 > dbl_mb(fdn(1))) 186 187 else if (gga.eq.11) then 188 call gen_BLYP_BW_restricted(n2ft3d, 189 > dbl_mb(rho(1)), 190 > dbl_mb(agr(1)), 191 > x_parameter,c_parameter, 192 > dbl_mb(xce(1)), !*** not used ??*** 193 > dbl_mb(fn(1)), !*** not used ??*** 194 > dbl_mb(fdn(1))) 195 else if (gga.eq.12) then 196 call gen_revPBE_BW_restricted(n2ft3d, 197 > dbl_mb(rho(1)), 198 > dbl_mb(agr(1)), 199 > x_parameter,c_parameter, 200 > dbl_mb(xce(1)), !*** not used ??*** 201 > dbl_mb(fn(1)), !*** not used ??*** 202 > dbl_mb(fdn(1))) 203 else if (gga.eq.13) then 204 call gen_PBEsol_BW_restricted(n2ft3d, 205 > dbl_mb(rho(1)), 206 > dbl_mb(agr(1)), 207 > x_parameter,c_parameter, 208 > dbl_mb(xce(1)), !*** not used ??*** 209 > dbl_mb(fn(1)), !*** not used ??*** 210 > dbl_mb(fdn(1))) 211 else if (gga.eq.14) then 212 call gen_HSE_BW_restricted(n2ft3d, 213 > dbl_mb(rho(1)), 214 > dbl_mb(agr(1)), 215 > x_parameter,c_parameter, 216 > dbl_mb(xce(1)), !*** not used ??*** 217 > dbl_mb(fn(1)), !*** not used ??*** 218 > dbl_mb(fdn(1))) 219 else if (gga.eq.15) then 220 call gen_B3LYP_BW_restricted(n2ft3d, 221 > dbl_mb(rho(1)), 222 > dbl_mb(agr(1)), 223 > x_parameter,c_parameter, 224 > dbl_mb(xce(1)), !*** not used ??*** 225 > dbl_mb(fn(1)), !*** not used ??*** 226 > dbl_mb(fdn(1))) 227 else if (gga.eq.16) then 228 call gen_BEEF_BW_restricted(n2ft3d, 229 > dbl_mb(rho(1)), 230 > dbl_mb(agr(1)), 231 > x_parameter,c_parameter,0.6001664769d0, 232 > dbl_mb(xce(1)), !*** not used ??*** 233 > dbl_mb(fn(1)), !*** not used ??*** 234 > dbl_mb(fdn(1))) 235 else if (gga.eq.17) then 236 call gen_BEEF_BW_restricted(n2ft3d, 237 > dbl_mb(rho(1)), 238 > dbl_mb(agr(1)), 239 > x_parameter,c_parameter,0.0d0, 240 > dbl_mb(xce(1)), !*** not used ??*** 241 > dbl_mb(fn(1)), !*** not used ??*** 242 > dbl_mb(fdn(1))) 243 else 244 call errquit('bad gga',0,0) 245 end if 246 endif 247 248 249 250c ***** calculate df/d|grad n| *(grad n)/|grad n| **** 251c call D3dB_rr_Divide(1,dbl_mb(grx(1)), 252c > dbl_mb(agr(1)), 253c > dbl_mb(grx(1))) 254c call D3dB_rr_Divide(1,dbl_mb(gry(1)), 255c > dbl_mb(agr(1)), 256c > dbl_mb(gry(1))) 257c call D3dB_rr_Divide(1,dbl_mb(grz(1)), 258c > dbl_mb(agr(1)), 259c > dbl_mb(grz(1))) 260 call D3dB_rr_Divide2(1,dbl_mb(agr(1)),dbl_mb(grx(1))) 261 call D3dB_rr_Divide2(1,dbl_mb(agr(1)),dbl_mb(gry(1))) 262 call D3dB_rr_Divide2(1,dbl_mb(agr(1)),dbl_mb(grz(1))) 263 264c call D3dB_rr_Mul(1,dbl_mb(grx(1)), 265c > dbl_mb(fdn(1)), 266c > dbl_mb(grx(1))) 267c call D3dB_rr_Mul(1,dbl_mb(gry(1)), 268c > dbl_mb(fdn(1)), 269c > dbl_mb(gry(1))) 270c call D3dB_rr_Mul(1,dbl_mb(grz(1)), 271c > dbl_mb(fdn(1)), 272c > dbl_mb(grz(1))) 273 call D3dB_rr_Mul2(1,dbl_mb(fdn(1)),dbl_mb(grx(1))) 274 call D3dB_rr_Mul2(1,dbl_mb(fdn(1)),dbl_mb(gry(1))) 275 call D3dB_rr_Mul2(1,dbl_mb(fdn(1)),dbl_mb(grz(1))) 276 277c call D3dB_r_SMul(1,scal1,dbl_mb(grx(1)), 278c > dbl_mb(grx(1))) 279c call D3dB_r_SMul(1,scal1,dbl_mb(gry(1)), 280c > dbl_mb(gry(1))) 281c call D3dB_r_SMul(1,scal1,dbl_mb(grz(1)), 282c > dbl_mb(grz(1))) 283 call D3dB_r_SMul1(1,scal1,dbl_mb(grx(1))) 284 call D3dB_r_SMul1(1,scal1,dbl_mb(gry(1))) 285 call D3dB_r_SMul1(1,scal1,dbl_mb(grz(1))) 286 287 call D3dB_r_Zero_Ends(1,dbl_mb(grx(1))) 288 call D3dB_r_Zero_Ends(1,dbl_mb(gry(1))) 289 call D3dB_r_Zero_Ends(1,dbl_mb(grz(1))) 290 call D3dB_rc_fft3f(1,dbl_mb(grx(1))) 291 call D3dB_rc_fft3f(1,dbl_mb(gry(1))) 292 call D3dB_rc_fft3f(1,dbl_mb(grz(1))) 293 294 295* **** W(u,s) = Sum(G) [i*G(u)*dcongj(rhog)*gr(s)] **** 296* **** where gr(1)=grx,gr(2)=gry,gr(3)=grz **** 297 298 !call Pack_c_pack(0,dbl_mb(rhog(1))) 299 call mask_C(0,dbl_mb(rhog(1))) 300 do u=1,3 301 302* **** agr = i*G(u)*grx **** 303 call D3dB_ic_Mul(1,dbl_mb(G_indx(u)), 304 > dbl_mb(grx(1)), 305 > dbl_mb(agr(1))) 306 call mask_C(0,dbl_mb(agr(1))) 307 call D3dB_cc_dot(1,dbl_mb(rhog(1)),dbl_mb(agr(1)),sum) 308 !call Pack_c_pack(0,dbl_mb(agr(1))) 309 !call Pack_cc_dot(0,dbl_mb(rhog(1)),dbl_mb(agr(1)),sum) 310 W(u,1) = sum*omega 311 312* **** agr = i*G(u)*gry **** 313 call D3dB_ic_Mul(1,dbl_mb(G_indx(u)), 314 > dbl_mb(gry(1)), 315 > dbl_mb(agr(1))) 316 call mask_C(0,dbl_mb(agr(1))) 317 call D3dB_cc_dot(1,dbl_mb(rhog(1)),dbl_mb(agr(1)),sum) 318 !call Pack_c_pack(0,dbl_mb(agr(1))) 319 !call Pack_cc_dot(0,dbl_mb(rhog(1)),dbl_mb(agr(1)),sum) 320 W(u,2) = sum*omega 321 322* **** agr = i*G(u)*grz **** 323 call D3dB_ic_Mul(1,dbl_mb(G_indx(u)), 324 > dbl_mb(grz(1)), 325 > dbl_mb(agr(1))) 326 call mask_C(0,dbl_mb(agr(1))) 327 call D3dB_cc_dot(1,dbl_mb(rhog(1)),dbl_mb(agr(1)),sum) 328 !call Pack_c_pack(0,dbl_mb(agr(1))) 329 !call Pack_cc_dot(0,dbl_mb(rhog(1)),dbl_mb(agr(1)),sum) 330 W(u,3) = sum*omega 331 end do 332 333 334 335 336* **** deallocate temporary memory **** 337 value = BA_pop_stack(xce(2)) 338 value = value.and.BA_pop_stack(rhog(2)) 339 value = value.and.BA_pop_stack(fdn(2)) 340 value = value.and.BA_pop_stack(fn(2)) 341 value = value.and.BA_pop_stack(agr(2)) 342 value = value.and.BA_pop_stack(grz(2)) 343 value = value.and.BA_pop_stack(gry(2)) 344 value = value.and.BA_pop_stack(grx(2)) 345 value = value.and.BA_pop_stack(rho(2)) 346 if (.not. value) call errquit('cannot pop stack memory',0,0) 347 348 349 350* ******************************************************* 351* ***** unrestricted calculation ***** 352* ******************************************************* 353 else 354 355c ***** tempory variables needed rho,grx,gry,grz ***** 356c ***** agr,fn,fdn ***** 357 value = BA_push_get(mt_dbl,n2ft3d,'rhoup', rhoup(2), rhoup(1)) 358 value = value.and. 359 > BA_push_get(mt_dbl, n2ft3d,'grupx',grupx(2),grupx(1)) 360 value = value.and. 361 > BA_push_get(mt_dbl, n2ft3d,'grupy',grupy(2),grupy(1)) 362 value = value.and. 363 > BA_push_get(mt_dbl, n2ft3d,'grupz',grupz(2),grupz(1)) 364 365 value = value.and. 366 > BA_push_get(mt_dbl,n2ft3d,'rhodn', rhodn(2), rhodn(1)) 367 value = value.and. 368 > BA_push_get(mt_dbl, n2ft3d,'grdnx',grdnx(2),grdnx(1)) 369 value = value.and. 370 > BA_push_get(mt_dbl, n2ft3d,'grdny',grdny(2),grdny(1)) 371 value = value.and. 372 > BA_push_get(mt_dbl, n2ft3d,'grdnz',grdnz(2),grdnz(1)) 373 374 value = value.and. 375 > BA_push_get(mt_dbl, n2ft3d,'grallx',grallx(2),grallx(1)) 376 value = value.and. 377 > BA_push_get(mt_dbl, n2ft3d,'grally',grally(2),grally(1)) 378 value = value.and. 379 > BA_push_get(mt_dbl, n2ft3d,'grallz',grallz(2),grallz(1)) 380 381 value = value.and. 382 > BA_push_get(mt_dbl, 3*n2ft3d,'xagr',xagr(2),xagr(1)) 383 agr(1) = xagr(1) 384 agr(2) = xagr(1) + n2ft3d 385 agr(3) = xagr(1) + 2*n2ft3d 386 value = value.and. 387 > BA_push_get(mt_dbl, 3*n2ft3d,'grad',grad(2),grad(1)) 388 value = value.and. 389 > BA_push_get(mt_dbl, 3*n2ft3d,'gtmp',gtmp(2),gtmp(1)) 390 value = value.and. 391 > BA_push_get(mt_dbl, 2*n2ft3d,'xfn',xfn(2),xfn(1)) 392 fn(1) = xfn(1) 393 fn(2) = xfn(1)+n2ft3d 394 value = value.and. 395 > BA_push_get(mt_dbl, 3*n2ft3d,'xfdn',xfdn(2),xfdn(1)) 396 fdn(1) = xfdn(1) 397 fdn(2) = xfdn(1) + n2ft3d 398 fdn(3) = xfdn(1) + 2*n2ft3d 399 tmp(1) = xfdn(1) 400 value = value.and. 401 > BA_push_get(mt_dbl, n2ft3d,'xce',xce(2),xce(1)) 402 if (.not. value) call errquit('out of stack memory',0,0) 403 call dcopy(n2ft3d,0.0d0,0,dbl_mb(rhoup(1)),1) 404 call dcopy(n2ft3d,0.0d0,0,dbl_mb(rhodn(1)),1) 405 call dcopy(3*n2ft3d,0.0d0,0,dbl_mb(xagr(1)),1) 406 407 408c ***** calculate rhoup **** 409 call D3dB_r_Copy(1,dn(1,1),dbl_mb(rhoup(1))) 410 call D3dB_r_Zero_Ends(1,dbl_mb(rhoup(1))) 411c call D3dB_r_SMul(1,scal1,dbl_mb(rhoup(1)), 412c > dbl_mb(rhoup(1))) 413 call D3dB_r_SMul1(1,scal1,dbl_mb(rhoup(1))) 414 call D3dB_rc_fft3f(1,dbl_mb(rhoup(1))) 415 call mask_C(0,dbl_mb(rhoup(1))) 416 417c ***** calculate grup= grad nup **** 418 call D3dB_ic_Mul(1,dbl_mb(G_indx(1)), 419 > dbl_mb(rhoup(1)), 420 > dbl_mb(grupx(1))) 421 call D3dB_ic_Mul(1,dbl_mb(G_indx(2)), 422 > dbl_mb(rhoup(1)), 423 > dbl_mb(grupy(1))) 424 call D3dB_ic_Mul(1,dbl_mb(G_indx(3)), 425 > dbl_mb(rhoup(1)), 426 > dbl_mb(grupz(1))) 427 call D3dB_cr_fft3b(1,dbl_mb(grupx(1))) 428 call D3dB_cr_fft3b(1,dbl_mb(grupy(1))) 429 call D3dB_cr_fft3b(1,dbl_mb(grupz(1))) 430 431c ***** calculate agrup = |grad nup| **** 432 call D3dB_rr_Sqr(1,dbl_mb(grupx(1)), 433 > dbl_mb(agr(1))) 434 call D3dB_rr_Sqr(1,dbl_mb(grupy(1)), 435 > dbl_mb(tmp(1))) 436 437c call D3dB_rr_Sum(1,dbl_mb(tmp(1)), 438c > dbl_mb(agr(1)), 439c > dbl_mb(agr(1))) 440 call D3dB_rr_Sum2(1,dbl_mb(tmp(1)),dbl_mb(agr(1))) 441 442 call D3dB_rr_Sqr(1,dbl_mb(grupz(1)), 443 > dbl_mb(tmp(1))) 444 445c call D3dB_rr_Sum(1,dbl_mb(tmp(1)), 446c > dbl_mb(agr(1)), 447c > dbl_mb(agr(1))) 448 call D3dB_rr_Sum2(1,dbl_mb(tmp(1)),dbl_mb(agr(1))) 449 450c call D3dB_rr_Sqrt(1,dbl_mb(agr(1)), 451c > dbl_mb(agr(1))) 452 call D3dB_rr_Sqrt1(1,dbl_mb(agr(1))) 453 454c ***** calculate rhodn **** 455 call D3dB_r_Copy(1,dn(1,2),dbl_mb(rhodn(1))) 456 call D3dB_r_Zero_Ends(1,dbl_mb(rhodn(1))) 457c call D3dB_r_SMul(1,scal1,dbl_mb(rhodn(1)), 458c > dbl_mb(rhodn(1))) 459 call D3dB_r_SMul1(1,scal1,dbl_mb(rhodn(1))) 460 call D3dB_rc_fft3f(1,dbl_mb(rhodn(1))) 461 call mask_C(0,dbl_mb(rhodn(1))) 462 463 464c ***** calculate grdn = grad ndn **** 465 call D3dB_ic_Mul(1,dbl_mb(G_indx(1)), 466 > dbl_mb(rhodn(1)), 467 > dbl_mb(grdnx(1))) 468 call D3dB_ic_Mul(1,dbl_mb(G_indx(2)), 469 > dbl_mb(rhodn(1)), 470 > dbl_mb(grdny(1))) 471 call D3dB_ic_Mul(1,dbl_mb(G_indx(3)), 472 > dbl_mb(rhodn(1)), 473 > dbl_mb(grdnz(1))) 474 call D3dB_cr_fft3b(1,dbl_mb(grdnx(1))) 475 call D3dB_cr_fft3b(1,dbl_mb(grdny(1))) 476 call D3dB_cr_fft3b(1,dbl_mb(grdnz(1))) 477 478c ***** calculate agrdn = |grad ndn| **** 479 call D3dB_rr_Sqr(1,dbl_mb(grdnx(1)), 480 > dbl_mb(agr(2))) 481 call D3dB_rr_Sqr(1,dbl_mb(grdny(1)), 482 > dbl_mb(tmp(1))) 483 484c call D3dB_rr_Sum(1,dbl_mb(tmp(1)), 485c > dbl_mb(agr(2)), 486c > dbl_mb(agr(2))) 487 call D3dB_rr_Sum2(1,dbl_mb(tmp(1)),dbl_mb(agr(2))) 488 489 call D3dB_rr_Sqr(1,dbl_mb(grdnz(1)), 490 > dbl_mb(tmp(1))) 491 492c call D3dB_rr_Sum(1,dbl_mb(tmp(1)), 493c > dbl_mb(agr(2)), 494c > dbl_mb(agr(2))) 495 call D3dB_rr_Sum2(1,dbl_mb(tmp(1)),dbl_mb(agr(2))) 496 497c call D3dB_rr_Sqrt(1,dbl_mb(agr(2)), 498c > dbl_mb(agr(2))) 499 call D3dB_rr_Sqrt1(1,dbl_mb(agr(2))) 500 501 502c ***** calculate agr = |grad nup +grad ndn| **** 503 call D3dB_rr_Sum(1,dbl_mb(grupx(1)), 504 > dbl_mb(grdnx(1)), 505 > dbl_mb(grallx(1))) 506 call D3dB_rr_Sum(1,dbl_mb(grupy(1)), 507 > dbl_mb(grdny(1)), 508 > dbl_mb(grally(1))) 509 call D3dB_rr_Sum(1,dbl_mb(grupz(1)), 510 > dbl_mb(grdnz(1)), 511 > dbl_mb(grallz(1))) 512 513 call D3dB_rr_Sqr(1,dbl_mb(grallx(1)), 514 > dbl_mb(agr(3))) 515 call D3dB_rr_Sqr(1,dbl_mb(grally(1)), 516 > dbl_mb(tmp(1))) 517 518c call D3dB_rr_Sum(1,dbl_mb(tmp(1)), 519c > dbl_mb(agr(3)), 520c > dbl_mb(agr(3))) 521 call D3dB_rr_Sum2(1,dbl_mb(tmp(1)),dbl_mb(agr(3))) 522 523 call D3dB_rr_Sqr(1,dbl_mb(grallz(1)), 524 > dbl_mb(tmp(1))) 525 526c call D3dB_rr_Sum(1,dbl_mb(tmp(1)), 527c > dbl_mb(agr(3)), 528c > dbl_mb(agr(3))) 529 call D3dB_rr_Sum2(1,dbl_mb(tmp(1)),dbl_mb(agr(3))) 530 531c call D3dB_rr_Sqrt(1,dbl_mb(agr(3)), 532c > dbl_mb(agr(3))) 533 call D3dB_rr_Sqrt1(1,dbl_mb(agr(3))) 534 if (use_nwpwxc) then 535c 536c Copy |grad nup|->grad(1), |grad n|->grad(2), and |grad ndn|->grad(3) 537c 538 call D3dB_r_Copy(1,dbl_mb(agr(1)),dbl_mb(grad(1))) 539 call D3dB_r_Copy(1,dbl_mb(agr(2)),dbl_mb(grad(1)+2*n2ft3d)) 540 call D3dB_r_Copy(1,dbl_mb(agr(3)),dbl_mb(grad(1)+n2ft3d)) 541c 542c Replace |grad x| with |grad x|^2 543c 544 call D3dB_rr_Sqr1(1,dbl_mb(grad(1))) 545 call D3dB_rr_Sqr1(1,dbl_mb(grad(1)+2*n2ft3d)) 546 call D3dB_rr_Sqr1(1,dbl_mb(grad(1)+n2ft3d)) 547c 548c Replace |grad n|^2 with |(grad nup|grad ndn)|^2 549c 550 call D3dB_rr_Sub2(1,dbl_mb(grad(1)),dbl_mb(grad(1)+n2ft3d)) 551 call D3dB_rr_Sub2(1,dbl_mb(grad(1)+2*n2ft3d), 552 + dbl_mb(grad(1)+n2ft3d)) 553 call D3dB_r_SMul1(1,0.5d0,dbl_mb(grad(1)+n2ft3d)) 554c 555c Evaluate the functional 556c 557 call nwpwxc_eval_df(2,n2ft3d,dn,dbl_mb(grad(1)), 558 > dumtau,xce, 559 > dbl_mb(fn(1)),dbl_mb(fdn(1)),dumtau) 560c 561c Replace f with the energy density f/n 562c 563 do u = 1, n2ft3d 564 xce(u) = xce(u)/(dn(u,1)+dn(u,2)+dncut) 565 enddo 566c 567c Replace (df/d|grad nup|^2) with (df/d|grad nup|) 568c 569 call D3dB_rr_daxpy(1,(-0.5d0),dbl_mb(fdn(2)),dbl_mb(fdn(1))) 570 call D3dB_rr_Mul2(1,dbl_mb(agr(1)),dbl_mb(fdn(1))) 571 call D3dB_r_SMul1(1,2.0d0,dbl_mb(fdn(1))) 572c 573c Replace (df/d|grad ndn|^2) with (df/d|grad ndn|) 574c 575 call D3dB_rr_daxpy(1,(-0.5d0),dbl_mb(fdn(2)),dbl_mb(fdn(3))) 576 call D3dB_rr_Mul(1,dbl_mb(agr(2)),dbl_mb(fdn(3)), 577 + dbl_mb(gtmp(1)+n2ft3d)) 578 call D3dB_r_SMul1(1,2.0d0,dbl_mb(gtmp(1)+n2ft3d)) 579c 580c Replace (df/d|(grad nup|grad ndn)|^2) with (df/d|grad n|) 581c 582 call D3dB_rr_Mul(1,dbl_mb(agr(3)),dbl_mb(fdn(2)), 583 + dbl_mb(gtmp(1)+2*n2ft3d)) 584c 585c Put the results back into fdn 586c 587 call dcopy(n2ft3d,dbl_mb(gtmp(1)+n2ft3d),1, 588 + dbl_mb(fdn(2)),1) 589 call dcopy(n2ft3d,dbl_mb(gtmp(1)+2*n2ft3d),1, 590 + dbl_mb(fdn(3)),1) 591c 592 else 593 594c ***** calculate **** 595c ***** fdn=(dfx/d|grad nup|,dfx/d|grad ndn|,dfc/d|grad n|) **** 596 if (gga.eq.10) then 597 call gen_PBE96_BW_unrestricted(n2ft3d,dn, 598 > dbl_mb(agr(1)), 599 > x_parameter,c_parameter, 600 > dbl_mb(xce(1)), !*** not used *** 601 > dbl_mb(fn(1)), !*** not used *** 602 > dbl_mb(fdn(1))) 603 else if (gga.eq.11) then 604 call gen_BLYP_BW_unrestricted(n2ft3d,dn, 605 > dbl_mb(agr(1)), 606 > x_parameter,c_parameter, 607 > dbl_mb(xce(1)), !*** not used *** 608 > dbl_mb(fn(1)), !*** not used *** 609 > dbl_mb(fdn(1))) 610 else if (gga.eq.12) then 611 call gen_revPBE_BW_unrestricted(n2ft3d,dn, 612 > dbl_mb(agr(1)), 613 > x_parameter,c_parameter, 614 > dbl_mb(xce(1)), !*** not used *** 615 > dbl_mb(fn(1)), !*** not used *** 616 > dbl_mb(fdn(1))) 617 else if (gga.eq.13) then 618 call gen_PBEsol_BW_unrestricted(n2ft3d,dn, 619 > dbl_mb(agr(1)), 620 > x_parameter,c_parameter, 621 > dbl_mb(xce(1)), !*** not used *** 622 > dbl_mb(fn(1)), !*** not used *** 623 > dbl_mb(fdn(1))) 624 else if (gga.eq.14) then 625 call gen_HSE_BW_unrestricted(n2ft3d,dn, 626 > dbl_mb(agr(1)), 627 > x_parameter,c_parameter, 628 > dbl_mb(xce(1)), !*** not used *** 629 > dbl_mb(fn(1)), !*** not used *** 630 > dbl_mb(fdn(1))) 631 else if (gga.eq.15) then 632 call gen_B3LYP_BW_unrestricted(n2ft3d,dn, 633 > dbl_mb(agr(1)), 634 > x_parameter,c_parameter, 635 > dbl_mb(xce(1)), !*** not used *** 636 > dbl_mb(fn(1)), !*** not used *** 637 > dbl_mb(fdn(1))) 638 else if (gga.eq.16) then 639 call gen_BEEF_BW_unrestricted(n2ft3d,dn, 640 > dbl_mb(agr(1)), 641 > x_parameter,c_parameter,0.6001664769d0, 642 > dbl_mb(xce(1)), !*** not used *** 643 > dbl_mb(fn(1)), !*** not used *** 644 > dbl_mb(fdn(1))) 645 else if (gga.eq.17) then 646 call gen_BEEF_BW_unrestricted(n2ft3d,dn, 647 > dbl_mb(agr(1)), 648 > x_parameter,c_parameter,0.0d0, 649 > dbl_mb(xce(1)), !*** not used *** 650 > dbl_mb(fn(1)), !*** not used *** 651 > dbl_mb(fdn(1))) 652 else 653 call errquit('bad gga',0,0) 654 end if 655 end if ! use_nwpwxc 656 657 658* **** calculate df/d|grad nup|* (grad nup)/|grad nup| **** 659* **** calculate df/d|grad ndn|* (grad ndn)/|grad ndn| **** 660* **** calculate df/d|grad n| * (grad n)/|grad n| **** 661c call D3dB_rr_Divide(1,dbl_mb(grupx(1)), 662c > dbl_mb(agr(1)), 663c > dbl_mb(grupx(1))) 664c call D3dB_rr_Divide(1,dbl_mb(grupy(1)), 665c > dbl_mb(agr(1)), 666c > dbl_mb(grupy(1))) 667c call D3dB_rr_Divide(1,dbl_mb(grupz(1)), 668c > dbl_mb(agr(1)), 669c > dbl_mb(grupz(1))) 670c call D3dB_rr_Divide(1,dbl_mb(grdnx(1)), 671c > dbl_mb(agr(2)), 672c > dbl_mb(grdnx(1))) 673c call D3dB_rr_Divide(1,dbl_mb(grdny(1)), 674c > dbl_mb(agr(2)), 675c > dbl_mb(grdny(1))) 676c call D3dB_rr_Divide(1,dbl_mb(grdnz(1)), 677c > dbl_mb(agr(2)), 678c > dbl_mb(grdnz(1))) 679c call D3dB_rr_Divide(1,dbl_mb(grallx(1)), 680c > dbl_mb(agr(3)), 681c > dbl_mb(grallx(1))) 682c call D3dB_rr_Divide(1,dbl_mb(grally(1)), 683c > dbl_mb(agr(3)), 684c > dbl_mb(grally(1))) 685c call D3dB_rr_Divide(1,dbl_mb(grallz(1)), 686c > dbl_mb(agr(3)), 687c > dbl_mb(grallz(1))) 688 call D3dB_rr_Divide2(1,dbl_mb(agr(1)),dbl_mb(grupx(1))) 689 call D3dB_rr_Divide2(1,dbl_mb(agr(1)),dbl_mb(grupy(1))) 690 call D3dB_rr_Divide2(1,dbl_mb(agr(1)),dbl_mb(grupz(1))) 691 call D3dB_rr_Divide2(1,dbl_mb(agr(2)),dbl_mb(grdnx(1))) 692 call D3dB_rr_Divide2(1,dbl_mb(agr(2)),dbl_mb(grdny(1))) 693 call D3dB_rr_Divide2(1,dbl_mb(agr(2)),dbl_mb(grdnz(1))) 694 call D3dB_rr_Divide2(1,dbl_mb(agr(3)),dbl_mb(grallx(1))) 695 call D3dB_rr_Divide2(1,dbl_mb(agr(3)),dbl_mb(grally(1))) 696 call D3dB_rr_Divide2(1,dbl_mb(agr(3)),dbl_mb(grallz(1))) 697 698 699c call D3dB_rr_Mul(1,dbl_mb(grupx(1)), 700c > dbl_mb(fdn(1)), 701c > dbl_mb(grupx(1))) 702c call D3dB_rr_Mul(1,dbl_mb(grupy(1)), 703c > dbl_mb(fdn(1)), 704c > dbl_mb(grupy(1))) 705c call D3dB_rr_Mul(1,dbl_mb(grupz(1)), 706c > dbl_mb(fdn(1)), 707c > dbl_mb(grupz(1))) 708c call D3dB_rr_Mul(1,dbl_mb(grdnx(1)), 709c > dbl_mb(fdn(2)), 710c > dbl_mb(grdnx(1))) 711c call D3dB_rr_Mul(1,dbl_mb(grdny(1)), 712c > dbl_mb(fdn(2)), 713c > dbl_mb(grdny(1))) 714c call D3dB_rr_Mul(1,dbl_mb(grdnz(1)), 715c > dbl_mb(fdn(2)), 716c > dbl_mb(grdnz(1))) 717c call D3dB_rr_Mul(1,dbl_mb(grallx(1)), 718c > dbl_mb(fdn(3)), 719c > dbl_mb(grallx(1))) 720c call D3dB_rr_Mul(1,dbl_mb(grally(1)), 721c > dbl_mb(fdn(3)), 722c > dbl_mb(grally(1))) 723c call D3dB_rr_Mul(1,dbl_mb(grallz(1)), 724c > dbl_mb(fdn(3)), 725c > dbl_mb(grallz(1))) 726 call D3dB_rr_Mul2(1,dbl_mb(fdn(1)),dbl_mb(grupx(1))) 727 call D3dB_rr_Mul2(1,dbl_mb(fdn(1)),dbl_mb(grupy(1))) 728 call D3dB_rr_Mul2(1,dbl_mb(fdn(1)),dbl_mb(grupz(1))) 729 call D3dB_rr_Mul2(1,dbl_mb(fdn(2)),dbl_mb(grdnx(1))) 730 call D3dB_rr_Mul2(1,dbl_mb(fdn(2)),dbl_mb(grdny(1))) 731 call D3dB_rr_Mul2(1,dbl_mb(fdn(2)),dbl_mb(grdnz(1))) 732 call D3dB_rr_Mul2(1,dbl_mb(fdn(3)),dbl_mb(grallx(1))) 733 call D3dB_rr_Mul2(1,dbl_mb(fdn(3)),dbl_mb(grally(1))) 734 call D3dB_rr_Mul2(1,dbl_mb(fdn(3)),dbl_mb(grallz(1))) 735 736* **** calculate (df/d|grad nup|* (grad nup)/|grad nup|) **** 737* **** + (df/d|grad n| * (grad n)/|grad n|) **** 738 739* **** calculate (df/d|grad ndn|* (grad ndn)/|grad ndn|) **** 740* **** + (df/d|grad n| * (grad n)/|grad n|) **** 741c call D3dB_rr_Sum(1,dbl_mb(grupx(1)), 742c > dbl_mb(grallx(1)), 743c > dbl_mb(grupx(1))) 744c call D3dB_rr_Sum(1,dbl_mb(grupy(1)), 745c > dbl_mb(grally(1)), 746c > dbl_mb(grupy(1))) 747c call D3dB_rr_Sum(1,dbl_mb(grupz(1)), 748c > dbl_mb(grallz(1)), 749c > dbl_mb(grupz(1))) 750c call D3dB_rr_Sum(1,dbl_mb(grdnx(1)), 751c > dbl_mb(grallx(1)), 752c > dbl_mb(grdnx(1))) 753c call D3dB_rr_Sum(1,dbl_mb(grdny(1)), 754c > dbl_mb(grally(1)), 755c > dbl_mb(grdny(1))) 756c call D3dB_rr_Sum(1,dbl_mb(grdnz(1)), 757c > dbl_mb(grallz(1)), 758c > dbl_mb(grdnz(1))) 759 call D3dB_rr_Sum2(1,dbl_mb(grallx(1)),dbl_mb(grupx(1))) 760 call D3dB_rr_Sum2(1,dbl_mb(grally(1)),dbl_mb(grupy(1))) 761 call D3dB_rr_Sum2(1,dbl_mb(grallz(1)),dbl_mb(grupz(1))) 762 call D3dB_rr_Sum2(1,dbl_mb(grallx(1)),dbl_mb(grdnx(1))) 763 call D3dB_rr_Sum2(1,dbl_mb(grally(1)),dbl_mb(grdny(1))) 764 call D3dB_rr_Sum2(1,dbl_mb(grallz(1)),dbl_mb(grdnz(1))) 765 766c call D3dB_r_SMul(1,scal1,dbl_mb(grupx(1)), 767c > dbl_mb(grupx(1))) 768c call D3dB_r_SMul(1,scal1,dbl_mb(grupy(1)), 769c > dbl_mb(grupy(1))) 770c call D3dB_r_SMul(1,scal1,dbl_mb(grupz(1)), 771c > dbl_mb(grupz(1))) 772c call D3dB_r_SMul(1,scal1,dbl_mb(grdnx(1)), 773c > dbl_mb(grdnx(1))) 774c call D3dB_r_SMul(1,scal1,dbl_mb(grdny(1)), 775c > dbl_mb(grdny(1))) 776c call D3dB_r_SMul(1,scal1,dbl_mb(grdnz(1)), 777c > dbl_mb(grdnz(1))) 778 call D3dB_r_SMul1(1,scal1,dbl_mb(grupx(1))) 779 call D3dB_r_SMul1(1,scal1,dbl_mb(grupy(1))) 780 call D3dB_r_SMul1(1,scal1,dbl_mb(grupz(1))) 781 call D3dB_r_SMul1(1,scal1,dbl_mb(grdnx(1))) 782 call D3dB_r_SMul1(1,scal1,dbl_mb(grdny(1))) 783 call D3dB_r_SMul1(1,scal1,dbl_mb(grdnz(1))) 784 785 call D3dB_r_Zero_Ends(1,dbl_mb(grupx(1))) 786 call D3dB_r_Zero_Ends(1,dbl_mb(grupy(1))) 787 call D3dB_r_Zero_Ends(1,dbl_mb(grupz(1))) 788 call D3dB_r_Zero_Ends(1,dbl_mb(grdnx(1))) 789 call D3dB_r_Zero_Ends(1,dbl_mb(grdny(1))) 790 call D3dB_r_Zero_Ends(1,dbl_mb(grdnz(1))) 791 792* **** put sums in k-space *** 793 call D3dB_rc_fft3f(1,dbl_mb(grupx(1))) 794 call D3dB_rc_fft3f(1,dbl_mb(grupy(1))) 795 call D3dB_rc_fft3f(1,dbl_mb(grupz(1))) 796 call D3dB_rc_fft3f(1,dbl_mb(grdnx(1))) 797 call D3dB_rc_fft3f(1,dbl_mb(grdny(1))) 798 call D3dB_rc_fft3f(1,dbl_mb(grdnz(1))) 799 800 801* **** W(u,s) = Sum(G) [i*G(u)*dcongj(rhoup)*grup(s)] **** 802* **** where grup(1)=grupx,grup(2)=grupy,grup(3)=grupz **** 803 !call Pack_c_pack(0,dbl_mb(rhoup(1))) 804 call mask_C(0,dbl_mb(rhoup(1))) 805 do u=1,3 806 807* **** agr = i*G(u)*grupx **** 808 call D3dB_ic_Mul(1,dbl_mb(G_indx(u)), 809 > dbl_mb(grupx(1)), 810 > dbl_mb(agr(1))) 811 call mask_C(0,dbl_mb(agr(1))) 812 call D3dB_cc_dot(1,dbl_mb(rhoup(1)),dbl_mb(agr(1)),sum) 813 !call Pack_c_pack(0,dbl_mb(agr(1))) 814 !call Pack_cc_dot(0,dbl_mb(rhoup(1)),dbl_mb(agr(1)),sum) 815 W(u,1) = sum*omega 816 817* **** agr = i*G(u)*grupy **** 818 call D3dB_ic_Mul(1,dbl_mb(G_indx(u)), 819 > dbl_mb(grupy(1)), 820 > dbl_mb(agr(1))) 821 call mask_C(0,dbl_mb(agr(1))) 822 call D3dB_cc_dot(1,dbl_mb(rhoup(1)),dbl_mb(agr(1)),sum) 823 !call Pack_c_pack(0,dbl_mb(agr(1))) 824 !call Pack_cc_dot(0,dbl_mb(rhoup(1)),dbl_mb(agr(1)),sum) 825 W(u,2) = sum*omega 826 827* **** agr = i*G(u)*grupz **** 828 call D3dB_ic_Mul(1,dbl_mb(G_indx(u)), 829 > dbl_mb(grupz(1)), 830 > dbl_mb(agr(1))) 831 call mask_C(0,dbl_mb(agr(1))) 832 call D3dB_cc_dot(1,dbl_mb(rhoup(1)),dbl_mb(agr(1)),sum) 833 !call Pack_c_pack(0,dbl_mb(agr(1))) 834 !call Pack_cc_dot(0,dbl_mb(rhoup(1)),dbl_mb(agr(1)),sum) 835 W(u,3) = sum*omega 836 end do 837 838 839* **** W(u,s) = Sum(G) [i*G(u)*dcongj(rhodn)*grup(s)] **** 840* **** where grdn(1)=grdnx,grup(2)=grdny,grup(3)=grdnz **** 841 !call Pack_c_pack(0,dbl_mb(rhodn(1))) 842 call mask_C(0,dbl_mb(rhodn(1))) 843 do u=1,3 844 845* **** agr = i*G(u)*grdnx **** 846 call D3dB_ic_Mul(1,dbl_mb(G_indx(u)), 847 > dbl_mb(grdnx(1)), 848 > dbl_mb(agr(1))) 849 call mask_C(0,dbl_mb(agr(1))) 850 call D3dB_cc_dot(1,dbl_mb(rhodn(1)),dbl_mb(agr(1)),sum) 851 !call Pack_c_pack(0,dbl_mb(agr(1))) 852 !call Pack_cc_dot(0,dbl_mb(rhodn(1)),dbl_mb(agr(1)),sum) 853 W(u,1) = W(u,1) + sum*omega 854 855* **** agr = i*G(u)*grdny **** 856 call D3dB_ic_Mul(1,dbl_mb(G_indx(u)), 857 > dbl_mb(grdny(1)), 858 > dbl_mb(agr(1))) 859 call mask_C(0,dbl_mb(agr(1))) 860 call D3dB_cc_dot(1,dbl_mb(rhodn(1)),dbl_mb(agr(1)),sum) 861 !call Pack_c_pack(0,dbl_mb(agr(1))) 862 !call Pack_cc_dot(0,dbl_mb(rhodn(1)),dbl_mb(agr(1)),sum) 863 W(u,2) = W(u,2) + sum*omega 864 865* **** agr = i*G(u)*grdnz **** 866 call D3dB_ic_Mul(1,dbl_mb(G_indx(u)), 867 > dbl_mb(grdnz(1)), 868 > dbl_mb(agr(1))) 869 call mask_C(0,dbl_mb(agr(1))) 870 call D3dB_cc_dot(1,dbl_mb(rhodn(1)),dbl_mb(agr(1)),sum) 871 !call Pack_c_pack(0,dbl_mb(agr(1))) 872 !call Pack_cc_dot(0,dbl_mb(rhodn(1)),dbl_mb(agr(1)),sum) 873 W(u,3) = W(u,3) + sum*omega 874 end do 875 876 877 878* **** deallocate temporary memory **** 879 value = BA_pop_stack(xce(2)) 880 value = value.and.BA_pop_stack(xfdn(2)) 881 value = value.and.BA_pop_stack(xfn(2)) 882 value = value.and.BA_pop_stack(gtmp(2)) 883 value = value.and.BA_pop_stack(grad(2)) 884 value = value.and.BA_pop_stack(xagr(2)) 885 886 value = value.and.BA_pop_stack(grallz(2)) 887 value = value.and.BA_pop_stack(grally(2)) 888 value = value.and.BA_pop_stack(grallx(2)) 889 value = value.and.BA_pop_stack(grdnz(2)) 890 value = value.and.BA_pop_stack(grdny(2)) 891 value = value.and.BA_pop_stack(grdnx(2)) 892 value = value.and.BA_pop_stack(rhodn(2)) 893 value = value.and.BA_pop_stack(grupz(2)) 894 value = value.and.BA_pop_stack(grupy(2)) 895 value = value.and.BA_pop_stack(grupx(2)) 896 value = value.and.BA_pop_stack(rhoup(2)) 897 if (.not. value) call errquit('cannot pop stack memory',0,0) 898 899 900 end if 901 902 903* **** stress(u,v) = Sum(s){W(u,s)*hm(s,v) } **** 904 do v=1,3 905 do u=1,3 906 stress(u,v) = 0.0d0 907 do s=1,3 908 stress(u,v) = stress(u,v) + W(u,s)*hm(s,v) 909 end do 910 end do 911 end do 912 913 914 call nwpw_timing_end(4) 915 916 return 917 end 918 919 920 921 922