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