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