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