1* 2* $Id$ 3* 4 5 6* *********************************** 7* * * 8* * electron_init * 9* * * 10* *********************************** 11 subroutine electron_init() 12 implicit none 13 14#include "bafdecls.fh" 15#include "electron_common.fh" 16#include "errquit.fh" 17 18 19* **** electron_counter common block **** 20 integer counter 21 common / electron_counter / counter 22 23 24* **** local variables **** 25 logical value 26 integer n2ft3d 27 28* **** external functions **** 29 logical ion_chargeexist,ion_mmexist,psp_pawexist 30 integer psi_ispin,psi_ne,psi_neq,control_version 31 external ion_chargeexist,ion_mmexist,psp_pawexist 32 external psi_ispin,psi_ne,psi_neq,control_version 33 34!$OMP MASTER 35 counter = 0 36!$OMP END MASTER 37 38 ispin = psi_ispin() 39 ne(1) = psi_ne(1) 40 ne(2) = psi_ne(2) 41 neq(1) = psi_neq(1) 42 neq(2) = psi_neq(2) 43 field_exist = ion_chargeexist().or.ion_mmexist() 44 confine_exist = .false. 45 46* **** get nfft3d, and n2ft3d **** 47 call Pack_npack(1,npack1) 48 call Pack_npack(0,npack0) 49 call D3dB_nfft3d(1,nfft3d) 50 n2ft3d = 2*nfft3d 51 52c paw_exist = psp_pawexist() 53 54* **** allocate memory **** 55 value = BA_alloc_get(mt_dcpl,npack1*(neq(1)+neq(2)), 56 > 'Hpsi_k',Hpsi_k(2),Hpsi_k(1)) 57 value = value.and. 58 > BA_alloc_get(mt_dbl,n2ft3d*(neq(1)+neq(2)), 59 > 'psi_r',psi_r(2),psi_r(1)) 60 value = value.and. 61 > BA_alloc_get(mt_dcpl,npack0, 62 > 'vl2',vl(2),vl(1)) 63 64 value = value.and. 65 > BA_alloc_get(mt_dbl,n2ft3d, 66 > 'vl_lr',vl_lr(2),vl_lr(1)) 67 value = value.and. 68 > BA_alloc_get(mt_dbl,n2ft3d, 69 > 'v_field',v_field(2),v_field(1)) 70 value = value.and. 71 > BA_alloc_get(mt_dbl,2*n2ft3d, 72 > 'vall',vall(2),vall(1)) 73 74 75 if (control_version().eq.3) then 76 value = value.and. 77 > BA_alloc_get(mt_dcpl,npack0, 78 > 'vc',vc(2),vc(1)) 79 end if 80 81 if (control_version().eq.4) then 82 value = value.and. 83 > BA_alloc_get(mt_dcpl,nfft3d, 84 > 'vc',vc(2),vc(1)) 85 end if 86 87 value = value.and. 88 > BA_alloc_get(mt_dbl,2*n2ft3d, 89 > 'xcp',xcp(2),xcp(1)) 90 value = value.and. 91 > BA_alloc_get(mt_dbl,2*n2ft3d, 92 > 'xce',xce(2),xce(1)) 93 if (.not. value) 94 > call errquit('electron_init: out of heap memory',0, MA_ERR) 95 call dcopy(n2ft3d*(neq(1)+neq(2)),0.0d0,0,dbl_mb(psi_r(1)),1) 96 call dcopy(2*n2ft3d,0.0d0,0,dbl_mb(vall(1)),1) 97 call dcopy(n2ft3d,0.0d0,0,dbl_mb(v_field(1)),1) 98 call dcopy(n2ft3d,0.0d0,0,dbl_mb(vl_lr(1)),1) 99 call dcopy(2*n2ft3d,0.0d0,0,dbl_mb(xcp(1)),1) 100 call dcopy(2*n2ft3d,0.0d0,0,dbl_mb(xce(1)),1) 101 102 return 103 end 104 105 106* *********************************** 107* * * 108* * electron_finalize * 109* * * 110* *********************************** 111 subroutine electron_finalize() 112 implicit none 113#include "errquit.fh" 114 115#include "bafdecls.fh" 116#include "electron_common.fh" 117 118 119* **** local variables **** 120 logical value 121 122* **** free heap memory **** 123 value = BA_free_heap(Hpsi_k(2)) 124 value = value.and. 125 > BA_free_heap(psi_r(2)) 126 value = value.and. 127 > BA_free_heap(vl(2)) 128 value = value.and. 129 > BA_free_heap(vl_lr(2)) 130 value = value.and. 131 > BA_free_heap(v_field(2)) 132 value = value.and. 133 > BA_free_heap(vall(2)) 134 value = value.and. 135 > BA_free_heap(vc(2)) 136 value = value.and. 137 > BA_free_heap(xcp(2)) 138 value = value.and. 139 > BA_free_heap(xce(2)) 140 if (.not. value) 141 > call errquit('electron_init: error freeing heap memory',0, 142 & MA_ERR) 143 144 return 145 end 146 147* *********************************** 148* * * 149* * electron_count * 150* * * 151* *********************************** 152 integer function electron_count() 153 implicit none 154 155* **** electron_counter common block **** 156 integer counter 157 common / electron_counter / counter 158 159 electron_count = counter 160 return 161 end 162 163 164* *********************************** 165* * * 166* * electron_run * 167* * * 168* *********************************** 169 subroutine electron_run(psi_k,dn,dng,dnall,fractional,occ) 170 implicit none 171 complex*16 psi_k(*) 172 real*8 dn(*) 173 complex*16 dng(*) 174 real*8 dnall(*) 175 logical fractional 176 real*8 occ(*) 177 178* **** electron_counter common block **** 179 integer counter 180 common / electron_counter / counter 181 182!$OMP MASTER 183 counter = counter+1 184!$OMP END MASTER 185 186 call electron_gen_psi_r(psi_k) 187 call electron_gen_densities(psi_k,dn,dng,dnall,fractional,occ) 188 call electron_gen_scf_potentials(dn,dng,dnall) 189 call electron_gen_Hpsi_k(psi_k) 190 191 return 192 end 193 194* *********************************** 195* * * 196* * electron_genrho * 197* * * 198* *********************************** 199 subroutine electron_genrho(psi_k,dn,fractional,occ) 200 implicit none 201 complex*16 psi_k(*) 202 real*8 dn(*) 203 logical fractional 204 real*8 occ(*) 205 206 call electron_gen_psi_r(psi_k) 207 call electron_gen_density(psi_k,dn,fractional,occ) 208 return 209 end 210 211 212* *********************************** 213* * * 214* * electron_run_orb * 215* * * 216* *********************************** 217 subroutine electron_run_orb(i,psi_k) 218 implicit none 219 integer i 220 complex*16 psi_k(*) 221 222 call electron_gen_psi_r_orb(i,psi_k) 223 call electron_gen_Hpsi_k_orb(i,psi_k) 224 225 return 226 end 227 228* *********************************** 229* * * 230* * electron_sd_update * 231* * * 232* *********************************** 233 subroutine electron_sd_update(psi1,psi2,dte) 234 implicit none 235 complex*16 psi1(*),psi2(*) 236 real*8 dte 237 238#include "bafdecls.fh" 239#include "electron_common.fh" 240 241 242 243c call ke_Precondition(npack1,(ne(1)+ne(2)),psi1,dcpl_mb(Hpsi_k(1))) 244c call electron_sd_subupdate(npack1,(ne(1)+ne(2)), 245c > psi1,psi2,dcpl_mb(Hpsi_k(1)), 246c > dte) 247 call dcopy(2*(npack1)*(neq(1)+neq(2)),psi1,1,psi2,1) 248 call daxpy(2*(npack1)*(neq(1)+neq(2)), 249 > (-dte), 250 > dcpl_mb(Hpsi_k(1)),1, 251 > psi2,1) 252 253 return 254 end 255* *********************************** 256* * * 257* * electron_cpmd_update * 258* * * 259* *********************************** 260 subroutine electron_cpmd_update(psi0,psi1,psi2,hml,dte) 261 implicit none 262 complex*16 psi0(*),psi1(*),psi2(*) 263 real*8 hml(*) 264 real*8 dte 265 266#include "bafdecls.fh" 267#include "electron_common.fh" 268 269 !**** psi2 = 2*psi1 - psi0 + dte*Hpsi **** 270 !**** - note that Hpsi = minus the gradient in electron **** 271 272 !**** rotate Hpsi **** 273c call electron_gen_Hpsi_k(psi1) 274 call Dneall_fmf_Multiply(0,dcpl_mb(Hpsi_k(1)),npack1, 275 > hml,(-dte), 276 > psi2,0.0d0) 277 call daxpy(2*(npack1)*(neq(1)+neq(2)), 278 > 1.0d0, 279 > psi1,1, 280 > psi2,1) 281 282c call daxpy(2*(npack1)*(neq(1)+neq(2)), 283c > -1.0d0, 284c > psi0,1, 285c > psi2,1) 286c call daxpy(2*(npack1)*(neq(1)+neq(2)), 287c > (-dte), 288c > dcpl_mb(Hpsi_k(1)),1, 289c > psi2,1) 290 291 return 292 end 293 294 295 296c subroutine electron_sd_subupdate(npack1,nn, 297c > psi1,psi2,Hpsi,dte) 298c implicit none 299c integer npack1,nn 300c complex*16 psi1(npack1,nn) 301c complex*16 psi2(npack1,nn) 302c complex*16 Hpsi(npack1,nn) 303c real*8 dte 304c 305c integer n 306c 307c* ************************************ 308c* **** do a steepest descent step **** 309c* ************************************ 310c do n=1,nn 311c call Pack_c_SMul(1,(-dte),Hpsi(1,n),psi2(1,n)) 312c call Pack_cc_Sum(1,psi2(1,n),psi1(1,n),psi2(1,n)) 313c end do 314c 315c return 316c end 317 318* *********************************** 319* * * 320* * electron_energy * 321* * * 322* *********************************** 323 real*8 function electron_energy(psi_k,dn,dng,dnall,fractional,occ) 324 implicit none 325 complex*16 psi_k(*) 326 real*8 dn(*) 327 complex*16 dng(*) 328 real*8 dnall(*) 329 logical fractional 330 real*8 occ(*) 331 332#include "bafdecls.fh" 333#include "electron_common.fh" 334 335 336* **** local variables **** 337 integer n2ft3d 338 integer ii,ms,n1(2),n2(2),nx,ny,nz 339 real*8 sum,eorbit,ehartr,exc,pxc,exc2,pxc2,dv,eion_core 340 real*8 ehartree_atom,ecmp_cmp,ecmp_pw,exc_atom,pxc_atom,eke_core 341 real*8 total_energy 342 343 real*8 e1,e2 344 common /eenergy_tmp_common/ e1,e2 345 346 347* **** external functions ***** 348 logical pspw_SIC,pspw_SIC_relaxed,psp_U_psputerm,psp_pawexist 349 logical pspw_HFX,pspw_HFX_relaxed,meta_found,nwpw_meta_gga_on 350 integer control_version 351 real*8 lattice_omega,coulomb_e,electron_ehartree2 352 real*8 nwpw_meta_gga_pxc,psp_kinetic_core,psp_ion_core 353 real*8 psp_hartree_atom,psp_hartree_cmp_cmp 354 real*8 psp_hartree_cmp_pw 355 external pspw_SIC,pspw_SIC_relaxed,psp_U_psputerm,psp_pawexist 356 external pspw_HFX,pspw_HFX_relaxed,meta_found,nwpw_meta_gga_on 357 external control_version 358 external lattice_omega,coulomb_e,electron_ehartree2 359 external nwpw_meta_gga_pxc,psp_kinetic_core,psp_ion_core 360 external psp_hartree_atom,psp_hartree_cmp_cmp 361 external psp_hartree_cmp_pw 362 363 call D3dB_nx(1,nx) 364 call D3dB_ny(1,ny) 365 call D3dB_nz(1,nz) 366 367 368 dv = lattice_omega()/dble(nx*ny*nz) 369 370 n2ft3d = 2*nfft3d 371 n1(1) = 1 372 n1(2) = neq(1) + 1 373 n2(1) = neq(1) 374 n2(2) = neq(1) + neq(2) 375 376 377* *** get orbital energies **** 378!$OMP MASTER 379 e1 = 0.0d0 380 e2 = 0.0d0 381!$OMP END MASTER 382!$OMP BARRIER 383 eorbit = 0.0d0 384 if (fractional) then 385 do ms=1,ispin 386 do ii=n1(ms),n2(ms) 387 call Pack_cc_idot(1, 388 > psi_k(1+(ii-1)*npack1), 389 > dcpl_mb(Hpsi_k(1)+(ii-1)*npack1), 390 > e1) 391!$OMP MASTER 392 !eorbit = eorbit + sum*occ(ii) 393 e2 = e2 + e1*occ(ii) 394!$OMP END MASTER 395 end do 396 end do 397 else 398 do ms=1,ispin 399 do ii=n1(ms),n2(ms) 400 call Pack_cc_idot(1, 401 > psi_k(1+(ii-1)*npack1), 402 > dcpl_mb(Hpsi_k(1)+(ii-1)*npack1), 403 > e1) 404!$OMP MASTER 405 !eorbit = eorbit + sum 406 e2 = e2 + e1 407!$OMP END MASTER 408 end do 409 end do 410 end if 411!$OMP BARRIER 412 call Parallel_SumAll(e2) 413 eorbit = e2 414 if (ispin.eq.1) eorbit = eorbit+eorbit 415 416 417 418* **** get coulomb energy **** 419 ehartr = 0.0d0 420 if (control_version().eq.3) ehartr = coulomb_e(dng) 421 if (control_version().eq.4) ehartr = electron_ehartree2(dn) 422 423 424 425* **** get exchange-correlation energy **** 426 call D3dB_rr_dot(1,dnall(1), 427 > dbl_mb(xce(1)), 428 > e1) 429 call D3dB_rr_dot(1,dn(1), 430 > dbl_mb(xcp(1)), 431 . e2) 432 exc = e1 433 pxc = e2 434 if (ispin.eq.1) then 435 exc= exc + exc 436 pxc= pxc + pxc 437 else 438 call D3dB_rr_dot(1,dnall(1+n2ft3d), 439 > dbl_mb(xce(1)), 440 > e1) 441 call D3dB_rr_dot(1,dn(1+n2ft3d), 442 > dbl_mb(xcp(1)+n2ft3d), 443 > e2) 444 exc2 = e1 445 pxc2 = e2 446 exc = exc + exc2 447 pxc = pxc + pxc2 448 end if 449 exc = exc*dv 450 pxc = pxc*dv 451 452* **** meta_GGA energy **** 453 if (nwpw_meta_gga_on()) then 454 pxc = pxc + nwpw_meta_gga_pxc(ispin,neq,psi_k) 455 end if 456 457 total_energy = eorbit + exc - ehartr - pxc 458 459 460 461* **** PAW ee terms **** 462 if (psp_pawexist()) then 463 eke_core = psp_kinetic_core() 464 eion_core = psp_ion_core() 465 ehartree_atom = psp_hartree_atom(ispin,neq,psi_k) 466 ecmp_cmp = psp_hartree_cmp_cmp(ispin) 467 ecmp_pw = psp_hartree_cmp_pw(ispin,dng,dn) 468 call psp_xc_atom(ispin,neq,psi_k,exc_atom,pxc_atom) 469 470 total_energy = total_energy + exc_atom - pxc_atom 471 > - ehartree_atom - ecmp_cmp - ecmp_pw 472c > + eke_core + eion_core 473 end if 474 475* **** SIC corrections **** 476 if (pspw_SIC()) then 477 if (pspw_SIC_relaxed()) then 478 call pspw_energy_SIC(ispin,dbl_mb(psi_r(1)), 479 > ehsic, 480 > phsic, 481 > exsic, 482 > pxsic) 483 total_energy = total_energy + ehsic + exsic - phsic - pxsic 484 end if 485 end if 486 487* **** HFX energy **** 488 if (pspw_HFX()) then 489 if (pspw_HFX_relaxed()) then 490 call pspw_energy_HFX(ispin,dbl_mb(psi_r(1)), 491 > ehfx, 492 > phfx) 493 total_energy = total_energy + ehfx - phfx 494 end if 495 end if 496 497* **** DFT+U energy **** 498 if (psp_U_psputerm()) then 499 call psp_U_psputerm_energy(edftu,pdftu) 500 total_energy = total_energy + edftu - pdftu 501 end if 502 503 504* **** metadynamics energy **** 505 if (meta_found()) then 506 call meta_energypotential(ispin,neq,psi_k,emeta,pmeta) 507 total_energy = total_energy + emeta - pmeta 508 end if 509 510 511* **** total energy **** 512 electron_energy = total_energy 513 514 return 515 end 516 517 518* *********************************** 519* * * 520* * electron_eorbit_noocc * 521* * * 522* *********************************** 523 real*8 function electron_eorbit_noocc(psi_k) 524 implicit none 525 complex*16 psi_k(*) 526 527#include "bafdecls.fh" 528#include "electron_common.fh" 529 530* **** local variables **** 531 integer ii,ms,n1(2),n2(2) 532 real*8 sum,eorbit 533 534 real*8 e1,e2 535 common /eenergy_tmp_common/ e1,e2 536 537 538 n1(1) = 1 539 n1(2) = neq(1) + 1 540 n2(1) = neq(1) 541 n2(2) = neq(1) + neq(2) 542 543* *** get orbital energies **** 544!$OMP MASTER 545 e1 = 0.0d0 546 e2 = 0.0d0 547!$OMP END MASTER 548!$OMP BARRIER 549 !eorbit = 0.0d0 550 do ms=1,ispin 551 do ii=n1(ms),n2(ms) 552 call Pack_cc_idot(1,psi_k(1+(ii-1)*npack1), 553 > dcpl_mb(Hpsi_k(1)+(ii-1)*npack1), 554 > e1) 555!$OMP MASTER 556 e2 = e2 + e1 557!$OMP END MASTER 558 end do 559 end do 560!$OMP BARRIER 561 call Parallel_SumAll(e2) 562 eorbit = e2 563 if (ispin.eq.1) eorbit = eorbit+eorbit 564 565 electron_eorbit_noocc = eorbit 566 return 567 end 568 569 570 571* *********************************** 572* * * 573* * electron_eorbit * 574* * * 575* *********************************** 576 real*8 function electron_eorbit(psi_k,fractional,occ) 577 implicit none 578 complex*16 psi_k(*) 579 logical fractional 580 real*8 occ(*) 581 582#include "bafdecls.fh" 583#include "electron_common.fh" 584 585 586* **** local variables **** 587 integer ii,ms,n1(2),n2(2) 588 real*8 sum,eorbit 589 590 common /eelectron_ejtmp/ sum,eorbit 591 592 n1(1) = 1 593 n1(2) = neq(1) + 1 594 n2(1) = neq(1) 595 n2(2) = neq(1) + neq(2) 596 597 598* *** get orbital energies **** 599 eorbit = 0.0d0 600 if (fractional) then 601 do ms=1,ispin 602 do ii=n1(ms),n2(ms) 603 call Pack_cc_idot(1,psi_k(1+(ii-1)*npack1), 604 > dcpl_mb(Hpsi_k(1)+(ii-1)*npack1), 605 > sum) 606!$OMP MASTER 607 eorbit = eorbit + sum*occ(ii) 608!$OMP END MASTER 609 end do 610 end do 611 else 612 do ms=1,ispin 613 do ii=n1(ms),n2(ms) 614 call Pack_cc_idot(1,psi_k(1+(ii-1)*npack1), 615 > dcpl_mb(Hpsi_k(1)+(ii-1)*npack1), 616 > sum) 617!$OMP MASTER 618 eorbit = eorbit + sum 619!$OMP END MASTER 620 end do 621 end do 622 end if 623!$OMP MASTER 624 if (ispin.eq.1) eorbit = eorbit+eorbit 625!$OMP END MASTER 626!$BARRIER 627 call Parallel_SumAll(eorbit) 628 629 electron_eorbit = eorbit 630 631 return 632 end 633 634 635* *********************************** 636* * * 637* * electron_ehartree * 638* * * 639* *********************************** 640 real*8 function electron_ehartree(dng) 641 implicit none 642 complex*16 dng(*) 643 644 645* **** external functions **** 646 real*8 coulomb_e 647 external coulomb_e 648 649 electron_ehartree = coulomb_e(dng) 650 651 return 652 end 653 654* *********************************** 655* * * 656* * electron_ehartree2 * 657* * * 658* *********************************** 659 real*8 function electron_ehartree2(dn) 660 implicit none 661 real*8 dn(*) 662 663#include "bafdecls.fh" 664#include "electron_common.fh" 665 666* **** local variables **** 667 real*8 ehartr, ehart1,ehart2,dv 668 integer nx,ny,nz 669 670 real*8 e1,e2 671 common /eenergy_tmp_common/ e1,e2 672 673* ***** external functions **** 674 real*8 lattice_omega 675 external lattice_omega 676 677 call D3dB_nx(1,nx) 678 call D3dB_ny(1,ny) 679 call D3dB_nz(1,nz) 680 dv = lattice_omega()/dble(nx*ny*nz) 681 682 call D3dB_rr_dot(1,dn(1), 683 > dcpl_mb(vc(1)), 684 > e1) 685 call D3dB_rr_dot(1,dn(1+(ispin-1)*2*nfft3d), 686 > dcpl_mb(vc(1)), 687 > e2) 688 ehartr = 0.5d0*(e1+e2)*dv 689 690 electron_ehartree2 = ehartr 691 692 return 693 end 694 695* *********************************** 696* * * 697* * electron_exc * 698* * * 699* *********************************** 700 real*8 function electron_exc(dnall) 701 implicit none 702 real*8 dnall(*) 703 704#include "bafdecls.fh" 705#include "electron_common.fh" 706 707 708* **** local variables **** 709 integer nx,ny,nz 710 real*8 exc,exc2,dv 711 712 common /eelectron_ejtmp/ exc,exc2 713 714* **** external functions **** 715 real*8 lattice_omega 716 external lattice_omega 717 718 719 call D3dB_nx(1,nx) 720 call D3dB_ny(1,ny) 721 call D3dB_nz(1,nz) 722 723 dv = lattice_omega()/dble(nx*ny*nz) 724 725* **** get exchange-correlation energy **** 726 call D3dB_rr_dot(1,dnall, 727 > dbl_mb(xce(1)), 728 > exc) 729 if (ispin.eq.1) then 730!$OMP MASTER 731 exc= exc + exc 732!$OMP END MASTER 733 else 734 call D3dB_rr_dot(1,dnall(1+2*nfft3d), 735 > dbl_mb(xce(1)), 736 > exc2) 737!$OMP MASTER 738 exc= exc + exc2 739!$OMP END MASTER 740 end if 741!$OMP MASTER 742 exc = exc*dv 743!$OMP END MASTER 744!$OMP BARRIER 745 746 electron_exc = exc 747 748 return 749 end 750 751 752* *********************************** 753* * * 754* * electron_pxc * 755* * * 756* *********************************** 757 real*8 function electron_pxc(dn) 758 implicit none 759 real*8 dn(*) 760 761#include "bafdecls.fh" 762#include "electron_common.fh" 763 764 765* **** local variables **** 766 integer nx,ny,nz 767 real*8 pxc,pxc2,dv 768 769 common /eelectron_ejtmp/ pxc,pxc2 770 771* **** external functions ***** 772 real*8 lattice_omega 773 external lattice_omega 774 775 call D3dB_nx(1,nx) 776 call D3dB_ny(1,ny) 777 call D3dB_nz(1,nz) 778 779 dv = lattice_omega()/dble(nx*ny*nz) 780 781* **** get exchange-correlation energy **** 782 call D3dB_rr_dot(1,dn(1), 783 > dbl_mb(xcp(1)), 784 > pxc) 785 if (ispin.eq.1) then 786!$OMP MASTER 787 pxc= pxc + pxc 788!$OMP END MASTER 789 else 790 call D3dB_rr_dot(1,dn(1+2*nfft3d), 791 > dbl_mb(xcp(1)+2*nfft3d), 792 > pxc2) 793!$OMP MASTER 794 pxc= pxc + pxc2 795!$OMP END MASTER 796 end if 797!$OMP MASTER 798 pxc = pxc*dv 799!$OMP END MASTER 800!$OMP BARRIER 801 802 electron_pxc = pxc 803 804 return 805 end 806 807* *********************************** 808* * * 809* * electron_pxc_rho * 810* * * 811* *********************************** 812 real*8 function electron_pxc_rho(rho) 813 implicit none 814 real*8 rho(*) 815 816#include "bafdecls.fh" 817#include "electron_common.fh" 818 819 820* **** local variables **** 821 integer nx,ny,nz 822 real*8 pxc,pxc2,dv 823 824* **** external functions ***** 825 real*8 lattice_omega 826 external lattice_omega 827 828 call D3dB_nx(1,nx) 829 call D3dB_ny(1,ny) 830 call D3dB_nz(1,nz) 831 832 dv = lattice_omega()/dble(nx*ny*nz) 833 834* **** get exchange-correlation energy **** 835 call D3dB_rr_dot(1,rho, 836 > dbl_mb(xcp(1)), 837 > pxc) 838 if (ispin.eq.1) then 839 pxc = pxc + pxc 840 else 841 call D3dB_rr_dot(1,rho, 842 > dbl_mb(xcp(1)+2*nfft3d), 843 > pxc2) 844 pxc = (pxc + pxc2) 845 end if 846 pxc = pxc*dv 847 848 electron_pxc_rho = pxc 849 return 850 end 851 852 853* *********************************** 854* * * 855* * electron_xcp_ptr * 856* * * 857* *********************************** 858 integer function electron_xcp_ptr() 859 implicit none 860 861#include "electron_common.fh" 862 863 electron_xcp_ptr = xcp(1) 864 return 865 end 866 867 868* *********************************** 869* * * 870* * electron_SIC_energies * 871* * * 872* *********************************** 873 874 subroutine electron_SIC_energies(ehsic0,phsic0,exsic0,pxsic0) 875 implicit none 876 real*8 ehsic0,phsic0 877 real*8 exsic0,pxsic0 878 879#include "bafdecls.fh" 880#include "electron_common.fh" 881 882 logical pspw_SIC_relaxed 883 external pspw_SIC_relaxed 884 885 if (.not.pspw_SIC_relaxed()) then 886 call pspw_energy_SIC(ispin,dbl_mb(psi_r(1)), 887 > ehsic, 888 > phsic, 889 > exsic, 890 > pxsic) 891 phsic = 0.0d0 892 pxsic = 0.0d0 893 end if 894 895 ehsic0 = ehsic 896 exsic0 = exsic 897 phsic0 = phsic 898 pxsic0 = pxsic 899 return 900 end 901 902 903 904* *********************************** 905* * * 906* * electron_SIC_stress * 907* * * 908* *********************************** 909 subroutine electron_SIC_stress(stress) 910 implicit none 911 real*8 stress(3,3) 912 913#include "bafdecls.fh" 914#include "electron_common.fh" 915 916 917 call pspw_SIC_euv(ispin,dbl_mb(psi_r(1)),stress) 918 return 919 end 920 921 922* *********************************** 923* * * 924* * electron_HFX_stress * 925* * * 926* *********************************** 927 subroutine electron_HFX_stress(stress) 928 implicit none 929 real*8 stress(3,3) 930 931#include "bafdecls.fh" 932#include "electron_common.fh" 933 934 935 call pspw_energy_euv_HFX(ispin,dbl_mb(psi_r(1)),stress) 936 return 937 end 938 939 940 941* *********************************** 942* * * 943* * electron_HFX_energies * 944* * * 945* *********************************** 946 947 subroutine electron_HFX_energies(ehfx0,phfx0) 948 implicit none 949 real*8 ehfx0,phfx0 950 951#include "bafdecls.fh" 952#include "electron_common.fh" 953 954 logical pspw_HFX_relaxed 955 external pspw_HFX_relaxed 956 957 if (.not.pspw_HFX_relaxed()) then 958 call pspw_energy_HFX(ispin,dbl_mb(psi_r(1)), 959 > ehfx, 960 > phfx) 961 phfx = 0.0d0 962 end if 963 964 ehfx0 = ehfx 965 phfx0 = phfx 966 return 967 end 968 969 970* *********************************** 971* * * 972* * electron_U_energies * 973* * * 974* *********************************** 975 976 subroutine electron_U_energies(edftu0,pdftu0) 977 implicit none 978 real*8 edftu0,pdftu0 979 980#include "bafdecls.fh" 981#include "electron_common.fh" 982 983 984 edftu0 = edftu 985 pdftu0 = pdftu 986 return 987 end 988 989 990* *********************************** 991* * * 992* * electron_meta_energies * 993* * * 994* *********************************** 995 996 subroutine electron_meta_energies(emeta0,pmeta0) 997 implicit none 998 real*8 emeta0,pmeta0 999 1000#include "bafdecls.fh" 1001#include "electron_common.fh" 1002 1003 1004 emeta0 = emeta 1005 pmeta0 = pmeta 1006 return 1007 end 1008 1009 1010 1011 1012* *********************************** 1013* * * 1014* * electron_get_Hpsi_k * 1015* * * 1016* *********************************** 1017 subroutine electron_get_Hpsi_k(Hpsi_k_new) 1018 implicit none 1019 complex*16 Hpsi_k_new(*) 1020 1021#include "bafdecls.fh" 1022#include "electron_common.fh" 1023 1024 1025 call dcopy(2*npack1*(neq(1)+neq(2)), 1026 > dcpl_mb(Hpsi_k(1)),1, 1027 > Hpsi_k_new,1) 1028 return 1029 end 1030 1031 1032 1033* *************************** 1034* * * 1035* * electron_ispin * 1036* * * 1037* *************************** 1038 integer function electron_ispin() 1039 implicit none 1040 1041#include "electron_common.fh" 1042 1043 electron_ispin = ispin 1044 return 1045 end 1046 1047 1048* *************************** 1049* * * 1050* * electron_ne * 1051* * * 1052* *************************** 1053 integer function electron_ne(ms) 1054 implicit none 1055 integer ms 1056 1057#include "electron_common.fh" 1058 1059 electron_ne = ne(ms) 1060 return 1061 end 1062 1063* *************************** 1064* * * 1065* * electron_neq * 1066* * * 1067* *************************** 1068 integer function electron_neq(ms) 1069 implicit none 1070 integer ms 1071 1072#include "electron_common.fh" 1073 1074 electron_neq = neq(ms) 1075 return 1076 end 1077 1078 1079* *********************************** 1080* * * 1081* * electron_get_Tgradient * 1082* * * 1083* *********************************** 1084 1085 subroutine electron_get_Tgradient(psi_k,hml,THpsi_k) 1086 implicit none 1087 complex*16 psi_k(*) 1088 real*8 hml(*) 1089 complex*16 THpsi_k(*) 1090 1091#include "bafdecls.fh" 1092#include "electron_common.fh" 1093 1094 1095* ***** local variables **** 1096 integer ms,n,shift1,shift2 1097 1098 call Dneall_fmf_Multiply(0, 1099 > psi_k,npack1, 1100 > hml, 1.0d0, 1101 > THpsi_k,0.0d0) 1102 call DAXPY_OMP(2*npack1*(neq(1)+neq(2)), 1103 > (-1.0d0), 1104 > dcpl_mb(Hpsi_k(1)),1, 1105 > THpsi_k,1) 1106 return 1107 end 1108 1109 1110* *********************************** 1111* * * 1112* * electron_gen_Tangent * 1113* * * 1114* *********************************** 1115 1116 subroutine electron_gen_Tangent(psi_k,hml,THpsi_k) 1117 implicit none 1118 complex*16 psi_k(*) 1119 real*8 hml(*) 1120 complex*16 THpsi_k(*) 1121 1122#include "bafdecls.fh" 1123#include "electron_common.fh" 1124 1125* ***** local variables **** 1126c integer ms,n,shift1,shift2 1127 1128 call Dneall_fmf_Multiply(0,psi_k,npack1, 1129 > hml,1.0d0, 1130 > THpsi_k,-1.0d0) 1131 1132c do ms=1,ispin 1133c n = ne(ms) 1134c if (n.le.0) go to 30 1135c shift1 = 1 + (ms-1)*ne(1)*npack1 1136c shift2 = 1 + (ms-1)*ne(1)*ne(1) 1137c call DGEMM('N','N',2*npack1,n,n, 1138c > (1.0d0), 1139c > psi_k(shift1), 2*npack1, 1140c > hml(shift2), n, 1141c > (-1.0d0), 1142c > THpsi_k(shift1),2*npack1) 1143c 30 continue 1144c end do 1145 1146 return 1147 end 1148 1149 1150* *********************************** 1151* * * 1152* * electron_get_Gradient * 1153* * * 1154* *********************************** 1155 1156 subroutine electron_get_Gradient(THpsi_k) 1157 implicit none 1158 complex*16 THpsi_k(*) 1159 1160#include "bafdecls.fh" 1161#include "electron_common.fh" 1162 1163 1164 call dcopy(2*npack1*(neq(1)+neq(2)), 1165 > dcpl_mb(Hpsi_k(1)),1, 1166 > THpsi_k,1) 1167 return 1168 end 1169 1170 1171* *********************************** 1172* * * 1173* * electron_get_gradient_orb * 1174* * * 1175* *********************************** 1176 1177 subroutine electron_get_gradient_orb(i,Horb) 1178 implicit none 1179 integer i 1180 complex*16 Horb(*) 1181 1182#include "bafdecls.fh" 1183#include "electron_common.fh" 1184 1185 call Pack_c_Copy(1,dcpl_mb(Hpsi_k(1)+(i-1)*npack1),Horb) 1186 1187 return 1188 end 1189 1190 1191 1192* *********************************** 1193* * * 1194* * electron_get_TMgradient * 1195* * * 1196* *********************************** 1197 1198 subroutine electron_get_TMgradient(psi_k,THpsi_k) 1199 implicit none 1200#include "errquit.fh" 1201 complex*16 psi_k(*) 1202 complex*16 THpsi_k(*) 1203 1204#include "bafdecls.fh" 1205#include "electron_common.fh" 1206 1207 1208* ***** local variables **** 1209 logical value 1210 integer ms,n,n1(2),shift,mhml(2) 1211 1212 1213 n1(1) = 1 1214 n1(2) = ne(1)+1 1215 1216 value = BA_push_get(mt_dbl,(2*ne(1)*ne(1)),'mhml',mhml(2),mhml(1)) 1217 if (.not. value) 1218 > call errquit('electron_get_Tgradient: push stack',0, MA_ERR) 1219 1220 1221 1222* **** generate M*H|psi> ***** 1223 call Grsm_gg_Copy(npack1,(neq(1)+neq(2)), 1224 > dcpl_mb(Hpsi_k(1)), 1225 > THpsi_k) 1226 call ke_Precondition(npack1,(neq(1)+neq(2)), 1227 > psi_k, 1228 > THpsi_k) 1229 1230 1231* **** generate mhml = <psi|M*H|psi> **** 1232 do ms=1,ispin 1233 shift = (ms-1)*ne(1)*ne(1) 1234 n = ne(ms) 1235 call Grsm_ggm_dot(npack1,n, 1236 > psi_k(1+(ms-1)*ne(1)*npack1), 1237 > dcpl_mb(Hpsi_k(1)+(ms-1)*ne(1)*npack1), 1238 > dbl_mb(mhml(1)+shift)) 1239 end do 1240 1241* **** mhml = -mhml **** 1242 call DSCAL_OMP(2*ne(1)*ne(1),(-1.0d0),dbl_mb(mhml(1)),1) 1243 1244 1245* **** generate TMG = M*H|psi> - |psi>*mhml **** 1246 do ms=1,ispin 1247 shift = (ms-1)*ne(1)*ne(1) 1248 n = ne(ms) 1249 call Grsm_gmg_daxpy(npack1,n, 1250 > psi_k(1+(ms-1)*ne(1)*npack1), 1251 > dbl_mb(mhml(1)+shift), 1252 > THpsi_k(1+(ms-1)*ne(1)*npack1)) 1253 end do 1254 1255 call Grsm_gg_dScale1(npack1,(neq(1)+neq(2)), 1256 > (-1.0d0), 1257 > THpsi_k) 1258 1259 1260 value = BA_pop_stack(mhml(2)) 1261 if (.not. value) 1262 > call errquit('electron_get_Tradient: popping stack',0, MA_ERR) 1263 1264 return 1265 end 1266 1267 1268* *************************** 1269* * * 1270* * electron_gen_hml * 1271* * * 1272* *************************** 1273 1274 subroutine electron_gen_hml(psi_k,hml) 1275 implicit none 1276 complex*16 psi_k(*) 1277 real*8 hml(*) 1278 1279#include "bafdecls.fh" 1280#include "electron_common.fh" 1281 1282 1283c* **** local variables **** 1284c integer ms,n,n1(2),shift 1285 1286 call Dneall_ffm_sym_Multiply(0,psi_k, 1287 > dcpl_mb(Hpsi_k(1)),npack1, 1288 > hml) 1289 1290c n1(1) = 1 1291c n1(2) = ne(1) + 1 1292 1293c do ms=1,ispin 1294c shift = (ms-1)*ne(1)*ne(1) 1295c n = ne(ms) 1296c if (n.le.0) go to 30 1297cc call Grsm_ggm_sym_dot(npack1,n, 1298cc > psi_k(1+(ms-1)*ne(1)*npack1), 1299cc > dcpl_mb(Hpsi_k(1)+(ms-1)*ne(1)*npack1), 1300cc > hml(shift+1)) 1301c call Pack_ccm_sym_dot(1,n, 1302c > psi_k(1+(ms-1)*ne(1)*npack1), 1303c > dcpl_mb(Hpsi_k(1)+(ms-1)*ne(1)*npack1), 1304c > hml(shift+1)) 1305c 30 continue 1306c end do 1307 1308 return 1309 end 1310 1311 1312 1313 1314* *************************** 1315* * * 1316* * electron_gen_hml_g * 1317* * * 1318* *************************** 1319 1320 subroutine electron_gen_hml_g(psi_k,hml) 1321 implicit none 1322 complex*16 psi_k(*) 1323 real*8 hml(*) 1324 1325#include "bafdecls.fh" 1326#include "electron_common.fh" 1327 1328 1329c* **** local variables **** 1330c integer ms,n,n1(2),shift 1331 1332 call Dneall_ffm_Multiply(0,psi_k, 1333 > dcpl_mb(Hpsi_k(1)),npack1, 1334 > hml) 1335 1336c n1(1) = 1 1337c n1(2) = ne(1) + 1 1338 1339c do ms=1,ispin 1340c n = ne(ms) 1341c if (n.le.0) go to 30 1342c shift = (ms-1)*ne(1)*ne(1) 1343c call Pack_ccm_dot(1,n, 1344c > psi_k(1+(ms-1)*ne(1)*npack1), 1345c > dcpl_mb(Hpsi_k(1)+(ms-1)*ne(1)*npack1), 1346c > hml(shift+1)) 1347c 30 continue 1348c end do 1349 return 1350 end 1351 1352 1353* *********************************** 1354* * * 1355* * electron_gen_hmlt * 1356* * * 1357* *********************************** 1358 1359 subroutine electron_gen_hmlt(psi_k,hmlt) 1360 implicit none 1361 complex*16 psi_k(*) 1362 real*8 hmlt(*) 1363 1364#include "bafdecls.fh" 1365#include "electron_common.fh" 1366 1367 1368c* **** local variables **** 1369c integer ms,n,n1(2),shift 1370 1371 call Dneall_ffm_Multiply(0,dcpl_mb(Hpsi_k(1)), 1372 > psi_k,npack1, 1373 > hmlt) 1374 1375c n1(1) = 1 1376c n1(2) = ne(1) + 1 1377 1378c do ms=1,ispin 1379c n = ne(ms) 1380c if (n.le.0) go to 30 1381c shift = (ms-1)*ne(1)*ne(1) 1382cc call Pack_ccm_sym_dot(1,n, 1383cc > dcpl_mb(Hpsi_k(1)+(ms-1)*ne(1)*npack1), 1384cc > psi_k(1+(ms-1)*ne(1)*npack1), 1385cc > hmlt(shift+1)) 1386c call Pack_ccm_dot(1,n, 1387c > dcpl_mb(Hpsi_k(1)+(ms-1)*ne(1)*npack1), 1388c > psi_k(1+(ms-1)*ne(1)*npack1), 1389c > hmlt(shift+1)) 1390c 30 continue 1391c end do 1392 1393 1394 return 1395 end 1396 1397 1398 1399* ************************************** 1400* * * 1401* * electron_gen_psiTangenthml * 1402* * * 1403* ************************************** 1404 1405 subroutine electron_gen_psiTangenthml(psi_k,THpsi_k,hml) 1406 implicit none 1407 complex*16 psi_k(*) 1408 complex*16 THpsi_k(*) 1409 real*8 hml(*) 1410 1411#include "bafdecls.fh" 1412#include "electron_common.fh" 1413 1414 1415c* **** local variables **** 1416c integer ms,n,shift 1417 1418 call Dneall_ffm_sym_Multiply(0,psi_k, 1419 > THpsi_k,npack1, 1420 > hml) 1421 1422c do ms=1,ispin 1423c n = ne(ms) 1424c if (n.le.0) go to 30 1425c shift = (ms-1)*ne(1)*ne(1) 1426c call Pack_ccm_sym_dot(1,n, 1427c > psi_k(1+(ms-1)*ne(1)*npack1), 1428c > THpsi_k(1+(ms-1)*ne(1)*npack1), 1429c > hml(shift+1)) 1430c 30 continue 1431c end do 1432 1433 return 1434 end 1435 1436 1437 1438 1439 1440 1441************************************************************************** 1442************************************************************************** 1443******* routines below this line are for internal use only ********* 1444************************************************************************** 1445************************************************************************** 1446 1447* *********************************** 1448* * * 1449* * electron_gen_Hpsi_k * 1450* * * 1451* *********************************** 1452 1453 subroutine electron_gen_Hpsi_k(psi_k) 1454 implicit none 1455 complex*16 psi_k(*) 1456 1457#include "bafdecls.fh" 1458#include "electron_common.fh" 1459cccc#include "frac_occ.fh" 1460 1461 1462* **** local variables **** 1463 logical move,fractional 1464 integer n 1465 real*8 fion(3,1) 1466 1467* **** external functions **** 1468 integer control_version 1469 external control_version 1470 1471 move = .false. 1472 fractional = .false. 1473* ****************** 1474* **** get Hpsi **** 1475* ****************** 1476 if (control_version().eq.3) 1477 > call psi_H(ispin,neq,psi_k, 1478 > dbl_mb(psi_r(1)), 1479 > dcpl_mb(vl(1)), 1480 > dbl_mb(v_field(1)),field_exist, 1481 > dcpl_mb(vc(1)), 1482 > dbl_mb(xcp(1)), 1483 > dcpl_mb(Hpsi_k(1)), 1484 > move, 1485 > fion, 1486 > fractional,fion) 1487 1488 if (control_version().eq.4) 1489 > call psi_Hv4(ispin,neq,psi_k, 1490 > dbl_mb(psi_r(1)), 1491 > dcpl_mb(vl(1)), 1492 > dbl_mb(vl_lr(1)), 1493 > dbl_mb(v_field(1)),field_exist, 1494 > dcpl_mb(vc(1)), 1495 > dbl_mb(xcp(1)), 1496 > dcpl_mb(Hpsi_k(1)), 1497 > move, 1498 > fion, 1499 > fractional,fion) 1500 1501 call Grsm_gg_dScale1(npack1,(neq(1)+neq(2)),(-1.0d0), 1502 > dcpl_mb(Hpsi_k(1))) 1503 1504 return 1505 end 1506 1507 1508 1509* *********************************** 1510* * * 1511* * electron_gen_Hpsi_k_orb * 1512* * * 1513* *********************************** 1514 1515 subroutine electron_gen_Hpsi_k_orb(n,psi_k) 1516 implicit none 1517 integer n 1518 complex*16 psi_k(*) 1519 1520#include "bafdecls.fh" 1521#include "electron_common.fh" 1522 1523 1524* **** local variables **** 1525 integer ms,index1,index1r,index2 1526 1527* **** external functions **** 1528 integer control_version 1529 external control_version 1530 1531 1532 if (n.le.neq(1)) then 1533 ms=1 1534 else 1535 ms=2 1536 end if 1537 index1 = (n-1)*nfft3d 1538 index1r = 2*index1 1539 index2 = (n-1)*npack1 1540 1541* ****************** 1542* **** get Hpsi **** 1543* ****************** 1544 if (control_version().eq.3) 1545 > call psi_Horb(.true.,ispin,ms, 1546 > dbl_mb(psi_r(1)), 1547 > dbl_mb(vall(1)), 1548 > psi_k(index2+1), 1549 > dbl_mb( psi_r(1)+index1r), 1550 > dcpl_mb(Hpsi_k(1)+index2)) 1551 1552 if (control_version().eq.4) 1553 > call psi_Horbv4(.true.,ispin,ms, 1554 > dbl_mb(psi_r(1)), 1555 > dbl_mb(vall(1)), 1556 > psi_k(index2+1), 1557 > dbl_mb( psi_r(1)+index1r), 1558 > dcpl_mb(Hpsi_k(1)+index2)) 1559 1560c call Pack_c_SMul(1,(-1.0d0), 1561c > dcpl_mb(Hpsi_k(1)+index2), 1562c > dcpl_mb(Hpsi_k(1)+index2)) 1563 call Pack_c_SMul1(1,(-1.0d0), 1564 > dcpl_mb(Hpsi_k(1)+index2)) 1565 1566 return 1567 end 1568 1569* *********************************** 1570* * * 1571* * electron_get_gradient_virtual * 1572* * * 1573* *********************************** 1574 1575 subroutine electron_get_gradient_virtual(ms,orb,Horb) 1576 implicit none 1577 integer ms 1578 complex*16 orb(*) 1579 complex*16 Horb(*) 1580 1581#include "bafdecls.fh" 1582#include "electron_common.fh" 1583#include "errquit.fh" 1584 1585 1586* **** local variables **** 1587 logical value,hfxon 1588 integer n2ft3d 1589 integer tmp_r(2) 1590 1591* **** external functions **** 1592 integer control_version 1593 external control_version 1594 logical pspw_HFX_virtual 1595 external pspw_HFX_virtual 1596 1597 1598 hfxon = pspw_HFX_virtual() 1599 n2ft3d = 2*nfft3d 1600 1601 value = BA_push_get(mt_dbl,(n2ft3d),'tmp_r',tmp_r(2),tmp_r(1)) 1602 if (.not. value) 1603 > call errquit('electron_get_gradient_virtual: push stack',0, 1604 & MA_ERR) 1605 1606 1607 call Pack_c_Copy(1,orb,dbl_mb(tmp_r(1))) 1608 call Pack_c_unpack(1, dbl_mb(tmp_r(1))) 1609 call D3dB_cr_pfft3b(1,1,dbl_mb(tmp_r(1))) 1610 1611* **** get Hpsi **** 1612 if (control_version().eq.3) 1613 > call psi_Horb_replicated(hfxon,ispin,ms, 1614 > dbl_mb(psi_r(1)), 1615 > dbl_mb(vall(1)), 1616 > orb, 1617 > dbl_mb(tmp_r(1)), 1618 > Horb) 1619 1620 if (control_version().eq.4) 1621 > call psi_Horbv4_replicated(hfxon,ispin,ms, 1622 > dbl_mb(psi_r(1)), 1623 > dbl_mb(vall(1)), 1624 > orb, 1625 > dbl_mb(tmp_r(1)), 1626 > Horb) 1627 1628 call Pack_c_SMul1(1,(-1.0d0),Horb) 1629 1630 value = BA_pop_stack(tmp_r(2)) 1631 if (.not. value) call errquit( 1632 > 'electron_get_gradient_virtual: poping stack',1, MA_ERR) 1633 1634 1635 return 1636 end 1637 1638* *********************************** 1639* * * 1640* * electron_get_H0psi_k_orb * 1641* * * 1642* *********************************** 1643 subroutine electron_get_H0psi_k_orb(orb,Horb) 1644 implicit none 1645 complex*16 orb(*) 1646 complex*16 Horb(*) 1647 1648#include "bafdecls.fh" 1649#include "electron_common.fh" 1650#include "errquit.fh" 1651 1652 1653* **** local variables **** 1654 logical value 1655 integer n2ft3d 1656 integer tmp_r(2) 1657 1658* **** external functions **** 1659 integer control_version 1660 external control_version 1661 1662 n2ft3d = 2*nfft3d 1663 1664 value = BA_push_get(mt_dbl,n2ft3d,'tmp_r',tmp_r(2),tmp_r(1)) 1665 if (.not.value) 1666 > call errquit('electron_get_H0psi_k_orb: push stack',0,MA_ERR) 1667 1668 1669 call Pack_c_Copy(1,orb,dbl_mb(tmp_r(1))) 1670 call Pack_c_unpack(1, dbl_mb(tmp_r(1))) 1671 call D3dB_cr_pfft3b(1,1,dbl_mb(tmp_r(1))) 1672 1673* **** get Hpsi **** 1674 call psi_H0orb_replicated( 1675 > dbl_mb(vall(1)), 1676 > orb, 1677 > dbl_mb(tmp_r(1)), 1678 > Horb) 1679 1680 call Pack_c_SMul1(1,(-1.0d0),Horb) 1681 1682 value = BA_pop_stack(tmp_r(2)) 1683 if (.not.value) call errquit( 1684 > 'electron_get_H0psi_k_orb: poping stack',1,MA_ERR) 1685 1686 return 1687 end 1688 1689 1690 1691* *************************** 1692* * * 1693* * electron_gen_psi_r * 1694* * * 1695* *************************** 1696 1697 subroutine electron_gen_psi_r(psi_k) 1698 implicit none 1699 complex*16 psi_k(*) 1700 1701#include "bafdecls.fh" 1702#include "electron_common.fh" 1703 1704* **** local variables **** 1705 integer n,nemax,n2ft3d 1706 1707 1708* ***** generate compensation charge **** 1709 n2ft3d = 2*nfft3d 1710 nemax = neq(1) + neq(2) 1711 1712c call Grsm_gg_Copy(npack1,nemax,psi_k,dbl_mb(psi_r(1))) 1713!$OMP DO private(n) 1714 do n=1,nemax 1715 call Pack_c_Copy0(1,psi_k(1+(n-1)*npack1), 1716 > dbl_mb(psi_r(1)+(n-1)*n2ft3d)) 1717 end do 1718!$OMP END DO 1719 1720 call Grsm_gh_fftb(nfft3d,nemax,dbl_mb(psi_r(1))) 1721 call Grsm_h_Zero_Ends(nfft3d,nemax,dbl_mb(psi_r(1))) !*** probably not neeeded! 1722 1723* **** generate tau functions **** 1724 call nwpw_meta_gga_gen_tau(ispin,neq,psi_k) 1725 1726 return 1727 end 1728 1729* *********************************** 1730* * * 1731* * electron_gen_psi_r_orb * 1732* * * 1733* *********************************** 1734 1735 subroutine electron_gen_psi_r_orb(n,psi_k) 1736 implicit none 1737 integer n 1738 complex*16 psi_k(*) 1739 1740#include "bafdecls.fh" 1741#include "electron_common.fh" 1742 1743* **** local variables **** 1744 integer n2ft3d 1745 1746 n2ft3d = 2*nfft3d 1747 1748 call Pack_c_Copy(1,psi_k(1+(n-1)*npack1), 1749 > dbl_mb(psi_r(1)+(n-1)*n2ft3d)) 1750 1751 call Pack_c_unpack(1, dbl_mb(psi_r(1)+(n-1)*n2ft3d)) 1752 !call D3dB_cr_fft3b(1, dbl_mb(psi_r(1)+(n-1)*n2ft3d)) 1753 call D3dB_cr_pfft3b(1,1,dbl_mb(psi_r(1)+(n-1)*n2ft3d)) 1754 call D3dB_r_Zero_Ends(1,dbl_mb(psi_r(1)+(n-1)*n2ft3d)) 1755 1756 return 1757 end 1758 1759* *************************** 1760* * * 1761* * electron_gen_density * 1762* * * 1763* *************************** 1764 1765 subroutine electron_gen_density(psi_k,dn,fractional,occ) 1766 implicit none 1767 complex*16 psi_k(*) 1768 real*8 dn(*) 1769 logical fractional 1770 real*8 occ(*) 1771 1772#include "bafdecls.fh" 1773#include "electron_common.fh" 1774cccc#include "frac_occ.fh" 1775 1776 1777* **** local variables **** 1778 integer i 1779 integer ms,n2ft3d 1780 integer n,n1(2),n2(2) 1781 real*8 scal2,wf 1782 integer tmp1(2) 1783 logical value 1784 1785* ***** external functions ***** 1786 logical psp_semicore 1787 real*8 lattice_omega 1788 external psp_semicore 1789 external lattice_omega 1790 1791 1792 n1(1) = 1 1793 n1(2) = neq(1) + 1 1794 n2(1) = neq(1) 1795 n2(2) = neq(1) + neq(2) 1796 1797 n2ft3d = 2*nfft3d 1798 scal2 = 1.0d0/lattice_omega() 1799 1800 1801* ********************* 1802* **** generate dn **** 1803* ********************* 1804 call dcopy(2*n2ft3d,0.0d0,0,dn,1) 1805 if (fractional) then 1806 do ms=1,ispin 1807 do n=n1(ms),n2(ms) 1808 wf = occ(n) 1809 do i=1,n2ft3d 1810 dn(i+(ms-1)*n2ft3d) 1811 > = dn(i+(ms-1)*n2ft3d) 1812 > + wf*scal2*(dbl_mb(psi_r(1)+i-1+(n-1)*n2ft3d)**2) 1813 end do 1814 end do 1815 call D3dB_r_Zero_Ends(1,dn(1+(ms-1)*n2ft3d)) 1816 call D1dB_Vector_SumAll(n2ft3d,dn(1+(ms-1)*n2ft3d)) 1817 end do 1818 1819 else 1820 do ms=1,ispin 1821 do n=n1(ms),n2(ms) 1822 do i=1,n2ft3d 1823 dn(i+(ms-1)*n2ft3d) 1824 > = dn(i+(ms-1)*n2ft3d) 1825 > + scal2*(dbl_mb(psi_r(1)+i-1+(n-1)*n2ft3d)**2) 1826 end do 1827 end do 1828 call D3dB_r_Zero_Ends(1,dn(1+(ms-1)*n2ft3d)) 1829 call D1dB_Vector_SumAll(n2ft3d,dn(1+(ms-1)*n2ft3d)) 1830 end do 1831 end if 1832 1833 return 1834 end 1835 1836 1837 1838* *************************** 1839* * * 1840* * electron_gen_densities * 1841* * * 1842* *************************** 1843 1844 subroutine electron_gen_densities(psi_k,dn,dng,dnall, 1845 > fractional,occ) 1846 implicit none 1847 complex*16 psi_k(*) 1848 real*8 dn(*) 1849 complex*16 dng(*) 1850 real*8 dnall(*) 1851 logical fractional 1852 real*8 occ(*) 1853 1854#include "bafdecls.fh" 1855#include "electron_common.fh" 1856cccccccccc#include "frac_occ.fh" 1857 1858 1859* **** local variables **** 1860 integer i 1861 integer ms,n2ft3d 1862 integer n,n1(2),n2(2) 1863 real*8 scal2,wf 1864 integer tmp1(2) 1865 logical value 1866 1867* ***** external functions ***** 1868 logical psp_semicore 1869 real*8 lattice_omega 1870 external psp_semicore 1871 external lattice_omega 1872 1873 1874 n1(1) = 1 1875 n1(2) = neq(1) + 1 1876 n2(1) = neq(1) 1877 n2(2) = neq(1) + neq(2) 1878 1879 n2ft3d = 2*nfft3d 1880 scal2 = 1.0d0/lattice_omega() 1881 1882 1883* ********************* 1884* **** generate dn **** 1885* ********************* 1886 1887 !call dcopy(2*n2ft3d,0.0d0,0,dn,1) 1888 call Parallel_shared_vector_zero(.true.,2*n2ft3d,dn) 1889 if (fractional) then 1890 do ms=1,ispin 1891 do n=n1(ms),n2(ms) 1892 wf = occ(n) 1893!$OMP DO private(i) 1894 do i=1,n2ft3d 1895 dn(i+(ms-1)*n2ft3d) 1896 > = dn(i+(ms-1)*n2ft3d) 1897 > + wf*scal2*(dbl_mb(psi_r(1)+i-1+(n-1)*n2ft3d)**2) 1898 end do 1899!$OMP END DO 1900 end do 1901 call D3dB_r_Zero_Ends(1,dn(1+(ms-1)*n2ft3d)) !*** probably not needed! 1902 call D1dB_Vector_SumAll(n2ft3d,dn(1+(ms-1)*n2ft3d)) 1903 end do 1904 else 1905 do ms=1,ispin 1906 do n=n1(ms),n2(ms) 1907!$OMP DO private(i) 1908 do i=1,n2ft3d 1909 dn(i+(ms-1)*n2ft3d) 1910 > = dn(i+(ms-1)*n2ft3d) 1911 > + scal2*(dbl_mb(psi_r(1)+i-1+(n-1)*n2ft3d)**2) 1912 end do 1913!$OMP END DO 1914 end do 1915 call D3dB_r_Zero_Ends(1,dn(1+(ms-1)*n2ft3d)) !*** probably not needed! 1916 call D1dB_Vector_SumAll(n2ft3d,dn(1+(ms-1)*n2ft3d)) 1917 end do 1918 end if 1919 1920 1921* **** generate dng and dnall **** 1922 call electron_gen_dng_dnall(dn,dng,dnall) 1923 1924 return 1925 end 1926 1927 1928 1929* *************************** 1930* * * 1931* * electron_gen_dng_dnall * 1932* * * 1933* *************************** 1934 1935 subroutine electron_gen_dng_dnall(dn,dng,dnall) 1936 implicit none 1937 real*8 dn(*) 1938 complex*16 dng(*) 1939 real*8 dnall(*) 1940 1941#include "bafdecls.fh" 1942#include "electron_common.fh" 1943#include "errquit.fh" 1944 1945 1946* **** local variables **** 1947 integer i 1948 integer ms,nx,ny,nz,n2ft3d 1949 integer n,n1(2),n2(2) 1950 real*8 scal1 1951 integer tmp1(2) 1952 logical value 1953 1954* ***** external functions ***** 1955 logical psp_semicore 1956 external psp_semicore 1957 1958 n2ft3d = 2*nfft3d 1959 call D3dB_nx(1,nx) 1960 call D3dB_ny(1,ny) 1961 call D3dB_nz(1,nz) 1962 scal1 = 1.0d0/dble(nx*ny*nz) 1963 1964* ********************** 1965* **** generate dng **** 1966* ********************** 1967 value = BA_push_get(mt_dbl,(n2ft3d),'tmp1',tmp1(2),tmp1(1)) 1968 if (.not. value) call errquit( 1969 > 'electron_gen_dng_dnall: out of stack memory',0, MA_ERR) 1970 1971 call D3dB_rr_Sum(1,dn,dn(1+(ispin-1)*n2ft3d),dbl_mb(tmp1(1))) 1972c call D3dB_r_SMul(1,scal1,dbl_mb(tmp1(1)),dbl_mb(tmp1(1))) 1973 call D3dB_r_SMul1(1,scal1,dbl_mb(tmp1(1))) 1974 !call D3dB_rc_fft3f(1,dbl_mb(tmp1(1))) 1975 call D3dB_rc_pfft3f(1,0,dbl_mb(tmp1(1))) 1976 call Pack_c_pack(0,dbl_mb(tmp1(1))) 1977 call Pack_c_Copy(0,dbl_mb(tmp1(1)),dng) 1978 1979* ******************************************************** 1980* **** generate dnall - used for semicore corrections **** 1981* ******************************************************** 1982 if (psp_semicore(0)) then 1983 call semicore_density(dbl_mb(tmp1(1))) 1984 call D3dB_r_SMul1(1,0.5d0,dbl_mb(tmp1(1))) 1985 else 1986 call D3dB_r_Zero(1,dbl_mb(tmp1(1))) 1987 end if 1988 do ms=1,ispin 1989 call D3dB_rr_Sum(1,dn(1+(ms-1)*n2ft3d), 1990 > dbl_mb(tmp1(1)), 1991 > dnall(1+(ms-1)*n2ft3d)) 1992 end do 1993 1994 1995 value = BA_pop_stack(tmp1(2)) 1996 if (.not. value) call errquit( 1997 > 'electron_gen_dng_dnall: poping stack',1, MA_ERR) 1998 1999 2000 return 2001 end 2002 2003 2004* *********************************** 2005* * * 2006* * electron_gen_vall * 2007* * * 2008* *********************************** 2009 2010 subroutine electron_gen_vall() 2011 implicit none 2012 2013 2014#include "bafdecls.fh" 2015#include "electron_common.fh" 2016 2017* **** local variables **** 2018 integer ms 2019 real*8 scal2 2020 2021* **** external functions **** 2022 integer control_version 2023 real*8 lattice_omega 2024 external control_version 2025 external lattice_omega 2026 2027 2028 scal2 = 1.0d0/lattice_omega() 2029 2030 if (control_version().eq.3) then 2031 2032* **** add up k-space potentials, vall = scal2*vl + vc **** 2033 call Pack_c_SMul(0,scal2,dcpl_mb(vl(1)), 2034 > dbl_mb(vall(1))) 2035 2036c call Pack_cc_Sum(0,dbl_mb(vall(1)), 2037c > dcpl_mb(vc(1)), 2038c > dbl_mb(vall(1))) 2039 call Pack_cc_Sum2(0,dcpl_mb(vc(1)),dbl_mb(vall(1))) 2040 2041* **** fourier transform k-space potentials **** 2042 call Pack_c_unpack(0,dbl_mb(vall(1))) 2043 !call D3dB_cr_fft3b(1,dbl_mb(vall(1))) 2044 call D3dB_cr_pfft3b(1,0,dbl_mb(vall(1))) 2045 2046* **** add v_field to vall **** 2047c if (field_exist) 2048c > call D3dB_rr_Sum(1,dbl_mb(vall(1)), 2049c > dbl_mb(v_field(1)), 2050c > dbl_mb(vall(1))) 2051 if (field_exist) 2052 > call D3dB_rr_Sum2(1,dbl_mb(v_field(1)),dbl_mb(vall(1))) 2053 2054 else 2055 2056* **** add up k-space potentials, vall = scal2*vsr_l **** 2057 call Pack_c_SMul(0,scal2,dcpl_mb(vl(1)), 2058 > dbl_mb(vall(1))) 2059 2060* **** fourier transform k-space potentials **** 2061 call Pack_c_unpack(0,dbl_mb(vall(1))) 2062 !call D3dB_cr_fft3b(1,dbl_mb(vall(1))) 2063 call D3dB_cr_pfft3b(1,0,dbl_mb(vall(1))) 2064 2065* **** add vlr_l, vc and v_field to vall **** 2066c call D3dB_rr_Sum(1,dbl_mb(vall(1)), 2067c > dbl_mb(vl_lr(1)), 2068c > dbl_mb(vall(1))) 2069c call D3dB_rr_Sum(1,dbl_mb(vall(1)), 2070c > dcpl_mb(vc(1)), 2071c > dbl_mb(vall(1))) 2072c if (field_exist) 2073c > call D3dB_rr_Sum(1,dbl_mb(vall(1)), 2074c > dbl_mb(v_field(1)), 2075c > dbl_mb(vall(1))) 2076 call D3dB_rr_Sum2(1,dbl_mb(vl_lr(1)),dbl_mb(vall(1))) 2077 call D3dB_rr_Sum2(1,dcpl_mb(vc(1)),dbl_mb(vall(1))) 2078 if (field_exist) 2079 > call D3dB_rr_Sum2(1,dbl_mb(v_field(1)),dbl_mb(vall(1))) 2080 2081 if (confine_exist) 2082 > call D3dB_rr_Sum2(1,dbl_mb(v_confine(1)),dbl_mb(vall(1))) 2083 2084 end if 2085 2086* **** add xcp to vall **** 2087c do ms=ispin,1,-1 2088c call D3dB_rr_Sum(1,dbl_mb(vall(1)), 2089c > dbl_mb(xcp(1)+(ms-1)*2*nfft3d), 2090c > dbl_mb(vall(1)+(ms-1)*2*nfft3d)) 2091c call D3dB_r_Zero_Ends(1,dbl_mb(vall(1)+(ms-1)*2*nfft3d)) 2092c end do 2093 2094 if (ispin.eq.2) then 2095 call D3dB_rr_Sum(1,dbl_mb(vall(1)), 2096 > dbl_mb(xcp(1) +2*nfft3d), 2097 > dbl_mb(vall(1)+2*nfft3d)) 2098 call D3dB_r_Zero_Ends(1,dbl_mb(vall(1)+2*nfft3d)) 2099 end if 2100 call D3dB_rr_Sum2(1,dbl_mb(xcp(1)),dbl_mb(vall(1))) 2101 call D3dB_r_Zero_Ends(1,dbl_mb(vall(1))) 2102 2103 return 2104 end 2105 2106 2107* *********************************** 2108* * * 2109* * electron_gen_vall0 * 2110* * * 2111* *********************************** 2112 2113* Generates non self-consistent potentials only. 2114 2115 subroutine electron_gen_vall0() 2116 implicit none 2117 2118#include "bafdecls.fh" 2119#include "electron_common.fh" 2120 2121* **** local variables **** 2122 integer ms 2123 real*8 scal2 2124 2125* **** external functions **** 2126 integer control_version 2127 real*8 lattice_omega 2128 external control_version 2129 external lattice_omega 2130 2131 2132 scal2 = 1.0d0/lattice_omega() 2133 2134 2135 if (control_version().eq.3) then 2136* **** add up k-space potentials, vall = scal2*vl **** 2137 call Pack_c_SMul(0,scal2,dcpl_mb(vl(1)), 2138 > dbl_mb(vall(1))) 2139 2140* **** fourier transform k-space potentials **** 2141 call Pack_c_unpack(0,dbl_mb(vall(1))) 2142 call D3dB_cr_pfft3b(1,0,dbl_mb(vall(1))) 2143 2144 if (field_exist) 2145 > call D3dB_rr_Sum2(1,dbl_mb(v_field(1)),dbl_mb(vall(1))) 2146 2147 else 2148* **** add up k-space potentials, vall = scal2*vsr_l **** 2149 call Pack_c_SMul(0,scal2,dcpl_mb(vl(1)), 2150 > dbl_mb(vall(1))) 2151 2152* **** fourier transform k-space potentials **** 2153 call Pack_c_unpack(0,dbl_mb(vall(1))) 2154 call D3dB_cr_pfft3b(1,0,dbl_mb(vall(1))) 2155 2156 call D3dB_rr_Sum2(1,dbl_mb(vl_lr(1)),dbl_mb(vall(1))) 2157 if (field_exist) 2158 > call D3dB_rr_Sum2(1,dbl_mb(v_field(1)),dbl_mb(vall(1))) 2159 2160 if (confine_exist) 2161 > call D3dB_rr_Sum2(1,dbl_mb(v_confine(1)),dbl_mb(vall(1))) 2162 2163 end if 2164 2165 return 2166 end 2167 2168 2169 2170 2171 2172* *********************************** 2173* * * 2174* * electron_add_oep_to_vall * 2175* * * 2176* *********************************** 2177 2178 subroutine electron_add_oep_to_vall(dn) 2179 implicit none 2180#include "errquit.fh" 2181 real*8 dn(*) 2182 2183#include "bafdecls.fh" 2184#include "electron_common.fh" 2185 2186 logical value 2187 integer ms,n2ft3d,v_oep(2) 2188 2189 call D3dB_n2ft3d(1,n2ft3d) 2190 value = BA_push_get(mt_dbl,(2*n2ft3d),'V_OEP',v_oep(2),v_oep(1)) 2191 if (.not. value) 2192 > call errquit('electron_add_oep_to_vall:out of stack memory',0,0) 2193 2194 2195 call pspw_potential_SIC_OEP(ispin,ne, 2196 > dn, 2197 > dbl_mb(psi_r(1)), 2198 > dbl_mb(v_oep(1))) 2199 2200* **** add v_oep to vall **** 2201 do ms=1,ispin 2202c call D3dB_rr_Sum(1,dbl_mb( vall(1)+(ms-1)*n2ft3d), 2203c > dbl_mb(v_oep(1)+(ms-1)*n2ft3d), 2204c > dbl_mb( vall(1)+(ms-1)*n2ft3d)) 2205 call D3dB_rr_Sum2(1,dbl_mb(v_oep(1)+(ms-1)*n2ft3d), 2206 > dbl_mb( vall(1)+(ms-1)*n2ft3d)) 2207 call D3dB_r_Zero_Ends(1,dbl_mb(vall(1)+(ms-1)*n2ft3d)) 2208 end do 2209 2210 value = BA_pop_stack(v_oep(2)) 2211 if (.not. value) 2212 > call errquit( 2213 > 'electron_add_oep_to_vall:popping stack memory',1,0) 2214 2215 return 2216 end 2217 2218 2219 2220* *********************************** 2221* * * 2222* * electron_get_vall * 2223* * * 2224* *********************************** 2225 2226 subroutine electron_get_vall(vall_out) 2227 implicit none 2228 real*8 vall_out(*) 2229 2230#include "bafdecls.fh" 2231#include "electron_common.fh" 2232 2233 call dcopy(4*nfft3d,dbl_mb(vall(1)),1,vall_out,1) 2234 return 2235 end 2236 2237 2238* *********************************** 2239* * * 2240* * electron_set_vall * 2241* * * 2242* *********************************** 2243 2244 subroutine electron_set_vall(vall_in) 2245 implicit none 2246 real*8 vall_in(*) 2247 2248#include "bafdecls.fh" 2249#include "electron_common.fh" 2250 2251 call dcopy(4*nfft3d,vall_in,1,dbl_mb(vall(1)),1) 2252 return 2253 end 2254 2255 2256* *********************************** 2257* * * 2258* * electron_gen_scf_potentials * 2259* * * 2260* *********************************** 2261 2262 subroutine electron_gen_scf_potentials(dn,dng,dnall) 2263 implicit none 2264 real*8 dn(*) 2265 complex*16 dng(*) 2266 real*8 dnall(*) 2267 2268#include "bafdecls.fh" 2269#include "errquit.fh" 2270#include "electron_common.fh" 2271 2272 2273* ***** local variables **** 2274 logical value 2275 integer n2ft3d,gga 2276 integer tmp1(2),tmp2(2) 2277 real*8 gmma 2278 2279* **** external functions **** 2280 integer control_gga,control_version 2281 external control_gga,control_version 2282 real*8 control_attenuation 2283 external control_attenuation 2284 2285 2286 n2ft3d = 2*nfft3d 2287 gga = control_gga() 2288 2289 2290 if (control_version().eq.3) then 2291 call coulomb_v(dng,dcpl_mb(vc(1))) 2292 end if 2293 2294 if (control_version().eq.4) then 2295 value = BA_push_get(mt_dbl,(n2ft3d),'tmp1',tmp1(2),tmp1(1)) 2296 if (.not. value) call errquit( 2297 > 'electron_gen_scf_potentials: out of stack memory',0,MA_ERR) 2298 2299 call D3dB_rr_Sum(1,dn(1),dn(1+(ispin-1)*n2ft3d), 2300 > dbl_mb(tmp1(1))) 2301 call coulomb2_v(dbl_mb(tmp1(1)),dcpl_mb(vc(1))) 2302 2303 value = BA_pop_stack(tmp1(2)) 2304 if (.not. value) call errquit( 2305 > 'electron_gen_scf_potentials: error popping stack memory',0, 2306 > MA_ERR) 2307 end if 2308 2309 2310* **** xc potential **** 2311 call v_bwexc_all(gga,n2ft3d,ispin,dnall, 2312 > dbl_mb(xcp(1)),dbl_mb(xce(1))) 2313 2314 return 2315 end 2316 2317 2318* *********************************** 2319* * * 2320* * electron_zero_scf_potentials * 2321* * * 2322* *********************************** 2323 2324 subroutine electron_zero_scf_potentials() 2325 implicit none 2326 2327#include "bafdecls.fh" 2328#include "errquit.fh" 2329#include "electron_common.fh" 2330 2331* ***** local variables **** 2332 integer n2ft3d 2333 2334* **** external functions **** 2335 integer control_version 2336 external control_version 2337 2338 n2ft3d = 2*nfft3d 2339 2340* **** zero out scf potentials **** 2341 if (control_version().eq.3) then 2342 call Parallel_shared_vector_zero(.true., 2343 > 2*npack0,dcpl_mb(vc(1))) 2344 else 2345 call Parallel_shared_vector_zero(.true.,n2ft3d,dcpl_mb(vc(1))) 2346 end if 2347 2348 call Parallel_shared_vector_zero(.true.,2*n2ft3d,dbl_mb(xce(1))) 2349 call Parallel_shared_vector_zero(.true.,2*n2ft3d,dbl_mb(xcp(1))) 2350 2351 return 2352 end 2353 2354 2355 2356* *********************************** 2357* * * 2358* * electron_gen_vl_potential * 2359* * * 2360* *********************************** 2361 2362 subroutine electron_gen_vl_potential() 2363 implicit none 2364#include "errquit.fh" 2365 2366#include "bafdecls.fh" 2367#include "electron_common.fh" 2368 2369 2370* **** local variables **** 2371 logical move,value 2372 integer n2ft3d 2373 integer tmp1(2) 2374 integer tmp2(2) 2375 integer r_grid(2) 2376 2377* **** external functions ***** 2378 logical pspw_charge_found 2379 integer control_version 2380 external pspw_charge_found 2381 external control_version 2382 2383 value = BA_push_get(mt_dcpl,(nfft3d),'tmp1',tmp1(2),tmp1(1)) 2384 value = value.and. 2385 > BA_push_get(mt_dbl,(3),'tmp2',tmp2(2),tmp2(1)) 2386 if (.not. value) call 2387 > errquit('electron_gen_vl_potential: out of stack memory',0, 2388 & MA_ERR) 2389 2390 move = .false. 2391 call v_local(dcpl_mb(vl(1)), 2392 > move, 2393 > dcpl_mb(tmp1(1)), 2394 > dbl_mb(tmp2(1))) 2395 2396 value = BA_pop_stack(tmp2(2)) 2397 value = value.and. 2398 > BA_pop_stack(tmp1(2)) 2399 if (.not. value) call errquit( 2400 > 'electron_gen_vl_potential: error popping stack memory',0, 2401 & MA_ERR) 2402 2403 2404* **** generate real-space fields **** 2405 if ((control_version().eq.4).or. 2406 > (pspw_charge_found())) then 2407 2408 value = BA_push_get(mt_dbl,(6*nfft3d),'r_grid', 2409 > r_grid(2),r_grid(1)) 2410 2411 call lattice_r_grid(dbl_mb(r_grid(1))) 2412 2413* **** generate long-range psp potential **** 2414 if (control_version().eq.4) then 2415 call v_lr_local(dbl_mb(r_grid(1)), 2416 > dbl_mb(vl_lr(1))) 2417 end if 2418 2419* **** zero out v_field **** 2420 call dcopy(2*nfft3d,0.0d0,0,dbl_mb(v_field(1)),1) 2421 2422 2423* **** generate charge potential **** 2424 if (pspw_charge_found()) then 2425 field_exist = .true. 2426 n2ft3d = 2*nfft3d 2427 call pspw_charge_Generate_V(n2ft3d, 2428 > dbl_mb(r_grid(1)), 2429 > dbl_mb(v_field(1))) 2430 end if 2431 2432 2433 value = BA_pop_stack(r_grid(2)) 2434 if (.not. value) call errquit( 2435 > 'electron_gen_vl_potential: error popping stack memory',0, 2436 & MA_ERR) 2437 2438 end if 2439 2440 2441 2442 return 2443 end 2444 2445 2446 2447* *********************************** 2448* * * 2449* * electron_psi_vl_ave * 2450* * * 2451* *********************************** 2452 2453 real*8 function electron_psi_vl_ave(psi1,dn) 2454 implicit none 2455 complex*16 psi1(*) 2456 real*8 dn(*) 2457 2458#include "bafdecls.fh" 2459#include "electron_common.fh" 2460#include "errquit.fh" 2461 2462 2463* **** local variables **** 2464 logical value 2465 integer n,ms,n1(2),n2(2) 2466 integer nx,ny,nz,n2ft3d,np 2467 real*8 elocal,sum,scal1,scal2,dv 2468 integer tmp1(2),tmp2(2) 2469 2470 common /eelectron_ejtmp/ sum,elocal 2471 2472 2473* **** external functions *** 2474 integer control_version 2475 real*8 lattice_omega 2476 external control_version 2477 external lattice_omega 2478 2479 2480 call Parallel_np(np) 2481 2482 n2ft3d = 2*nfft3d 2483 value = BA_push_get(mt_dbl,(n2ft3d),'tmp1',tmp1(2),tmp1(1)) 2484 value = value.and. 2485 > BA_push_get(mt_dbl,(n2ft3d),'tmp2',tmp2(2),tmp2(1)) 2486 if (.not. value) call errquit( 2487 > 'electron_psi_vl_ave: out of stack memory',0, MA_ERR) 2488 2489 call D3dB_nx(1,nx) 2490 call D3dB_ny(1,ny) 2491 call D3dB_nz(1,nz) 2492 n1(1) = 1 2493 n2(1) = neq(1) 2494 n1(2) = neq(1) + 1 2495 n2(2) = neq(1) + neq(2) 2496 2497 scal1 = 1.0d0/dble(nx*ny*nz) 2498 scal2 = 1.0d0/lattice_omega() 2499 dv = scal1/scal2 2500 2501!$OMP MASTER 2502 elocal = 0.0d0 2503!$OMP END MASTER 2504 2505* **** average Kohn-Sham v_local energy **** 2506 call Pack_c_Copy(0,dcpl_mb(vl(1)),dbl_mb(tmp1(1))) 2507 call Pack_c_unpack(0,dbl_mb(tmp1(1))) 2508 !call D3dB_cr_fft3b(1,dbl_mb(tmp1(1))) 2509 call D3dB_cr_pfft3b(1,0,dbl_mb(tmp1(1))) 2510 do ms=1,ispin 2511 do n=n1(ms),n2(ms) 2512 call D3dB_rr_Mul(1, 2513 > dbl_mb(tmp1(1)), 2514 > dbl_mb(psi_r(1)+(n-1)*n2ft3d), 2515 > dbl_mb(tmp2(1))) 2516 2517c call D3dB_rc_fft3f(1,dbl_mb(tmp2(1))) 2518c call Pack_c_pack(1,dbl_mb(tmp2(1))) 2519c call Pack_cc_dot(1,psi1(1+(n-1)*npack1), 2520c > dbl_mb(tmp2(1)), 2521c > sum) 2522 2523 call D3dB_rr_idot(1, 2524 > dbl_mb(psi_r(1)+(n-1)*n2ft3d), 2525 > dbl_mb(tmp2(1)), 2526 > sum) 2527 2528!$OMP MASTER 2529 elocal = elocal + sum*scal1*scal2 2530!$OMP END MASTER 2531 end do 2532 end do 2533!$OMP BARRIER 2534 if (np.gt.1) call Parallel_SumAll(elocal) 2535!$OMP MASTER 2536 if (ispin.eq.1) elocal = 2.0d0*elocal 2537!$OMP END MASTER 2538 2539* *** add in long range part of psp **** 2540 if (control_version().eq.4) then 2541 call D3dB_rr_dot(1,dn(1),dbl_mb(vl_lr(1)),sum) 2542!$OMP MASTER 2543 elocal = elocal + sum*dv 2544!$OMP END MASTER 2545 call D3dB_rr_dot(1,dn(1+(ispin-1)*n2ft3d), 2546 > dbl_mb(vl_lr(1)),sum) 2547!$OMP MASTER 2548 elocal = elocal + sum*dv 2549!$OMP END MASTER 2550 2551 end if 2552 2553* **** add in other real-space fields **** 2554 if (field_exist) then 2555 call D3dB_rr_dot(1,dn(1),dbl_mb(v_field(1)),sum) 2556!$OMP MASTER 2557 elocal = elocal + sum*dv 2558!$OMP END MASTER 2559 call D3dB_rr_dot(1,dn(1+(ispin-1)*n2ft3d), 2560 > dbl_mb(v_field(1)),sum) 2561!$OMP MASTER 2562 elocal = elocal + sum*dv 2563!$OMP END MASTER 2564 end if 2565!$OMP BARRIER 2566 2567 2568* ***** ncmp*Vl+ncmp_smooth*vlpaw terms **** 2569c if (paw_exist) then 2570c call nwpw_compcharge_gen_dn_cmp2(ispin, 2571c > dbl_mb(tmp1(1)), 2572c > dbl_mb(tmp2(1))) 2573c call Pack_cc_dot(0, 2574c > dbl_mb(tmp1(1)), 2575c > dcpl_mb(vl(1)), 2576c > sum) 2577c elocal = elocal + sum 2578c call Pack_cc_dot(0, 2579c > dbl_mb(tmp2(1)), 2580c > dcpl_mb(vlpaw(1)), 2581c > sum) 2582c elocal = elocal + sum 2583c 2584c if ((control_version().eq.4).or.(field_exist)) then 2585c call Pack_c_unpack(0,dbl_mb(tmp1(1))) 2586c call D3dB_cr_pfft3b(1,0,dbl_mb(tmp1(1))) 2587c call Pack_c_unpack(0,dbl_mb(tmp2(1))) 2588c call D3dB_cr_pfft3b(1,0,dbl_mb(tmp2(1))) 2589c end if 2590c if (control_version().eq.4) then 2591c call D3dB_rr_dot(1, 2592c > dbl_mb(tmp1(1)), 2593c > dbl_mb(vl_lr(1)), 2594c > sum) 2595c elocal = elocal + sum*dv 2596c call D3dB_rr_dot(1, 2597c > dbl_mb(tmp2(1)), 2598c > dbl_mb(vl_lr_paw(1)), 2599c > sum) 2600c elocal = elocal + sum*dv 2601c end if 2602c if (field_exist) then 2603c call D3dB_rr_dot(1,dbl_mb(tmp1(1)), 2604c > dbl_mb(v_field(1)),sum) 2605c elocal = elocal + sum*dv 2606c end if 2607c end if 2608 2609 value = BA_pop_stack(tmp2(2)) 2610 value = value.and. 2611 > BA_pop_stack(tmp1(2)) 2612 if (.not. value) call errquit( 2613 > 'electron_psi_vl_ave: error popping stack memory',0, 2614 & MA_ERR) 2615 2616 electron_psi_vl_ave = elocal 2617 return 2618 end 2619 2620 2621 2622* *********************************** 2623* * * 2624* * electron_psi_vnl_ave * 2625* * * 2626* *********************************** 2627 2628 real*8 function electron_psi_vnl_ave(psi1,fractional,occ) 2629 implicit none 2630 complex*16 psi1(*) 2631 logical fractional 2632 real*8 occ(*) 2633 2634#include "bafdecls.fh" 2635#include "electron_common.fh" 2636#include "errquit.fh" 2637 2638 real*8 E_vnonlocal 2639 external E_vnonlocal 2640 2641c* **** local variables **** 2642c logical value 2643c integer i,n,ms,n1(2),n2(2),np 2644c integer nee(2) 2645c integer n2ft3d 2646c real*8 enlocal,sum 2647c integer tmp1(2),tmp2(2) 2648c 2649c 2650c call Parallel_np(np) 2651c 2652c n2ft3d = 2*nfft3d 2653c value = BA_push_get(mt_dbl,(n2ft3d),'tmp1',tmp1(2),tmp1(1)) 2654c value = value.and. 2655c > BA_push_get(mt_dbl,(n2ft3d),'tmp2',tmp2(2),tmp2(1)) 2656c if (.not. value) call errquit( 2657c > 'electron_psi_vl_ave: out of stack memory',0, MA_ERR) 2658c 2659c n1(1) = 1 2660c n2(1) = neq(1) 2661c n1(2) = neq(1) + 1 2662c n2(2) = neq(1) + neq(2) 2663c 2664c 2665c 2666c* **** average Kohn-Sham v_nonlocal energy **** 2667c nee(1) = 1 2668c nee(2) = 0 2669c enlocal = 0.0d0 2670c do ms=1,ispin 2671c do n=n1(ms),n2(ms) 2672c call dcopy(n2ft3d,0.0d0,0,dbl_mb(tmp1(1)),1) 2673c call v_nonlocal(ispin,nee,psi1(1+(n-1)*npack1), 2674c > dbl_mb(tmp1(1)), 2675c > .false.,dbl_mb(tmp2(1)),fractional,occ) 2676c call Pack_cc_idot(1,psi1(1+(n-1)*npack1), 2677c > dbl_mb(tmp1(1)), 2678c > sum) 2679c if (fractional) then 2680c call Dneall_qton(n,i) 2681c sum=sum*occ(i) 2682c end if 2683c enlocal = enlocal - sum 2684c end do 2685c end do 2686c if (np.gt.1) call Parallel_SumAll(enlocal) 2687c if (ispin.eq.1) enlocal = enlocal+enlocal 2688c 2689c 2690c value = BA_pop_stack(tmp2(2)) 2691c value = value.and. 2692c > BA_pop_stack(tmp1(2)) 2693c if (.not. value) call errquit( 2694c > 'electron_psi_vl_ave: error popping stack memory',0, 2695c & MA_ERR) 2696c 2697c electron_psi_vnl_ave = enlocal 2698 2699 electron_psi_vnl_ave = E_Vnonlocal(ispin,neq,fractional,occ) 2700 return 2701 end 2702 2703* *********************************** 2704* * * 2705* * electron_psi_v_field_ave * 2706* * * 2707* *********************************** 2708 2709 real*8 function electron_psi_v_field_ave(psi1,dn) 2710 implicit none 2711 complex*16 psi1(*) 2712 real*8 dn(*) 2713 2714#include "bafdecls.fh" 2715#include "electron_common.fh" 2716 2717 2718* **** local variables **** 2719 integer nx,ny,nz,n2ft3d 2720 real*8 elocal,sum,scal1,scal2,dv 2721 2722* **** external functions *** 2723 real*8 lattice_omega 2724 external lattice_omega 2725 2726 2727 n2ft3d = 2*nfft3d 2728 2729 call D3dB_nx(1,nx) 2730 call D3dB_ny(1,ny) 2731 call D3dB_nz(1,nz) 2732 2733 scal1 = 1.0d0/dble(nx*ny*nz) 2734 scal2 = 1.0d0/lattice_omega() 2735 dv = scal1/scal2 2736 2737 2738 elocal = 0.0d0 2739 2740* **** add in other real-space fields **** 2741 if (field_exist) then 2742 call D3dB_rr_dot(1,dn(1),dbl_mb(v_field(1)),sum) 2743 elocal = elocal + sum*dv 2744 call D3dB_rr_dot(1,dn(1+(ispin-1)*n2ft3d), 2745 > dbl_mb(v_field(1)),sum) 2746 elocal = elocal + sum*dv 2747 end if 2748 2749 electron_psi_v_field_ave = elocal 2750 return 2751 end 2752 2753 2754 2755* *********************************** 2756* * * 2757* * electron_semicoreforce * 2758* * * 2759* *********************************** 2760 2761 subroutine electron_semicoreforce(fion) 2762 implicit none 2763 real*8 fion(3,*) 2764 2765#include "bafdecls.fh" 2766#include "electron_common.fh" 2767 2768 2769 call semicore_xc_F(ispin,dbl_mb(xcp(1)),fion) 2770 2771 return 2772 end 2773 2774* *********************************** 2775* * * 2776* * electron_create_confine * 2777* * * 2778* *********************************** 2779 subroutine electron_create_confine() 2780 implicit none 2781 2782#include "bafdecls.fh" 2783#include "electron_common.fh" 2784#include "errquit.fh" 2785 2786 !*** local variables **** 2787 logical value 2788 integer n2ft3d,r_grid(2) 2789 2790 n2ft3d = 2*nfft3d 2791 confine_exist = .true. 2792 value = BA_alloc_get(mt_dbl,n2ft3d, 2793 > 'v_confine',v_confine(2),v_confine(1)) 2794 value = value.and.BA_push_get(mt_dbl,(6*nfft3d),'r_grid', 2795 > r_grid(2),r_grid(1)) 2796 if (.not.value) 2797 > call errquit('electron_create_confine: error allocate heap',0, 2798 > MA_ERR) 2799 2800 call Parallel_shared_vector_zero(.true.,n2ft3d, 2801 > dbl_mb(v_confine(1))) 2802 call lattice_r_grid(dbl_mb(r_grid(1))) 2803 call electron_atom_Generate_V(1,n2ft3d, 2804 > dbl_mb(r_grid(1)), 2805 > dbl_mb(v_confine(1))) 2806 value = BA_pop_stack(r_grid(2)) 2807 if (.not. value) call errquit( 2808 > 'electron_create_confine: error popping stack memory',0, 2809 > MA_ERR) 2810 return 2811 end 2812 2813* *********************************** 2814* * * 2815* * electron_destroy_confine * 2816* * * 2817* *********************************** 2818 subroutine electron_destroy_confine() 2819 implicit none 2820 2821#include "bafdecls.fh" 2822#include "electron_common.fh" 2823#include "errquit.fh" 2824 2825 !*** local variables **** 2826 logical value 2827 2828 confine_exist = .false. 2829 if (.not.BA_free_heap(v_confine(2))) 2830 > call errquit('electron_destroy_confine: error freeing heap',0, 2831 > MA_ERR) 2832 return 2833 end 2834 2835 2836 2837* ********************************** 2838* * * 2839* * electron_atom_Generate_V * 2840* * * 2841* ********************************** 2842 2843 subroutine electron_atom_Generate_V(ctype,n2ft3d,rgrid,Vqm) 2844 implicit none 2845 integer ctype 2846 integer n2ft3d 2847 real*8 rgrid(3,*) 2848 real*8 Vqm(*) 2849 2850#include "stdio.fh" 2851 2852* ***** local variables **** 2853 integer taskid,MASTER 2854 parameter(MASTER=0) 2855 2856 integer ii,ia,k 2857 real*8 x1,y1,z1,q1,r,rs,re 2858 real*8 r0,dr0,beta 2859 2860* **** external functions **** 2861 integer ion_katm_qm,ion_nion_qm,ion_nkatm_qm 2862 external ion_katm_qm,ion_nion_qm,ion_nkatm_qm 2863 real*8 ion_rion 2864 external ion_rion 2865 character*4 ion_atom_qm 2866 external ion_atom_qm 2867 2868 if (ion_nion_qm().gt.0) then 2869 call Parallel_taskid(taskid) 2870 2871* **** point charges **** 2872 if (ctype.eq.1) then 2873 if (taskid.eq.MASTER) then 2874 write(luout,*) 2875 write(luout,*) " - Confining virtual space" 2876 write(luout,*) 2877 > " - Parameters: Symbol Ve rs re" 2878 end if 2879 do ia = 1,ion_nkatm_qm() 2880 call control_fetch_confine(ion_atom_qm(ia),q1,rs,re) 2881 if (taskid.eq.MASTER) 2882 > write(luout,'(15x,A7,F8.3,F8.3,F8.3)') 2883 > ion_atom_qm(ia),q1,rs,re 2884 end do 2885 if (taskid.eq.MASTER) write(luout,*) 2886 2887 do ii=1,ion_nion_qm() 2888 ia = ion_katm_qm(ii) 2889 call control_fetch_confine(ion_atom_qm(ia),q1,rs,re) 2890 r0 = 0.5d0*(rs+re) 2891 dr0 = (re-rs) 2892 beta = dlog(999.0d0)/dr0 2893 2894 x1 = ion_rion(1,ii) 2895 y1 = ion_rion(2,ii) 2896 z1 = ion_rion(3,ii) 2897!$OMP DO 2898 do k=1,n2ft3d 2899 r = (rgrid(1,k)-x1)**2 2900 > + (rgrid(2,k)-y1)**2 2901 > + (rgrid(3,k)-z1)**2 2902 r = dsqrt(r) 2903 Vqm(k) = Vqm(k) 2904 > + q1*(1.0d0-1.0d0/(1.0d0+dexp(beta*(r-r0)))) 2905 end do 2906!$OMP END DO 2907 end do 2908 end if 2909 end if 2910 return 2911 end 2912 2913 2914 2915 2916 2917c* *********************************** 2918c* * * 2919c* * electron_dn_cmp_coulomb * 2920c* * * 2921c* *********************************** 2922c 2923c real*8 function electron_dn_cmp_coulomb() 2924c implicit none 2925c#include "bafdecls.fh" 2926c#include "errquit.fh" 2927c#include "electron_common.fh" 2928c 2929c* **** local variables **** 2930c real*8 E 2931c 2932c E = 0.0d0 2933c if (paw_exist) then 2934c call nwpw_compcharge_gen_dn_cmp2(ispin, 2935c > dcpl_mb(dng_cmp(1)), 2936c > dcpl_mb(dng_cmp_smooth(1))) 2937c end if 2938c 2939c electron_dn_cmp_coulomb = E 2940c return 2941c end 2942 2943