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