TABLE OF CONTENTS
- m_wfd/copy_kdata_0D
- m_wfd/copy_kdata_1D
- m_wfd/kdata_free_0D
- m_wfd/kdata_free_1D
- m_wfd/kdata_init
- m_wfd/kdata_t
- m_wfd/kpt_store_t
- m_wfd/spin_store_t
- m_wfd/test_charge
- m_wfd/wave_copy
- m_wfd/wave_free
- m_wfd/wave_init
- m_wfd/wave_t
- m_wfd/wfd_change_ngfft
- m_wfd/wfd_copy_cg
- m_wfd/wfd_dump_errinfo
- m_wfd/wfd_extract_cgblock
- m_wfd/wfd_free
- m_wfd/wfd_get_cprj
- m_wfd/wfd_get_many_ur
- m_wfd/wfd_get_socpert
- m_wfd/wfd_get_ug
- m_wfd/wfd_get_ur
- m_wfd/wfd_get_wave_prt
- m_wfd/wfd_ihave_ug
- m_wfd/wfd_init
- m_wfd/wfd_mybands
- m_wfd/wfd_norm2
- m_wfd/wfd_paw_get_aeur
- m_wfd/wfd_print
- m_wfd/wfd_push_ug
- m_wfd/wfd_read_wfk
- m_wfd/wfd_reset_ur_cprj
- m_wfd/wfd_sym_ug_kg
- m_wfd/wfd_sym_ur
- m_wfd/wfd_t
- m_wfd/wfd_test_ortho
- m_wfd/wfd_ug2cprj
- m_wfd/wfd_wave_free
- m_wfd/wfd_xdotc
- m_wfd/wfdgw_bands_of_rank
- m_wfd/wfdgw_bks_distrb
- m_wfd/wfdgw_copy
- m_wfd/wfdgw_distribute_bands
- m_wfd/wfdgw_distribute_bbp
- m_wfd/wfdgw_distribute_kb_kpbp
- m_wfd/wfdgw_get_nl_me
- m_wfd/wfdgw_iterator_bks
- m_wfd/wfdgw_mkrho
- m_wfd/wfdgw_pawrhoij
- m_wfd/wfdgw_plot_ur
- m_wfd/wfdgw_rank_has_ug
- m_wfd/wfdgw_rotate
- m_wfd/wfdgw_sanity_check
- m_wfd/wfdgw_show_bkstab
- m_wfd/wfdgw_t
- m_wfd/wfdgw_update_bkstab
- m_wfd/wfdgw_who_has_ug
- m_wfd/wfdgw_write_wfk
m_wfd/copy_kdata_0D [ Functions ]
NAME
copy_kdata_0D
FUNCTION
Copy object
SOURCE
745 subroutine copy_kdata_0D(Kdata_in, Kdata_out) 746 747 !Arguments ------------------------------------ 748 class(kdata_t),intent(in) :: Kdata_in 749 class(kdata_t),intent(inout) :: Kdata_out 750 751 !************************************************************************ 752 753 !@kdata_t 754 Kdata_out%istwfk = Kdata_in%istwfk 755 Kdata_out%npw = Kdata_in%npw 756 Kdata_out%useylm = Kdata_in%useylm 757 Kdata_out%has_ylm = Kdata_in%has_ylm 758 759 call alloc_copy(Kdata_in%kg_k, Kdata_out%kg_k) 760 call alloc_copy(Kdata_in%gbound, Kdata_out%gbound) 761 762 call alloc_copy(Kdata_in%ph3d,Kdata_out%ph3d) 763 call alloc_copy(Kdata_in%phkxred,Kdata_out%phkxred) 764 call alloc_copy(Kdata_in%fnl_dir0der0,Kdata_out%fnl_dir0der0) 765 call alloc_copy(Kdata_in%ylm,Kdata_out%ylm) 766 767 end subroutine copy_kdata_0D
m_wfd/copy_kdata_1D [ Functions ]
NAME
copy_kdata_1D
FUNCTION
Deallocate memory.
SOURCE
781 subroutine copy_kdata_1D(Kdata_in, Kdata_out) 782 783 !Arguments ------------------------------------ 784 !scalars 785 type(kdata_t),intent(in) :: Kdata_in(:) 786 type(kdata_t),intent(inout) :: Kdata_out(:) 787 788 !Local variables ------------------------------ 789 !scalars 790 integer :: ik 791 792 !************************************************************************ 793 794 if (size(Kdata_in,DIM=1) /= size(Kdata_out,DIM=1)) then 795 ABI_ERROR("copy_kdata_1D: wrong sizes !") 796 end if 797 798 do ik=LBOUND(Kdata_in,DIM=1),UBOUND(Kdata_in,DIM=1) 799 call copy_kdata_0d(Kdata_in(ik), Kdata_out(ik)) 800 end do 801 802 end subroutine copy_kdata_1D
m_wfd/kdata_free_0D [ Functions ]
NAME
kdata_free_0D
FUNCTION
Deallocate memory
SOURCE
686 subroutine kdata_free_0D(Kdata) 687 688 !Arguments ------------------------------------ 689 class(kdata_t),intent(inout) :: Kdata 690 691 !************************************************************************ 692 693 ABI_SFREE(Kdata%kg_k) 694 ABI_SFREE(Kdata%gbound) 695 696 ABI_SFREE(Kdata%ph3d) 697 ABI_SFREE(Kdata%phkxred) 698 ABI_SFREE(Kdata%fnl_dir0der0) 699 ABI_SFREE(Kdata%ylm) 700 701 end subroutine kdata_free_0D
m_wfd/kdata_free_1D [ Functions ]
NAME
kdata_free_1D
FUNCTION
Deallocate memory.
SOURCE
715 subroutine kdata_free_1D(Kdata) 716 717 !Arguments ------------------------------------ 718 !scalars 719 type(kdata_t),intent(inout) :: Kdata(:) 720 721 !Local variables ------------------------------ 722 !scalars 723 integer :: ik 724 725 !************************************************************************ 726 727 do ik=LBOUND(Kdata,DIM=1),UBOUND(Kdata,DIM=1) 728 call kdata_free_0D(Kdata(ik)) 729 end do 730 731 end subroutine kdata_free_1D
m_wfd/kdata_init [ Functions ]
NAME
kdata_init
FUNCTION
Main creation method for the kdata_t datatype.
SOURCE
553 subroutine kdata_init(Kdata, Cryst, Psps, kpoint, istwfk, ngfft, MPI_enreg, ecut, kg_k) 554 555 !Arguments ------------------------------------ 556 !scalars 557 class(kdata_t),intent(inout) :: Kdata 558 integer,intent(in) :: istwfk 559 real(dp),optional,intent(in) :: ecut 560 type(crystal_t),intent(in) :: Cryst 561 type(pseudopotential_type),intent(in) :: Psps 562 type(MPI_type),intent(in) :: MPI_enreg 563 !arrays 564 integer,optional,target,intent(in) :: kg_k(:,:) 565 integer,intent(in) :: ngfft(18) 566 real(dp),intent(in) :: kpoint(3) 567 568 !Local variables ------------------------------ 569 !scalars 570 integer,parameter :: ider0 = 0, idir0 = 0 571 integer :: mpw_, npw_k, dimffnl, useylmgr, nkpg, iatom, mkmem_, nkpt_, optder, mgfft, iatm, matblk 572 real(dp) :: arg 573 !arrays 574 integer :: nband_(1), npwarr_(1) 575 real(dp),allocatable :: ylmgr_k(:,:,:),kpg_k(:,:),ph1d(:,:) 576 577 !************************************************************************ 578 579 !@kdata_t 580 Kdata%istwfk = istwfk 581 Kdata%useylm = Psps%useylm 582 583 if (present(ecut)) then 584 ! Calculate G-sphere from input ecut. 585 ABI_CHECK(.not.allocated(Kdata%kg_k), "Kdata%kg_k is allocated!") 586 call get_kg(kpoint,istwfk,ecut,Cryst%gmet,npw_k,Kdata%kg_k) 587 588 else if (present(kg_k)) then 589 ! Use input g-vectors. 590 npw_k = SIZE(kg_k,DIM=2) 591 ABI_MALLOC(Kdata%kg_k,(3,npw_k)) 592 Kdata%kg_k = kg_k 593 else 594 ABI_ERROR("Either ecut or kg_k must be present") 595 end if 596 Kdata%npw = npw_k 597 598 mgfft = MAXVAL(ngfft(1:3)) 599 600 ! Finds the boundary of the basis sphere of G vectors (for this k point) 601 ! for use in improved zero padding of ffts in 3 dimensions. 602 ABI_MALLOC(Kdata%gbound,(2*mgfft+8, 2)) 603 call sphereboundary(Kdata%gbound, istwfk, Kdata%kg_k, mgfft, npw_k) 604 605 ! Compute e^{ik.Ra} for each atom. Packed according to the atom type (atindx). 606 ABI_MALLOC(Kdata%phkxred,(2, Cryst%natom)) 607 do iatom=1,Cryst%natom 608 iatm=Cryst%atindx(iatom) 609 arg=two_pi*(DOT_PRODUCT(kpoint,Cryst%xred(:,iatom))) 610 Kdata%phkxred(1,iatm)=DCOS(arg) 611 Kdata%phkxred(2,iatm)=DSIN(arg) 612 end do 613 614 ! TODO: Should avoid storing all this stuff in memory (risky if lots of k-points) 615 ! Write method to prepare kdata inside loop 616 617 ! Calculate 1-dim structure factor phase information. 618 mgfft = MAXVAL(ngfft(1:3)) 619 ABI_MALLOC(ph1d,(2, 3*(2*mgfft+1)*Cryst%natom)) 620 call getph(Cryst%atindx,Cryst%natom,ngfft(1),ngfft(2),ngfft(3),ph1d,Cryst%xred) 621 622 ! Calculate 3-dim structure factor phase information. 623 matblk = 0 624 if(psps%usepaw == 1 .or. Kdata%use_fnl_dir0der0) then 625 matblk=Cryst%natom 626 end if 627 ABI_MALLOC(Kdata%ph3d,(2, npw_k, matblk)) 628 if(psps%usepaw == 1 .or. Kdata%use_fnl_dir0der0) then 629 call ph1d3d(1,Cryst%natom,Kdata%kg_k,matblk,Cryst%natom,npw_k,ngfft(1),ngfft(2),& 630 &ngfft(3),Kdata%phkxred,ph1d,Kdata%ph3d) 631 end if 632 ABI_FREE(ph1d) 633 634 ! Compute spherical harmonics. 635 Kdata%has_ylm = 1 636 ABI_MALLOC(Kdata%ylm, (npw_k, Psps%mpsang**2*Psps%useylm)) 637 useylmgr=0 638 ABI_MALLOC(ylmgr_k,(npw_k, 3, Psps%mpsang**2*useylmgr)) 639 640 if (Kdata%useylm == 1) then 641 mkmem_=1; mpw_=npw_k; nband_=0; nkpt_=1; npwarr_(1)=npw_k 642 optder=0 ! only Ylm(K) are computed. 643 644 call initylmg(Cryst%gprimd, Kdata%kg_k, kpoint, mkmem_, MPI_enreg, Psps%mpsang, mpw_, nband_, nkpt_,& 645 npwarr_, 1, optder, Cryst%rprimd, Kdata%ylm, ylmgr_k) 646 647 Kdata%has_ylm = 2 648 end if 649 650 ! Compute (k+G) vectors. 651 nkpg = 0 652 ABI_MALLOC(kpg_k,(npw_k, nkpg)) 653 if (nkpg>0) call mkkpg(Kdata%kg_k, kpg_k, kpoint, nkpg, npw_k) 654 655 ! Compute nonlocal form factors fnl_dir0der0 for all (k+G). 656 dimffnl = 0 657 if(psps%usepaw == 1 .or. Kdata%use_fnl_dir0der0) then 658 dimffnl = 1+3*ider0 659 end if 660 ABI_MALLOC(Kdata%fnl_dir0der0,(npw_k, dimffnl, Psps%lmnmax, Cryst%ntypat)) 661 662 if (dimffnl /=0 ) then 663 call mkffnl(Psps%dimekb,dimffnl,Psps%ekb,Kdata%fnl_dir0der0,Psps%ffspl,& 664 Cryst%gmet,Cryst%gprimd,ider0,idir0,Psps%indlmn,Kdata%kg_k,kpg_k,kpoint,Psps%lmnmax,& 665 Psps%lnmax,Psps%mpsang,Psps%mqgrid_ff,nkpg,npw_k,Cryst%ntypat,& 666 Psps%pspso,Psps%qgrid_ff,Cryst%rmet,Psps%usepaw,Psps%useylm,Kdata%ylm,ylmgr_k) 667 end if 668 669 ABI_FREE(kpg_k) 670 ABI_FREE(ylmgr_k) 671 672 end subroutine kdata_init
m_wfd/kdata_t [ Types ]
NAME
kdata_t
FUNCTION
Datatype storing k-dependent quantities and tables needed for performing the zero-padded FFT of wavefunctions.
SOURCE
98 type,public :: kdata_t 99 100 logical :: use_fnl_dir0der0 = .False. 101 ! Decide if we need to use fnl_dir0der0. 102 103 integer :: istwfk = -1 104 ! Storage mode for this k point. 105 106 integer :: npw = -1 107 ! Number of plane-waves for this k-point. 108 109 integer :: useylm = -1 110 ! 1 if nonlocal part is applied using real spherical Harmonics. 0 for Legendre polynomial. 111 112 integer :: has_ylm = -1 113 ! 0 if ylm is not used. 114 ! 1 if ylm is allocated. 115 ! 2 if ylm is already computed. 116 117 integer,allocatable :: kg_k(:,:) 118 ! kg_k(3,npw) 119 ! G vector coordinates in reduced cordinates. 120 121 integer,allocatable :: gbound(:,:) 122 ! gbound(2*mgfft+8,2)) 123 ! The boundary of the basis sphere of G vectors at a given k point. 124 ! for use in improved zero padding of FFTs in 3 dimensions. 125 126 real(dp),allocatable :: ph3d(:,:,:) 127 ! ph3d(2, npw, natom) 128 ! 3-dim structure factors, for each atom and each plane wave. 129 ! Available only for PAW or use_fnl_dir0der0 is true. 130 131 real(dp),allocatable :: phkxred(:,:) 132 ! phkxred(2,natom)) 133 ! e^{ik.Ra} for each atom. Packed according to the atom type (atindx). 134 135 real(dp),allocatable :: fnl_dir0der0(:,:,:,:) 136 ! fnl_dir0der0(npw, 1, lmnmax,ntypat) 137 ! nonlocal form factors. Computed only if usepaw == 1 or use_fnl_dir0der0 is true. 138 ! fnl(k+G).ylm(k+G) if PAW 139 ! f_ln(k+G)/|k+G|^l if NC 140 141 real(dp),allocatable :: ylm(:,:) 142 ! ylm(npw, mpsang**2*useylm) 143 ! Real spherical harmonics for each k+G 144 145 contains 146 147 procedure :: init => kdata_init 148 ! Init object 149 150 procedure :: free => kdata_free_0D 151 ! Free memory 152 153 end type kdata_t 154 155 interface kdata_free 156 module procedure kdata_free_0D 157 module procedure kdata_free_1D 158 end interface kdata_free 159 160 public :: kdata_copy 161 162 interface kdata_copy 163 module procedure copy_kdata_0D 164 module procedure copy_kdata_1D 165 end interface kdata_copy
m_wfd/kpt_store_t [ Types ]
NAME
kpt_store_t
FUNCTION
Used to build ragged arrays of wave_t in compact form.
SOURCE
237 type :: kpt_store_t 238 type(wave_t),allocatable :: b(:) 239 end type kpt_store_t
m_wfd/spin_store_t [ Types ]
NAME
spin_store_t
FUNCTION
Used to build ragged arrays of wave_t in compact form.
SOURCE
253 type :: spin_store_t 254 type(kpt_store_t),allocatable :: k(:) 255 end type spin_store_t
m_wfd/test_charge [ Functions ]
NAME
test_charge
FUNCTION
Reports info on the electronic charge as well as Drude plasma frequency. Mainly used in the GW part.
INPUTS
nelectron_exp=Expected total number of electrons (used to normalize the charge)
OUTPUT
SOURCE
5690 subroutine test_charge(nfftf,nelectron_exp,nspden,rhor,ucvol,& 5691 & usepaw,usexcnhat,usefinegrid,compch_sph,compch_fft,omegaplasma) 5692 5693 !Arguments ------------------------------------ 5694 !scalars 5695 integer,intent(in) :: nfftf,nspden,usefinegrid,usepaw,usexcnhat 5696 real(dp),intent(in) :: compch_fft,compch_sph,ucvol,nelectron_exp 5697 real(dp),intent(out) :: omegaplasma 5698 !arrays 5699 real(dp),intent(inout) :: rhor(nfftf,nspden) 5700 5701 !Local variables ------------------------------ 5702 !scalars 5703 real(dp) :: nelectron_tot,nelectron_fft 5704 real(dp) :: nelectron_pw,nelectron_sph,rhoav,rs,nratio 5705 character(len=500) :: msg 5706 5707 !************************************************************************* 5708 5709 ! ABI_UNUSED(usexcnhat) 5710 if (usexcnhat==0)then 5711 end if 5712 5713 ! === For PAW output of compensation charges === 5714 if (usepaw==1) then 5715 !if (usepaw==1.and.usexcnhat>0) then ! TODO I still dont understand this if! 5716 write(msg,'(4a)')ch10,' PAW TEST:',ch10,' ==== Compensation charge inside spheres ============' 5717 if (compch_sph<greatest_real.and.compch_fft<greatest_real) & 5718 & write(msg,'(3a)')TRIM(msg),ch10,' The following values must be close...' 5719 if (compch_sph<greatest_real) & 5720 & write(msg,'(3a,f22.15)')TRIM(msg),ch10,' Compensation charge over spherical meshes = ',compch_sph 5721 if (compch_fft<greatest_real) then 5722 if (usefinegrid==1) then 5723 write(msg,'(3a,f22.15)')TRIM(msg),ch10,' Compensation charge over fine fft grid = ',compch_fft 5724 else 5725 write(msg,'(3a,f22.15)')TRIM(msg),ch10,' Compensation charge over fft grid = ',compch_fft 5726 end if 5727 end if 5728 call wrtout([std_out, ab_out], msg) 5729 write(msg,'(a)')ch10 5730 call wrtout([std_out, ab_out], msg) 5731 end if !PAW 5732 5733 nelectron_pw =SUM(rhor(:,1))*ucvol/nfftf 5734 nelectron_tot=nelectron_pw 5735 nratio =nelectron_exp/nelectron_tot 5736 5737 if (usepaw==1) then 5738 nelectron_sph=nelectron_pw+compch_sph 5739 nelectron_fft=nelectron_pw+compch_fft 5740 nelectron_tot=nelectron_sph 5741 nratio=(nelectron_exp-nelectron_sph)/nelectron_pw 5742 end if 5743 5744 rhoav=nelectron_tot/ucvol ; rs=(three/(four_pi*rhoav))**third 5745 if (usepaw==0) then 5746 write(msg,'(2(a,f9.4))')& 5747 ' Number of electrons calculated from density = ',nelectron_tot,'; Expected = ',nelectron_exp 5748 else 5749 write(msg,'(2(a,f9.4),a)')& 5750 ' Total number of electrons per unit cell = ',nelectron_sph,' (Spherical mesh), ',nelectron_fft,' (FFT mesh)' 5751 end if 5752 call wrtout([std_out, ab_out], msg) 5753 5754 !$write(msg,'(a,f9.4)')' Renormalizing smooth charge density using nratio = ',nratio 5755 !! rhor(:,:)=nratio*rhor(:,:) 5756 5757 write(msg,'(a,f9.6)')' average of density, n = ',rhoav 5758 call wrtout([std_out, ab_out], msg) 5759 write(msg,'(a,f9.4)')' r_s = ',rs 5760 call wrtout([std_out, ab_out], msg) 5761 omegaplasma=SQRT(four_pi*rhoav) 5762 write(msg,'(a,f9.4,2a)')' omega_plasma = ',omegaplasma*Ha_eV,' [eV]',ch10 5763 call wrtout([std_out, ab_out], msg) 5764 5765 end subroutine test_charge
m_wfd/wave_copy [ Functions ]
NAME
wave_copy
FUNCTION
Copy method for the wave_t datatype.
SOURCE
2057 type(wave_t) function wave_copy(Wave_in) result(Wave_out) 2058 2059 !Arguments ------------------------------------ 2060 !scalars 2061 class(wave_t),intent(in) :: Wave_in 2062 2063 !Local variables ------------------------------ 2064 integer :: natom,nspinor,iatom,ispinor 2065 2066 !************************************************************************ 2067 2068 Wave_out%has_ug = Wave_in%has_ug 2069 Wave_out%has_ur = Wave_in%has_ur 2070 Wave_out%has_cprj = Wave_in%has_cprj 2071 Wave_out%cprj_order = Wave_in%cprj_order 2072 2073 ABI_MALLOC(Wave_out%ug, (SIZE(Wave_in%ug))) 2074 Wave_out%ug = Wave_in%ug 2075 ABI_MALLOC(Wave_out%ur, (SIZE(Wave_in%ur))) 2076 Wave_out%ur = Wave_in%ur 2077 2078 natom = size(Wave_in%Cprj,dim=1) 2079 nspinor = size(Wave_in%Cprj,dim=2) 2080 if ((size(Wave_out%Cprj,dim=1) .ne. natom) .or. (size(Wave_out%Cprj,dim=2) .ne. nspinor)) then 2081 if (allocated(Wave_out%Cprj)) then 2082 ABI_FREE(Wave_out%Cprj) 2083 end if 2084 ABI_MALLOC(Wave_out%Cprj,(natom,nspinor)) 2085 end if 2086 2087 do ispinor=1,nspinor 2088 do iatom=1,natom 2089 Wave_out%Cprj(iatom,ispinor)%ncpgr=Wave_in%Cprj(iatom,ispinor)%ncpgr 2090 Wave_out%Cprj(iatom,ispinor)%nlmn=Wave_in%Cprj(iatom,ispinor)%nlmn 2091 call alloc_copy(Wave_in%Cprj(iatom,ispinor)%cp,Wave_out%Cprj(iatom,ispinor)%cp) 2092 call alloc_copy(Wave_in%Cprj(iatom,ispinor)%dcp,Wave_out%Cprj(iatom,ispinor)%dcp) 2093 end do 2094 end do 2095 2096 end function wave_copy
m_wfd/wave_free [ Functions ]
NAME
wave_free
FUNCTION
Main destruction method for the wave_t datatype.
INPUTS
[what]=String defining what has to be freed.
"A" =Both ug and ur and Cprj. Default.
"G" =Only ug.
"R" =Only ur
"C" =Only PAW Cprj.
SIDE EFFECTS
Memory in Wave is deallocated depending on what
SOURCE
2006 subroutine wave_free(Wave, what) 2007 2008 !Arguments ------------------------------------ 2009 !scalars 2010 class(wave_t),intent(inout) :: Wave 2011 character(len=*),optional,intent(in) :: what 2012 2013 !Local variables ------------------------------ 2014 !scalars 2015 character(len=10) :: my_what 2016 2017 !************************************************************************ 2018 2019 my_what="ALL"; if (present(what)) my_what=toupper(what) 2020 2021 if (.not.firstchar(my_what, ["A", "G", "R", "C"] )) then 2022 ABI_ERROR(sjoin("Unknow what:", what)) 2023 end if 2024 2025 if (firstchar(my_what, ["A", "G"])) then 2026 ABI_SFREE(Wave%ug) 2027 Wave%has_ug = WFD_NOWAVE 2028 end if 2029 2030 if (firstchar(my_what, ["A", "R"])) then 2031 ABI_SFREE(Wave%ur) 2032 Wave%has_ur = WFD_NOWAVE 2033 end if 2034 2035 if (firstchar(my_what, ["A", "C"])) then 2036 if (allocated(Wave%Cprj)) then 2037 call pawcprj_free(Wave%Cprj) 2038 ABI_FREE(Wave%Cprj) 2039 end if 2040 Wave%has_cprj = WFD_NOWAVE 2041 end if 2042 2043 end subroutine wave_free
m_wfd/wave_init [ Functions ]
NAME
wave_init
FUNCTION
Main creation method for the wave_t data type
INPUTS
usepaw=1 if PAW is used. npw =Number of plane-waves for ug nfft=Number of FFT points for the real space wavefunction. nspinor=Number of spinor components. natom=Number of atoms in cprj matrix elements. nlmn_size(natom)=Number of (n,l,m) channel for each atom. Ordering of atoms depends on cprj_order cprj_order=Flag defining the ordering of the atoms in the cprj matrix elements (CPR_RANDOM|CPR_SORTED). Use to know if we have to reorder the matrix elements when wfd_get_cprj is called.
OUTPUT
Wave<wave_t>=The structure fully initialized.
SOURCE
1947 subroutine wave_init(Wave, usepaw, npw, nfft, nspinor, natom, nlmn_size, cprj_order) 1948 1949 !Arguments ------------------------------------ 1950 !scalars 1951 integer,intent(in) :: npw,nfft,nspinor,usepaw,natom 1952 integer(c_int8_t),intent(in) :: cprj_order 1953 type(wave_t),intent(inout) :: Wave 1954 !arrays 1955 integer,intent(in) :: nlmn_size(:) 1956 1957 !Local variables ------------------------------ 1958 !scalars 1959 integer,parameter :: ncpgr0=0 ! For the time being, no derivatives 1960 1961 !************************************************************************ 1962 1963 !@wave_t 1964 if (npw >0) then 1965 ABI_MALLOC(Wave%ug, (npw*nspinor)) 1966 Wave%has_ug = WFD_ALLOCATED 1967 Wave%ug = huge(one_gw) 1968 if (usepaw == 1) then 1969 ABI_MALLOC(Wave%Cprj, (natom,nspinor)) 1970 call pawcprj_alloc(Wave%Cprj,ncpgr0,nlmn_size) 1971 Wave%has_cprj = WFD_ALLOCATED 1972 Wave%cprj_order = cprj_order 1973 end if 1974 end if 1975 1976 if (nfft > 0) then 1977 ABI_MALLOC(Wave%ur, (nfft*nspinor)) 1978 Wave%ur = huge(one_gw) 1979 Wave%has_ur = WFD_ALLOCATED 1980 end if 1981 1982 end subroutine wave_init
m_wfd/wave_t [ Types ]
NAME
wave_t
FUNCTION
Object storing a single wavefunction in G-space and, optionally, its r-space representation.
SOURCE
179 type, public :: wave_t 180 181 !! integer :: cplex 182 ! 1 for real wavefunctions u(r) 183 ! 2 for complex wavefunctions u(r). 184 ! At gamma we always have real u(r) provided that time-reversal can be used. 185 ! In systems with both time-reversal and spatial inversion, wavefunctions can be chosen to be real. 186 ! One might use this to reduce memory in wave_t. 187 188 integer(c_int8_t) :: has_ug = WFD_NOWAVE 189 ! Flag giving the status of ug. 190 191 integer(c_int8_t) :: has_ur = WFD_NOWAVE 192 ! Flag giving the status of ur. 193 194 integer(c_int8_t) :: has_cprj = WFD_NOWAVE 195 ! Flag giving the status of cprj. 196 197 integer(c_int8_t) :: cprj_order = CPR_RANDOM 198 ! Flag defining whether cprj are sorted by atom type or ordered according 199 ! to the typat variable used in the input file. 200 201 complex(gwpc),allocatable :: ug(:) 202 ! ug(npw_k*nspinor) 203 ! The periodic part of the Bloch wavefunction in G-space. 204 205 complex(gwpc),allocatable :: ur(:) 206 ! ur(nfft*nspinor) 207 ! The periodic part of the Bloch wavefunction in real space. 208 209 type(pawcprj_type),allocatable :: Cprj(:,:) 210 ! Cprj(natom,nspinor) 211 ! PAW projected wave function <Proj_i|Cnk> with all NL projectors. 212 213 contains 214 215 procedure :: free => wave_free 216 ! Free memory 217 218 procedure :: copy => wave_copy 219 ! Copy object. 220 221 end type wave_t 222 223 public :: wave_init
m_wfd/wfd_change_ngfft [ Functions ]
NAME
wfd_change_ngfft
FUNCTION
Reallocate and reinitialize internal tables for performing FFTs of wavefunctions.
INPUTS
Cryst<crystal_t>=Info on unit cell. Psps<pseudopotential_type>=Pseudopotential info. new_ngfft(18)=FFT descriptor for the new FFT mesh.
SIDE EFFECTS
Wfd<wfd_t>=Wavefunction descriptor with new internal tables for FFT defined by new_ngfft.
SOURCE
3679 subroutine wfd_change_ngfft(Wfd, Cryst, Psps, new_ngfft) 3680 3681 !Arguments ------------------------------------ 3682 !scalars 3683 integer,intent(in) :: new_ngfft(18) 3684 type(crystal_t),intent(in) :: Cryst 3685 type(pseudopotential_type),intent(in) :: Psps 3686 class(wfd_t),intent(inout) :: Wfd 3687 3688 !Local variables ------------------------------ 3689 !scalars 3690 integer,parameter :: npw0=0 3691 integer :: npw_k, ik_ibz, istwf_k, is, ik, ib 3692 logical :: iscompatibleFFT 3693 !character(len=500) :: msg 3694 !arrays 3695 integer,allocatable :: kg_k(:,:) 3696 3697 !************************************************************************ 3698 3699 if (all(Wfd%ngfft(1:3) == new_ngfft(1:3)) ) RETURN ! Nothing to do. 3700 3701 if (Wfd%prtvol > 0) then 3702 call wrtout(std_out, & 3703 sjoin(" Changing FFT mesh for wavefunctions: ",ltoa(Wfd%ngfft(1:3)), " ==> ", ltoa(new_ngfft(1:3)))) 3704 end if 3705 3706 ! Change FFT dimensions. 3707 Wfd%ngfft = new_ngfft 3708 Wfd%mgfft = MAXVAL(new_ngfft(1:3)) 3709 Wfd%nfftot = PRODUCT(new_ngfft(1:3)) 3710 Wfd%nfft = Wfd%nfftot ! No FFT parallelism. 3711 3712 ! Re-initialize fft distribution 3713 call destroy_distribfft(Wfd%MPI_enreg%distribfft) 3714 call init_distribfft(Wfd%MPI_enreg%distribfft,'c',Wfd%MPI_enreg%nproc_fft,new_ngfft(2),new_ngfft(3)) 3715 3716 ABI_REMALLOC(Wfd%ph1d,(2,3*(2*Wfd%mgfft+1)*Cryst%natom)) 3717 call getph(Cryst%atindx,Cryst%natom,Wfd%ngfft(1),Wfd%ngfft(2),Wfd%ngfft(3),Wfd%ph1d,Cryst%xred) 3718 3719 ! Recalculate FFT tables. 3720 ! Calculate the FFT index of $ R^{-1} (r-\tau) $ used to symmetrize u_Rk. 3721 #ifdef FC_LLVM 3722 ! LLVM 16 doesn't recognize this macro here 3723 ABI_REMALLOC(Wfd%irottb, (Wfd%nfftot,Cryst%nsym) ) 3724 #else 3725 ABI_REMALLOC(Wfd%irottb, (Wfd%nfftot,Cryst%nsym)) 3726 #endif 3727 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,Wfd%ngfft,Wfd%irottb,iscompatibleFFT) 3728 3729 if (.not. iscompatibleFFT) then 3730 ABI_WARNING("FFT mesh not compatible with symmetries. Wavefunction symmetrization should not be done in r-space!") 3731 end if 3732 3733 ! Is the new real space FFT mesh compatible with the rotational part? 3734 Wfd%rfft_is_symok = check_rot_fft(Cryst%nsym,Cryst%symrel,Wfd%ngfft(1),Wfd%ngfft(2),Wfd%ngfft(3)) 3735 3736 ! Reallocate ur buffers with correct dimensions. 3737 do is=1,size(wfd%s) 3738 do ik=1,size(wfd%s(is)%k) 3739 do ib=1,size(wfd%s(is)%k(ik)%b) 3740 call wfd%s(is)%k(ik)%b(ib)%free(what="R") 3741 end do 3742 end do 3743 end do 3744 3745 ! Reinit Kdata_t 3746 do ik_ibz=1,Wfd%nkibz 3747 if (any(wfd%bks2wfd(1, :, ik_ibz, :) /= 0)) then 3748 istwf_k = Wfd%istwfk(ik_ibz) 3749 npw_k = Wfd%Kdata(ik_ibz)%npw 3750 ABI_MALLOC(kg_k, (3,npw_k)) 3751 kg_k = Wfd%Kdata(ik_ibz)%kg_k 3752 call Wfd%Kdata(ik_ibz)%free() 3753 call Wfd%Kdata(ik_ibz)%init(Cryst,Psps,Wfd%kibz(:,ik_ibz),istwf_k,new_ngfft,Wfd%MPI_enreg,kg_k=kg_k) 3754 ABI_FREE(kg_k) 3755 end if 3756 end do 3757 3758 end subroutine wfd_change_ngfft
m_wfd/wfd_copy_cg [ Functions ]
NAME
wfd_copy_cg
FUNCTION
Return a copy u(g) in a real(dp) array. Useful if we have to interface the wavefunction descriptor with Abinit code expecting cg(2,npw_k*nspinor) arrays The routine takes also into account the fact that the ug in wfs could be stored in single-precision.
INPUTS
wfd<wfd_t>=the wavefunction descriptor. band=Band index. ik_ibz=Index of the k-point in the IBZ. spin=Spin index
OUTPUT
cg(npw_k*nspinor)=The wavefunction in real space in the Abinit cg convention.
SOURCE
1551 subroutine wfd_copy_cg(wfd, band, ik_ibz, spin, cg) 1552 1553 !Arguments ------------------------------------ 1554 !scalars 1555 integer,intent(in) :: band,ik_ibz,spin 1556 class(wfd_t),intent(in) :: wfd 1557 !arrays 1558 real(dp),intent(out) :: cg(2,*) ! npw_k*wfd%nspinor) 1559 1560 !Local variables ------------------------------ 1561 !scalars 1562 integer :: siz 1563 type(wave_t),pointer :: wave 1564 character(len=500) :: msg 1565 !************************************************************************ 1566 1567 ABI_CHECK(wfd%get_wave_ptr(band, ik_ibz, spin, wave, msg) == 0, msg) 1568 1569 if (.not. wave%has_ug == WFD_STORED) then 1570 write(msg,'(a,3(i0,1x),a)')" ug for (band, ik_ibz, spin): ",band, ik_ibz, spin," is not stored in memory!" 1571 ABI_ERROR(msg) 1572 end if 1573 1574 siz = wfd%npwarr(ik_ibz) * wfd%nspinor 1575 #ifdef HAVE_GW_DPC 1576 call zcopy(siz, wave%ug, 1, cg, 1) 1577 #else 1578 cg(1,1:siz) = dble(wave%ug) 1579 cg(2,1:siz) = aimag(wave%ug) 1580 #endif 1581 1582 end subroutine wfd_copy_cg
m_wfd/wfd_dump_errinfo [ Functions ]
NAME
wfd_dump_errinfo
FUNCTION
INPUTS
Wfd<wfd_t>=
OUTPUT
SOURCE
3318 subroutine wfd_dump_errinfo(Wfd,onfile) 3319 3320 !Arguments ------------------------------------ 3321 !scalars 3322 logical,optional,intent(in) :: onfile 3323 class(wfd_t),intent(in) :: Wfd 3324 3325 !Local variables ------------------------------ 3326 !scalars 3327 integer :: ik_ibz,spin,band,how_manyb,unt_dbg 3328 character(len=10) :: strank 3329 character(len=500) :: msg 3330 character(len=fnlen) :: fname_dbg 3331 !arrays 3332 integer :: my_band_list(Wfd%mband) 3333 3334 !************************************************************************ 3335 3336 unt_dbg=std_out 3337 3338 if (present(onfile)) then 3339 if (onfile) then 3340 call int2char10(Wfd%my_rank,strank) 3341 fname_dbg = "WFD_DEBUG_RANK"//TRIM(strank) 3342 if (open_file(fname_dbg,msg,newunit=unt_dbg,form="formatted") /= 0) then 3343 ABI_ERROR(msg) 3344 end if 3345 end if 3346 end if 3347 3348 write(unt_dbg,*)" (k,b,s) states owned by rank: ",Wfd%my_rank 3349 do spin=1,Wfd%nsppol 3350 do ik_ibz=1,Wfd%nkibz 3351 write(unt_dbg,*)" ug stored at (ik_ibz, spin) ",ik_ibz,spin 3352 call wfd%mybands(ik_ibz, spin, how_manyb, my_band_list, how="Stored") 3353 write(unt_dbg,*) (my_band_list(band),band=1,how_manyb) 3354 end do 3355 end do 3356 3357 end subroutine wfd_dump_errinfo
m_wfd/wfd_extract_cgblock [ Functions ]
NAME
wfd_extract_cgblock
FUNCTION
This routine extract a block of wavefunctions for a given spin and k-points. The wavefunctions are stored in a real(dp) array with the same convention as the one used in the GS part of Abinit, i.e cg_block(2,nspinor*npw_k*num_bands)
INPUTS
Wfd<wfd_t>=Wavefunction descriptor. band_list(:)=List of bands to extract ik_ibz=k-point index spin=Spin index.
OUTPUT
cgblock(nspinor*npw_k*num_bands)=A contiguous block of memory with the set of u(g)
SOURCE
2272 subroutine wfd_extract_cgblock(Wfd,band_list,ik_ibz,spin,cgblock) 2273 2274 !Arguments ------------------------------------ 2275 !scalars 2276 integer,intent(in) :: ik_ibz,spin 2277 class(wfd_t),intent(in) :: Wfd 2278 !arrays 2279 integer,intent(in) :: band_list(:) 2280 real(dp),intent(out) :: cgblock(:,:) 2281 2282 !Local variables ------------------------------ 2283 !scalars 2284 integer :: ii,band,start,istop,npw_k 2285 character(len=500) :: msg 2286 type(wave_t),pointer :: wave 2287 2288 !************************************************************************ 2289 2290 npw_k = Wfd%npwarr(ik_ibz) 2291 2292 if (size(cgblock, dim=1) /= 2) then 2293 ABI_ERROR("Wrong size(1) in assumed shape array") 2294 end if 2295 2296 if (size(cgblock, dim=2) /= Wfd%nspinor* npw_k * size(band_list)) then 2297 ABI_ERROR("Wrong size in assumed shape array") 2298 end if 2299 2300 start = 1 2301 do ii=1,size(band_list) 2302 band = band_list(ii) 2303 ABI_CHECK(wfd%get_wave_ptr(band, ik_ibz, spin, wave, msg) == 0, msg) 2304 if (wave%has_ug /= WFD_STORED) then 2305 write(msg,"(3(a,i0),a)")"u(g) for band: ",band,", ik_ibz: ",ik_ibz,", spin: ",spin," is not stored!" 2306 ABI_ERROR(msg) 2307 end if 2308 istop = start + Wfd%nspinor*npw_k - 1 2309 cgblock(1,start:istop) = REAL(wave%ug) 2310 cgblock(2,start:istop) = AIMAG(wave%ug) 2311 start = start + Wfd%nspinor * npw_k 2312 end do 2313 2314 end subroutine wfd_extract_cgblock
m_wfd/wfd_free [ Functions ]
NAME
wfd_free
FUNCTION
Free the memory allocated in the wfd_t data type.
SOURCE
1106 subroutine wfd_free(Wfd) 1107 1108 !Arguments ------------------------------------ 1109 !scalars 1110 class(wfd_t),intent(inout) :: Wfd 1111 1112 !Local variables ------------------------------ 1113 !scalars 1114 integer :: ib, ik, is 1115 1116 !************************************************************************ 1117 1118 ! integer. 1119 ABI_SFREE(Wfd%irottb) 1120 ABI_SFREE(Wfd%istwfk) 1121 ABI_SFREE(Wfd%nband) 1122 ABI_SFREE(Wfd%indlmn) 1123 ABI_SFREE(Wfd%nlmn_atm) 1124 ABI_SFREE(Wfd%nlmn_sort) 1125 ABI_SFREE(Wfd%nlmn_type) 1126 ABI_SFREE(Wfd%npwarr) 1127 1128 select type (wfd) 1129 class is (wfdgw_t) 1130 ABI_SFREE(Wfd%bks_tab) 1131 end select 1132 1133 if (allocated(wfd%s)) then 1134 do is=1,size(wfd%s) 1135 do ik=1,size(wfd%s(is)%k) 1136 do ib=1,size(wfd%s(is)%k(ik)%b) 1137 call wfd%s(is)%k(ik)%b(ib)%free() 1138 end do 1139 ABI_FREE(wfd%s(is)%k(ik)%b) 1140 end do 1141 ABI_FREE(wfd%s(is)%k) 1142 end do 1143 ABI_FREE(wfd%s) 1144 end if 1145 1146 ABI_SFREE(wfd%my_nkspin) 1147 ABI_SFREE(wfd%bks2wfd) 1148 1149 ! real arrays. 1150 ABI_SFREE(Wfd%kibz) 1151 ABI_SFREE(Wfd%ph1d) 1152 1153 ! logical arrays. 1154 ABI_SFREE(Wfd%keep_ur) 1155 1156 ! datatypes. 1157 if (allocated(Wfd%Kdata)) then 1158 call kdata_free(Wfd%Kdata) 1159 ABI_FREE(Wfd%Kdata) 1160 end if 1161 1162 call destroy_mpi_enreg(Wfd%MPI_enreg) 1163 1164 end subroutine wfd_free
m_wfd/wfd_get_cprj [ Functions ]
NAME
wfd_get_cprj
FUNCTION
Return a copy of Cprj either by calculating it on-the-fly or by just retrieving the data already stored in the data type.
INPUTS
Wfd<wfd_t>=the wavefunction descriptor. band=Band index. ik_ibz=Index of the k-point in the IBZ. spin=Spin index sorted=.TRUE. if the output cprj matrix elements have to be sorted by atom type.
OUTPUT
Cprj_out(Wfd%natom,Wfd%nspinor) <type(pawcprj_type)>=Unsorted matrix elements.
SOURCE
3571 subroutine wfd_get_cprj(Wfd, band, ik_ibz, spin, Cryst, Cprj_out, sorted) 3572 3573 !Arguments ------------------------------------ 3574 !scalars 3575 integer,intent(in) :: band,ik_ibz,spin 3576 logical,intent(in) :: sorted 3577 class(wfd_t),intent(inout) :: Wfd 3578 type(crystal_t),intent(in) :: Cryst 3579 !arrays 3580 type(pawcprj_type),intent(inout) :: Cprj_out(Wfd%natom,Wfd%nspinor) 3581 3582 !Local variables ------------------------------ 3583 !scalars 3584 integer,parameter :: choice1=1,idir0=0 3585 integer :: want_order,iatom,sidx 3586 character(len=500) :: msg 3587 type(wave_t),pointer :: wave 3588 3589 !************************************************************************ 3590 3591 want_order=CPR_RANDOM; if (sorted) want_order=CPR_SORTED 3592 3593 ABI_CHECK(wfd%get_wave_ptr(band, ik_ibz, spin, wave, msg) == 0, msg) 3594 3595 select case (wave%has_cprj) 3596 3597 case (WFD_NOWAVE, WFD_ALLOCATED) 3598 ! Have to calculate it! 3599 if (.not. wave%has_ug == WFD_STORED) then 3600 write(msg,'(a,3(i0,1x),a)')" ug for (band, ik_ibz, spin): ",band,ik_ibz,spin," is not stored in memory!" 3601 ABI_ERROR(msg) 3602 end if 3603 ! Get cprj. 3604 call wfd%ug2cprj(band,ik_ibz,spin,choice1,idir0,Wfd%natom,Cryst,Cprj_out,sorted=sorted) 3605 3606 if (wave%has_cprj == WFD_ALLOCATED) then 3607 ! Store it. 3608 if (want_order == wave%cprj_order) then 3609 call pawcprj_copy(Cprj_out, wave%Cprj) 3610 wave%has_cprj = WFD_STORED 3611 3612 else 3613 ! Have to reorder cprj_out 3614 select case (want_order) 3615 case (CPR_SORTED) 3616 do iatom=1,Cryst%natom 3617 sidx = Cryst%atindx(iatom) ! random --> sorted table. 3618 call pawcprj_copy(Cprj_out(sidx:sidx,:), wave%Cprj(iatom:iatom,:)) 3619 end do 3620 case (CPR_RANDOM) 3621 do sidx=1,Cryst%natom 3622 iatom = Cryst%atindx1(sidx) ! sorted --> random table. 3623 call pawcprj_copy(Cprj_out(iatom:iatom,:), wave%Cprj(sidx:sidx,:)) 3624 end do 3625 case default 3626 ABI_ERROR(sjoin("Wrong value for want_order:", itoa(want_order))) 3627 end select 3628 end if 3629 end if 3630 3631 case (WFD_STORED) 3632 ! copy it back. 3633 if (want_order == wave%cprj_order) then 3634 call pawcprj_copy(wave%Cprj,Cprj_out) 3635 3636 else 3637 select case (want_order) 3638 case (CPR_SORTED) 3639 do iatom=1,Cryst%natom 3640 sidx = Cryst%atindx(iatom) ! random --> sorted table. 3641 call pawcprj_copy(wave%Cprj(iatom:iatom,:),Cprj_out(sidx:sidx,:)) 3642 end do 3643 case (CPR_RANDOM) 3644 do sidx=1,Cryst%natom 3645 iatom = Cryst%atindx1(sidx) ! sorted --> random table. 3646 call pawcprj_copy(wave%Cprj(sidx:sidx,:),Cprj_out(iatom:iatom,:)) 3647 end do 3648 case default 3649 ABI_ERROR(sjoin("Wrong value for want_order:", itoa(want_order))) 3650 end select 3651 end if 3652 3653 case default 3654 ABI_BUG(sjoin("Wrong has_cprj: ", itoa(wave%has_cprj))) 3655 end select 3656 3657 end subroutine wfd_get_cprj
m_wfd/wfd_get_many_ur [ Functions ]
NAME
wfd_get_many_ur
FUNCTION
Get many wave functions in real space, either by doing a G-->R FFT or by just retrieving the data already stored in Wfd.
INPUTS
Wfd<wfd_t>=the wavefunction descriptor. ndat=Number of wavefunctions required bands(:)=Band indices. ik_ibz=Index of the k-point in the IBZ. spin=Spin index
OUTPUT
ur(Wfd%nfft*Wfd%nspinor*SIZE(bands))=The wavefunction in real space.
SOURCE
1505 subroutine wfd_get_many_ur(Wfd, bands, ik_ibz, spin, ur) 1506 1507 !Arguments ------------------------------------ 1508 !scalars 1509 integer,intent(in) :: ik_ibz,spin 1510 class(wfd_t),intent(inout) :: Wfd 1511 !arrays 1512 integer,intent(in) :: bands(:) 1513 complex(gwpc),intent(out) :: ur(Wfd%nfft*Wfd%nspinor*SIZE(bands)) 1514 1515 !Local variables ------------------------------ 1516 !scalars 1517 integer :: dat,ptr,band 1518 !************************************************************************ 1519 1520 do dat=1,SIZE(bands) 1521 band = bands(dat) 1522 ptr = 1 + (dat-1)*Wfd%nfft*Wfd%nspinor 1523 call wfd%get_ur(band,ik_ibz,spin,ur(ptr)) 1524 end do 1525 1526 end subroutine wfd_get_many_ur
m_wfd/wfd_get_socpert [ Functions ]
NAME
wfd_get_socpert
FUNCTION
INPUTS
cryst<crystal_t>= data type gathering info on symmetries and unit cell psps<pseudopotential_type>=variables related to pseudopotentials pawtab(psps%ntypat) <type(pawtab_type)>=paw tabulated starting data paw_ij(natom)<type(paw_ij_type)>=data structure containing PAW arrays given on (i,j) channels.
OUTPUT
SOURCE
5258 !!! subroutine wfd_get_socpert(wfd, cryst, psps, pawtab, bks_mask, soc_bks) 5259 !!! 5260 !!! !use m_pawcprj 5261 !!! use m_hamiltonian, only : destroy_hamiltonian, init_hamiltonian, & 5262 !!! load_spin_hamiltonian,load_k_hamiltonian, gs_hamiltonian_type 5263 !!! 5264 !!! implicit none 5265 !!! 5266 !!! !Arguments ------------------------------------ 5267 !!! !scalars 5268 !!! type(wfd_t),target,intent(inout) :: wfd 5269 !!! type(crystal_t),intent(in) :: cryst 5270 !!! type(pseudopotential_type),intent(in) :: psps 5271 !!! ! arrays 5272 !!! logical,intent(in) :: bks_mask(wfd%mband, wfd%nkibz, wfd%nsppol) 5273 !!! real(dp),allocatable,intent(out) :: osoc_bks(:, :, :) 5274 !!! type(Pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 5275 !!! !type(paw_ij_type),intent(in) :: paw_ij(cryst%natom*psps%usepaw) 5276 !!! 5277 !!! !Local variables ------------------------------ 5278 !!! !scalars 5279 !!! integer,parameter :: nspinor2=2,nspden4=4,nsppol1=1,spin1=1 5280 !!! integer,parameter :: ndat1=1,nnlout0=0,tim_nonlop0=0,idir0=0 !,ider0=0, 5281 !!! integer :: natom,band,spin,ik_ibz,npw_k,istwf_k,nkpg !,ig,optder,matblk,mkmem_,nkpg,dimffnl,nspinortot 5282 !!! integer :: choice,cpopt,cp_dim,paw_opt,signs,ierr 5283 !!! !character(len=500) :: msg 5284 !!! type(gs_hamiltonian_type) :: ham_k 5285 !!! !arrays 5286 !!! integer :: bks_distrb(wfd%mband, wfd%nkibz, wfd%nsppol) 5287 !!! integer, ABI_CONTIGUOUS pointer :: kg_k(:,:) 5288 !!! !real(dp) :: kptns_(3,1),ylmgr_dum(1,1,1),shifts(3) 5289 !!! !real(dp),allocatable :: ylm_k(:,:),dum_ylm_gr_k(:,:,:) 5290 !!! !real(dp),pointer :: ffnl_k(:,:,:,:) 5291 !!! real(dp) :: kpoint(3),dum_enlout(0),dummy_lambda(1),soc(2) 5292 !!! real(dp),allocatable :: kpg_k(:,:),vnl_psi(:,:),vectin(:,:) !,s_psi(:,:) 5293 !!! real(dp),allocatable :: opaw_psi(:,:) !2, npw_k*wfd%nspinor*wfd%usepaw) ! <G|1+S|Cnk> 5294 !!! real(dp),ABI_CONTIGUOUS pointer :: ffnl_k(:,:,:,:),ph3d_k(:,:,:) 5295 !!! type(pawcprj_type),allocatable :: cprj(:,:) 5296 !!! 5297 !!! !************************************************************************ 5298 !!! 5299 !!! ABI_CHECK(wfd%paral_kgb == 0, "paral_kgb not coded") 5300 !!! 5301 !!! natom = cryst%natom 5302 !!! 5303 !!! signs = 2 ! => apply the non-local operator to a function in G-space. 5304 !!! choice = 1 ! => <G|V_nonlocal|vectin>. 5305 !!! cpopt =-1; paw_opt= 0 5306 !!! if (wfd%usepaw==1) then 5307 !!! paw_opt=4 ! both PAW nonlocal part of H (Dij) and overlap matrix (Sij) 5308 !!! cpopt=3 ! <p_lmn|in> are already in memory 5309 !!! 5310 !!! cp_dim = ((cpopt+5) / 5) 5311 !!! ABI_MALLOC(cprj, (natom, nspinor2*cp_dim)) 5312 !!! call pawcprj_alloc(cprj, 0, wfd%nlmn_sort) 5313 !!! end if 5314 !!! 5315 !!! ! Initialize the Hamiltonian on the coarse FFT mesh. 5316 !!! call init_hamiltonian(ham_k, psps, pawtab, nspinor2, nsppol1, nspden4, natom, cryst%typat, cryst%xred, & 5317 !!! wfd%nfft, wfd%mgfft, wfd%ngfft, cryst%rprimd, wfd%nloalg) 5318 !!! !ham_k%ekb(:,:,1) = zero 5319 !!! 5320 !!! ! Continue to prepare the GS Hamiltonian (note spin1) 5321 !!! call load_spin_hamiltonian(ham_k, spin1, with_nonlocal=.True.) 5322 !!! 5323 !!! ! Distribute (b, k, s) states. 5324 !!! call wfd%bks_distrb(bks_distrb, bks_mask=bks_mask) 5325 !!! 5326 !!! ABI_CALLOC(osoc_bks, (wfd%mband, wfd%nkibz, wfd%nsppol)) 5327 !!! osoc_bks = zero 5328 !!! 5329 !!! do spin=1,wfd%nsppol 5330 !!! do ik_ibz=1,wfd%nkibz 5331 !!! if (all(bks_distrb(:, ik_ibz, spin) /= wfd%my_rank)) cycle 5332 !!! 5333 !!! kpoint = wfd%kibz(:, ik_ibz) 5334 !!! npw_k = wfd%Kdata(ik_ibz)%npw; istwf_k = wfd%istwfk(ik_ibz) 5335 !!! ABI_CHECK(istwf_k == 1, "istwf_k must be 1 if SOC term is computed with perturbation theory.") 5336 !!! kg_k => wfd%kdata(ik_ibz)%kg_k 5337 !!! ffnl_k => wfd%Kdata(ik_ibz)%fnl_dir0der0 5338 !!! ph3d_k => wfd%Kdata(ik_ibz)%ph3d 5339 !!! 5340 !!! ABI_MALLOC(vectin, (2, npw_k * nspinor2)) 5341 !!! ABI_MALLOC(vnl_psi, (2, npw_k * nspinor2)) 5342 !!! !ABI_MALLOC(cvnl_psi, (npw_k * nspinor2)) 5343 !!! !ABI_MALLOC(s_psi, (2, npw_k * nspinor2 * psps%usepaw)) 5344 !!! 5345 !!! ! Compute (k+G) vectors (only if psps%useylm=1) 5346 !!! nkpg = 3 * wfd%nloalg(3) 5347 !!! ABI_MALLOC(kpg_k, (npw_k, nkpg)) 5348 !!! if (nkpg > 0) then 5349 !!! call mkkpg(kg_k, kpg_k, kpoint, nkpg, npw_k) 5350 !!! end if 5351 !!! 5352 !!! ! Load k-dependent part in the Hamiltonian datastructure 5353 !!! !matblk = min(NLO_MINCAT, maxval(ham_k%nattyp)); if (wfd%nloalg(2) > 0) matblk = natom 5354 !!! !ABI_MALLOC(ph3d_k,(2, npw_k, matblk)) 5355 !!! call load_k_hamiltonian(ham_k, kpt_k=kpoint, npw_k=npw_k, istwf_k=istwf_k, kg_k=kg_k, & 5356 !!! kpg_k=kpg_k, ffnl_k=ffnl_k, ph3d_k=ph3d_k, compute_ph3d=(wfd%paral_kgb/=1)) 5357 !!! 5358 !!! ! THIS PART IS NEEDED FOR THE CALL TO opernl although some quantities won't be used. 5359 !!! ! Now I do things cleanly then we try to pass zero-sized arrays! 5360 !!! !ABI_MALLOC(ylm_k, (npw_k, psps%mpsang**2 * psps%useylm)) 5361 !!! !if (psps%useylm == 1) then 5362 !!! ! kptns_(:,1) = k4intp; optder = 0; mkmem_ = 1 5363 !!! ! ABI_MALLOC(dum_ylm_gr_k,(npw_k,3+6*(optder/2),psps%mpsang**2)) 5364 !!! ! ! Here mband is not used if paral_compil_kpt=0 5365 !!! ! call initylmg(cryst%gprimd, kg_k, kptns_, mkmem_, wfd%MPI_enreg, psps%mpsang, npw_k, [1], 1,& 5366 !!! ! [npw_k], 1, optder, cryst%rprimd, ylm_k, dum_ylm_gr_k) 5367 !!! ! ABI_FREE(dum_ylm_gr_k) 5368 !!! !end if 5369 !!! 5370 !!! ! ======================================================== 5371 !!! ! ==== Compute nonlocal form factors ffnl at all (k+G) ==== 5372 !!! ! ======================================================== 5373 !!! !dimffnl = 1 + ider0 ! Derivatives are not needed. 5374 !!! !ABI_MALLOC(ffnl_k, (npw_k, dimffnl, psps%lmnmax, psps%ntypat)) 5375 !!! !! ffnl_k => Kdata%fnl_dir0der0 5376 !!! !call mkffnl(psps%dimekb, dimffnl, psps%ekb, ffnl_k, psps%ffspl, cryst%gmet, cryst%gprimd, ider0, idir0, psps%indlmn,& 5377 !!! ! kg_k, kpg_k, k4intp, psps%lmnmax, psps%lnmax, psps%mpsang, psps%mqgrid_ff, nkpg, npw_k, & 5378 !!! ! psps%ntypat, psps%pspso, psps%qgrid_ff, cryst%rmet, psps%usepaw, psps%useylm, ylm_k, ylmgr_dum) 5379 !!! !ABI_FREE(ylm_k) 5380 !!! 5381 !!! ! Calculate <G|Vnl|psi> for this k-point 5382 !!! do band=1,wfd%nband(ik_ibz, spin) 5383 !!! if (bks_distrb(band, ik_ibz, spin) /= wfd%my_rank) cycle 5384 !!! 5385 !!! ! Input wavefunction coefficients <G|Cnk>. 5386 !!! ! vectin, (2, npw_k * nspinor2)) 5387 !!! if (spin == 1) then 5388 !!! vectin(1, 1:npw_k) = dble(wfd%wave(band, ik_ibz, spin)%ug) 5389 !!! vectin(2, 1:npw_k) = aimag(wfd%wave(band, ik_ibz, spin)%ug) 5390 !!! vectin(:, npw_k+1:) = zero 5391 !!! else 5392 !!! vectin(:, 1:npw_k) = zero 5393 !!! vectin(1, npw_k+1:) = dble(wfd%wave(band, ik_ibz, spin)%ug) 5394 !!! vectin(2, npw_k+1:) = aimag(wfd%wave(band, ik_ibz, spin)%ug) 5395 !!! end if 5396 !!! 5397 !!! if (wfd%usepaw == 1) call wfd%get_cprj(band, ik_ibz, spin, cryst, cprj, sorted=.True.) 5398 !!! 5399 !!! ! TODO: consistency check for only_SO 5400 !!! call nonlop(choice, cpopt, cprj, dum_enlout, ham_k, idir0, dummy_lambda, wfd%mpi_enreg, ndat1, nnlout0, & 5401 !!! paw_opt, signs, opaw_psi, tim_nonlop0, vectin, vnl_psi, only_SO=1) 5402 !!! 5403 !!! soc = cg_zdotc(npw_k * nspinor2, vectin, vnl_psi) 5404 !!! write(std_out,*)soc * Ha_eV, "for (b, k, s)",band, ik_ibz, spin 5405 !!! osoc_bks(band, ik_ibz, spin) = soc(1) 5406 !!! end do ! band 5407 !!! 5408 !!! !ABI_FREE(ffnl_k) 5409 !!! !ABI_FREE(ph3d_k) 5410 !!! ABI_FREE(vectin) 5411 !!! ABI_FREE(vnl_psi) 5412 !!! ABI_FREE(kpg_k) 5413 !!! !ABI_FREE(cvnl_psi) 5414 !!! !ABI_FREE(s_psi) 5415 !!! end do ! ik_ibz 5416 !!! end do ! spin 5417 !!! 5418 !!! call xmpi_sum(osoc_bks, wfd%comm, ierr) 5419 !!! 5420 !!! call destroy_hamiltonian(ham_k) 5421 !!! 5422 !!! if (wfd%usepaw == 1) then 5423 !!! call pawcprj_free(cprj) 5424 !!! ABI_FREE(cprj) 5425 !!! end if 5426 !!! 5427 !!! end subroutine wfd_get_socpert
m_wfd/wfd_get_ug [ Functions ]
NAME
wfd_get_ug
FUNCTION
Get a **copy** of a wave function in G-space.
INPUTS
Wfd<wfd_t>=the data type band=the index of the band. ik_ibz=Index of the k-point in the IBZ spin=spin index
OUTPUT
ug(npw_k*Wfd%nspinor)=The required wavefunction in G-space
SOURCE
2618 subroutine wfd_get_ug(Wfd, band, ik_ibz, spin, ug) 2619 2620 !Arguments ------------------------------------ 2621 !scalars 2622 integer,intent(in) :: band,ik_ibz,spin 2623 class(wfd_t),intent(inout) :: Wfd 2624 !arrays 2625 complex(gwpc),intent(out) :: ug(Wfd%npwarr(ik_ibz)*Wfd%nspinor) 2626 2627 !Local variables ------------------------------ 2628 !scalars 2629 integer :: npw_k 2630 character(len=500) :: msg 2631 type(wave_t),pointer :: wave 2632 !************************************************************************ 2633 2634 ABI_CHECK(wfd%get_wave_ptr(band, ik_ibz, spin, wave, msg) == 0, msg) 2635 2636 if (.not. wave%has_ug == WFD_STORED) then 2637 write(msg,'(a,i0,a,3i0)')" Node ",Wfd%my_rank," doesn't have (band,ik_ibz,spin): ",band,ik_ibz,spin 2638 ABI_BUG(msg) 2639 end if 2640 2641 npw_k = Wfd%npwarr(ik_ibz) 2642 call xcopy(npw_k*Wfd%nspinor, wave%ug, 1, ug, 1) 2643 2644 end subroutine wfd_get_ug
m_wfd/wfd_get_ur [ Functions ]
NAME
wfd_get_ur
FUNCTION
Get a wave function in real space, either by doing a G-->R FFT or by just retrieving the data already stored in Wfd.
INPUTS
Wfd<wfd_t>=the wavefunction descriptor. band=Band index. ik_ibz=Index of the k-point in the IBZ. spin=Spin index
OUTPUT
ur(Wfd%nfft*Wfd%nspinor)=The wavefunction in real space.
SOURCE
1606 subroutine wfd_get_ur(Wfd, band, ik_ibz, spin, ur) 1607 1608 !Arguments ------------------------------------ 1609 !scalars 1610 integer,intent(in) :: band,ik_ibz,spin 1611 class(wfd_t),target,intent(inout) :: Wfd 1612 !arrays 1613 complex(gwpc),intent(out) :: ur(Wfd%nfft*Wfd%nspinor) 1614 1615 !Local variables ------------------------------ 1616 !scalars 1617 integer,parameter :: npw0=0,ndat1=1 1618 integer :: npw_k, nfft, nspinor 1619 character(len=500) :: msg 1620 type(wave_t),pointer :: wave 1621 !arrays 1622 integer,ABI_CONTIGUOUS pointer :: kg_k(:,:),gbound(:,:) 1623 complex(gwpc),ABI_CONTIGUOUS pointer :: ug(:) 1624 !************************************************************************ 1625 1626 npw_k = Wfd%npwarr(ik_ibz) 1627 nfft = Wfd%nfft 1628 nspinor= Wfd%nspinor 1629 1630 ABI_CHECK(wfd%get_wave_ptr(band, ik_ibz, spin, wave, msg) == 0, msg) 1631 1632 select case (wave%has_ur) 1633 1634 case (WFD_NOWAVE, WFD_ALLOCATED) 1635 ! FFT is required. 1636 if (.not. wave%has_ug == WFD_STORED) then 1637 write(msg,'(a,3(i0,1x),a)')" ug for (band, ik_ibz, spin): ",band,ik_ibz,spin," is not stored in memory!" 1638 ABI_ERROR(msg) 1639 end if 1640 1641 ug => wave%ug 1642 kg_k => Wfd%Kdata(ik_ibz)%kg_k 1643 gbound => Wfd%Kdata(ik_ibz)%gbound(:,:) 1644 1645 call fft_ug(npw_k,nfft,nspinor,ndat1,Wfd%mgfft,Wfd%ngfft,Wfd%istwfk(ik_ibz),kg_k,gbound,ug,ur) 1646 1647 if (Wfd%keep_ur(band,ik_ibz,spin)) then 1648 ! Store results 1649 if (wave%has_ur == WFD_NOWAVE) then 1650 ! Alloc buffer for ur. 1651 call wave_init(wave, Wfd%usepaw,npw0,nfft,nspinor,Wfd%natom,Wfd%nlmn_atm,CPR_RANDOM) 1652 end if 1653 call xcopy(nfft*nspinor, ur, 1, wave%ur, 1) 1654 wave%has_ur = WFD_STORED 1655 end if 1656 1657 case (WFD_STORED) 1658 ! copy it back. 1659 call xcopy(nfft*nspinor, wave%ur, 1, ur, 1) 1660 1661 case default 1662 ABI_BUG(sjoin("Wrong has_ur:", itoa(wave%has_ur))) 1663 end select 1664 1665 end subroutine wfd_get_ur
m_wfd/wfd_get_wave_prt [ Functions ]
NAME
wfd_get_wave_prt
FUNCTION
Return pointer to the wave object corresponding to the given (band, ik_ibz, spin) indices. If the state is not treated ...
INPUTS
band=Band index. ik_ibz=k-point index spin=Spin index.
SOURCE
2116 integer function wfd_get_wave_ptr(wfd, band, ik_ibz, spin, wave_ptr, msg) result(ierr) 2117 2118 !Arguments ------------------------------------ 2119 !scalars 2120 integer,intent(in) :: ik_ibz, spin, band 2121 class(wfd_t),target,intent(in) :: wfd 2122 type(wave_t),pointer :: wave_ptr 2123 character(len=*),intent(out) :: msg 2124 2125 !Local variables ------------------------------ 2126 !scalars 2127 integer :: ib, ik, is 2128 2129 !************************************************************************ 2130 2131 ierr = 1 2132 if (any(wfd%bks2wfd(:, band, ik_ibz, spin) == 0)) then 2133 write(msg,'(a,i0,a,3(i0,1x))')" MPI rank ",Wfd%my_rank," doesn't have ug for (band, ik_ibz, spin): ",band,ik_ibz,spin 2134 wave_ptr => null(); return 2135 end if 2136 2137 ib = wfd%bks2wfd(1, band, ik_ibz, spin) 2138 ik = wfd%bks2wfd(2, band, ik_ibz, spin) 2139 is = wfd%bks2wfd(3, band, ik_ibz, spin) 2140 wave_ptr => wfd%s(is)%k(ik)%b(ib) 2141 !if (wave_ptr%has_ug /= WFD_STORED) 2142 2143 ierr = 0 2144 2145 end function wfd_get_wave_ptr
m_wfd/wfd_ihave_ug [ Functions ]
NAME
wfd_ihave_ug
FUNCTION
This function is used to ask the processor whether it has a particular ug and with which status.
INPUTS
band=Band index.
ik_ibz=k-point index
spin=Spin index.
[how]=string defining which status is checked.
Possible mutually exclusive values: "Allocated", "Stored".
Only the first character is checked (no case-sensitive)
By default the function returns .TRUE. if the wave is either WFD_ALLOCATED or WFD_STORED.
NOTES
A zero index can be used to inquire the status of a bunch of states. Thus (band,ik_ibz,spin) = (0,1,1) means: Do you have at least one band for the first k-point and the first spin.
SOURCE
2408 pure function wfd_ihave_ug(Wfd, band, ik_ibz, spin, how) 2409 2410 !Arguments ------------------------------------ 2411 !scalars 2412 integer,intent(in) :: band,ik_ibz,spin 2413 logical :: wfd_ihave_ug 2414 character(len=*),optional,intent(in) :: how 2415 class(wfd_t),intent(in) :: Wfd 2416 2417 !Local variables ------------------------------ 2418 !scalars 2419 integer :: ib, ik, is 2420 integer(c_int8_t) :: check2(2) 2421 2422 !************************************************************************ 2423 2424 check2 = [WFD_ALLOCATED, WFD_STORED] 2425 if (present(how)) then 2426 if (firstchar(how, ["A", "a"])) check2 = [WFD_ALLOCATED, WFD_ALLOCATED] 2427 if (firstchar(how, ["S", "s"])) check2 = [WFD_STORED, WFD_STORED] 2428 end if 2429 ib = wfd%bks2wfd(1, band, ik_ibz, spin) 2430 ik = wfd%bks2wfd(2, band, ik_ibz, spin) 2431 is = wfd%bks2wfd(3, band, ik_ibz, spin) 2432 wfd_ihave_ug = .False. 2433 if (ib /= 0) wfd_ihave_ug = any(wfd%s(is)%k(ik)%b(ib)%has_ug == check2) 2434 2435 end function wfd_ihave_ug
m_wfd/wfd_init [ Functions ]
NAME
wfd_init
FUNCTION
Initialize the object.
INPUTS
Cryst<crystal_t>=Object defining the unit cell and its symmetries. Pawtab(ntypat*usepaw)<type(pawtab_type)>=PAW tabulated starting data. Psps<Pseudopotential_type>=datatype storing data on the pseudopotentials. ngfft(18)=All needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft nkibz=Number of irreducible k-points. nsppol=Number of independent spin polarizations. nspden=Number of density components. nspinor=Number of spinorial components. ecut=Cutoff energy in Hartree ecutsm=Smearing for kinetic energy dilatmx mband nband(nkibz,nsppol) keep_ur(mband,nkibz,nsppol)=Option for memory storage of u(r). istwfk(nkibz)=Storage mode. kibz(3,nkibz)=Reduced coordinates of the k-points. nloalg(3)=Governs the choice of the algorithm for nonlocal operator. See doc. prtvol=Verbosity level. comm=MPI communicator.
OUTPUT
Initialize the object with basic dimensions, allocate also memory for u(g) and u(r) according to keep_ur %ug in G-space are always allocated. %ur in r-space only if keep_ur.
SOURCE
842 subroutine wfd_init(Wfd,Cryst,Pawtab,Psps,keep_ur,mband,nband,nkibz,nsppol,bks_mask,& 843 & nspden,nspinor,ecut,ecutsm,dilatmx,istwfk,kibz,ngfft,nloalg,prtvol,pawprtvol,comm,& 844 & use_fnl_dir0der0) ! optional 845 846 !Arguments ------------------------------------ 847 !scalars 848 integer,intent(in) :: mband,comm,prtvol,pawprtvol 849 integer,intent(in) :: nkibz,nsppol,nspden,nspinor 850 real(dp),intent(in) :: ecut,ecutsm,dilatmx 851 type(crystal_t),intent(in) :: Cryst 852 type(pseudopotential_type),intent(in) :: Psps 853 class(wfd_t),intent(inout) :: Wfd 854 !array 855 integer,intent(in) :: ngfft(18),istwfk(nkibz),nband(nkibz,nsppol) 856 integer,intent(in) :: nloalg(3) 857 real(dp),intent(in) :: kibz(3,nkibz) 858 logical,intent(in) :: bks_mask(mband,nkibz,nsppol) 859 logical,intent(in) :: keep_ur(mband,nkibz,nsppol) 860 logical,intent(in),optional :: use_fnl_dir0der0 861 type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Psps%usepaw) 862 863 !Local variables ------------------------------ 864 !scalars 865 integer,parameter :: nfft0=0,mpw0=0,ikg0=0 866 integer :: ik_ibz,spin,band,mpw,exchn2n3d,istwf_k,npw_k,iatom,itypat,iat 867 integer :: cnt_b, cnt_k, cnt_s, ierr 868 real(dp) :: ug_size,ur_size,cprj_size, bks_size 869 real(dp) :: cpu, wall, gflops 870 logical :: iscompatibleFFT 871 character(len=500) :: msg 872 !arrays 873 integer :: dum_kg(3,0) 874 real(dp) :: kpoint(3) 875 !integer :: my_band_list(Wfd%mband) 876 877 !************************************************************************ 878 879 call cwtime(cpu, wall, gflops, "start") 880 881 ! MPI info 882 Wfd%comm = comm 883 Wfd%my_rank = xmpi_comm_rank(Wfd%comm) 884 Wfd%nproc = xmpi_comm_size(Wfd%comm) 885 Wfd%master = 0 886 887 ! Sequential MPI datatype to be passed to abinit routines. 888 call initmpi_seq(Wfd%MPI_enreg) 889 call init_distribfft(Wfd%MPI_enreg%distribfft,'c',Wfd%MPI_enreg%nproc_fft,ngfft(2),ngfft(3)) 890 891 ! TODO: To simply high-level API. 892 !wfd%cryst => cryst 893 !wfd%psps => psps 894 !wfd%pawtab => pawtab 895 896 ! Basic dimensions 897 Wfd%nkibz = nkibz 898 Wfd%nsppol = nsppol 899 Wfd%nspden = nspden 900 Wfd%nspinor = nspinor 901 Wfd%paral_kgb = 0 902 Wfd%nloalg = nloalg 903 904 Wfd%usepaw = Psps%usepaw 905 Wfd%usewvl = 0 ! wavelets are not supported. 906 Wfd%use_fnl_dir0der0 = .false. 907 if(present(use_fnl_dir0der0)) Wfd%use_fnl_dir0der0 = use_fnl_dir0der0 908 Wfd%natom = Cryst%natom 909 Wfd%ntypat = Cryst%ntypat 910 Wfd%lmnmax = Psps%lmnmax 911 Wfd%prtvol = prtvol 912 Wfd%pawprtvol = pawprtvol 913 914 Wfd%ecutsm = ecutsm 915 Wfd%dilatmx = dilatmx 916 917 ABI_MALLOC(Wfd%indlmn,(6, Wfd%lmnmax, Wfd%ntypat)) 918 Wfd%indlmn = Psps%indlmn 919 920 if (Wfd%usepaw==1) then 921 ABI_MALLOC(Wfd%nlmn_atm, (Cryst%natom)) 922 ABI_MALLOC(Wfd%nlmn_type, (Cryst%ntypat)) 923 do iatom=1,Cryst%natom 924 Wfd%nlmn_atm(iatom) = Pawtab(Cryst%typat(iatom))%lmn_size 925 end do 926 927 do itypat=1,Cryst%ntypat 928 Wfd%nlmn_type(itypat)=Pawtab(itypat)%lmn_size 929 end do 930 931 ABI_MALLOC(Wfd%nlmn_sort,(Cryst%natom)) 932 iat=0 ! nlmn dims sorted by atom type. 933 do itypat=1,Cryst%ntypat 934 Wfd%nlmn_sort(iat+1:iat+Cryst%nattyp(itypat))=Pawtab(itypat)%lmn_size 935 iat=iat+Cryst%nattyp(itypat) 936 end do 937 end if 938 939 ABI_MALLOC(Wfd%keep_ur, (mband, nkibz, nsppol)) 940 Wfd%keep_ur = keep_ur 941 942 ! Setup of the FFT mesh 943 Wfd%ngfft = ngfft 944 Wfd%mgfft = MAXVAL (Wfd%ngfft(1:3)) 945 Wfd%nfftot = PRODUCT(Wfd%ngfft(1:3)) 946 Wfd%nfft = Wfd%nfftot ! At present no FFT parallelism. 947 948 Wfd%ecut = ecut 949 950 ! Precalculate the FFT index of $ R^{-1} (r-\tau) $ used to symmetrize u_Rk. 951 ABI_MALLOC(Wfd%irottb,(Wfd%nfftot, Cryst%nsym)) 952 call rotate_FFT_mesh(Cryst%nsym,Cryst%symrel,Cryst%tnons,Wfd%ngfft,Wfd%irottb,iscompatibleFFT) 953 954 if (.not. iscompatibleFFT) then 955 msg = "FFT mesh is not compatible with symmetries. Wavefunction symmetrization might be affected by large errors!" 956 ABI_WARNING(msg) 957 end if 958 959 ! Is the real space mesh compatible with the rotational part? 960 Wfd%rfft_is_symok = check_rot_fft(Cryst%nsym,Cryst%symrel,Wfd%ngfft(1),Wfd%ngfft(2),Wfd%ngfft(3)) 961 962 ABI_MALLOC(Wfd%kibz, (3, Wfd%nkibz)) 963 Wfd%kibz = kibz 964 ABI_MALLOC(Wfd%istwfk, (Wfd%nkibz)) 965 Wfd%istwfk = istwfk 966 967 ! Get the number of planewaves npw_k 968 ! TODO Here we should use ecut_eff instead of ecut 969 ABI_ICALLOC(Wfd%npwarr, (Wfd%nkibz)) 970 exchn2n3d = 0 971 do ik_ibz=1,Wfd%nkibz 972 if (mod(ik_ibz, wfd%nproc) /= wfd%my_rank) cycle ! MPI parallelism. 973 istwf_k = Wfd%istwfk(ik_ibz) 974 kpoint = Wfd%kibz(:,ik_ibz) 975 call kpgsph(Wfd%ecut,exchn2n3d,Cryst%gmet,ikg0,ik_ibz,istwf_k,dum_kg,kpoint,0,Wfd%MPI_enreg,mpw0,npw_k) 976 Wfd%npwarr(ik_ibz)= npw_k 977 end do 978 call xmpi_sum(wfd%npwarr, wfd%comm, ierr) 979 980 mpw = maxval(Wfd%npwarr) 981 982 ABI_MALLOC(Wfd%nband, (nkibz,nsppol)) 983 Wfd%nband = nband; Wfd%mband = mband 984 ABI_CHECK_IEQ(maxval(Wfd%nband), mband, "Wrong mband") 985 986 ! Allocate u(g) and, if required, also u(r) 987 ug_size = one*nspinor*mpw*COUNT(bks_mask) 988 write(msg,'(a,f8.1,a)')' Memory needed for Fourier components u(G): ',two*gwpc*ug_size*b2Mb, ' [Mb] <<< MEM' 989 call wrtout(std_out, msg) 990 #ifdef HAVE_GW_DPC 991 call wrtout(std_out, ' Storing wavefunctions in double precision as `enable_gw_dpc="no"`') 992 call wrtout(std_out, ' Recompile the code with `enable_gw_dpc="no"` to halve memory requirements for the WFs') 993 #else 994 call wrtout(std_out, ' Storing wavefunctions in single precision as `enable_gw_dpc="no"`') 995 #endif 996 997 if (Wfd%usepaw==1) then 998 cprj_size = one * nspinor*SUM(Wfd%nlmn_atm)*COUNT(bks_mask) 999 write(msg,'(a,f8.1,a)')' Memory needed for PAW projections cprj: ',dp*cprj_size*b2Mb,' [Mb] <<< MEM' 1000 call wrtout(std_out, msg) 1001 end if 1002 1003 ur_size = one*nspinor*Wfd%nfft*COUNT(Wfd%keep_ur) 1004 write(msg,'(a,f8.1,a)')' Memory needed for real-space u(r): ',two*gwpc*ur_size*b2Mb,' [Mb] <<< MEM' 1005 call wrtout(std_out, msg) 1006 1007 ! Count the number of spins treated by this proc. 1008 wfd%my_nspins = 0 1009 do spin=1,wfd%nsppol 1010 if (any(bks_mask(:,:,spin))) wfd%my_nspins = wfd%my_nspins + 1 1011 end do 1012 !write(std_out, *)"my_nspins", wfd%my_nspins 1013 ABI_MALLOC(wfd%s, (wfd%my_nspins)) 1014 1015 ! Count the number of kpts in the IBZ treated by this proc. may be spin-dependent. 1016 ABI_ICALLOC(wfd%my_nkspin, (wfd%nsppol)) 1017 cnt_s = 0 1018 do spin=1,wfd%nsppol 1019 do ik_ibz=1,wfd%nkibz 1020 if (any(bks_mask(:,ik_ibz,spin))) wfd%my_nkspin(spin) = wfd%my_nkspin(spin) + 1 1021 end do 1022 if (wfd%my_nkspin(spin) > 0) then 1023 cnt_s = cnt_s + 1 1024 ABI_MALLOC(wfd%s(cnt_s)%k, (wfd%my_nkspin(spin))) 1025 end if 1026 end do 1027 1028 ! Allocate bands in packed format and use bks2wfd to go from global (b,k,s) index to local index. 1029 ABI_ICALLOC(wfd%bks2wfd, (3, wfd%mband, wfd%nkibz, wfd%nsppol)) 1030 cnt_s = 0 1031 do spin=1,wfd%nsppol 1032 if (wfd%my_nkspin(spin) == 0) cycle 1033 cnt_s = cnt_s + 1 1034 cnt_k = 0 1035 do ik_ibz=1,wfd%nkibz 1036 cnt_b = count(bks_mask(:, ik_ibz, spin)) 1037 if (cnt_b == 0) cycle 1038 cnt_k = cnt_k + 1 1039 ABI_MALLOC(wfd%s(cnt_s)%k(cnt_k)%b, (cnt_b)) 1040 cnt_b = 0 1041 npw_k = Wfd%npwarr(ik_ibz) 1042 do band=1,Wfd%nband(ik_ibz, spin) 1043 if (bks_mask(band, ik_ibz, spin)) then 1044 cnt_b = cnt_b + 1 1045 call wave_init(wfd%s(cnt_s)%k(cnt_k)%b(cnt_b), & 1046 Wfd%usepaw, npw_k, nfft0, Wfd%nspinor, Wfd%natom, Wfd%nlmn_atm, CPR_RANDOM) 1047 wfd%bks2wfd(:, band, ik_ibz, spin) = [cnt_b, cnt_k, cnt_s] 1048 end if 1049 end do 1050 end do 1051 end do 1052 ! 1053 ! =================================================== 1054 ! ==== Precalculate nonlocal form factors for PAW ==== 1055 ! =================================================== 1056 ! 1057 ! Calculate 1-dim structure factor phase information. 1058 ABI_MALLOC(Wfd%ph1d, (2, 3*(2*Wfd%mgfft+1)*Wfd%natom)) 1059 call getph(Cryst%atindx, Wfd%natom, Wfd%ngfft(1), Wfd%ngfft(2), Wfd%ngfft(3), Wfd%ph1d, Cryst%xred) 1060 1061 ! TODO: This one will require some memory if nkibz is large. 1062 ABI_MALLOC(Wfd%Kdata, (Wfd%nkibz)) 1063 Wfd%Kdata%use_fnl_dir0der0 = Wfd%use_fnl_dir0der0 1064 1065 do ik_ibz=1,Wfd%nkibz 1066 kpoint = Wfd%kibz(:,ik_ibz) 1067 istwf_k = Wfd%istwfk(ik_ibz) 1068 npw_k = Wfd%npwarr(ik_ibz) 1069 if (any(wfd%bks2wfd(1, :, ik_ibz, :) /= 0)) then 1070 call Wfd%Kdata(ik_ibz)%init(Cryst, Psps, kpoint, istwf_k, ngfft, Wfd%MPI_enreg, ecut=Wfd%ecut) 1071 end if 1072 end do 1073 1074 select type (wfd) 1075 class is (wfdgw_t) 1076 ! Allocate the global table used to keep track of the bks distribution, including possible duplication. 1077 bks_size = one * wfd%mband * wfd%nkibz * wfd%nsppol * wfd%nproc 1078 write(msg,'(a,f8.1,a)')' Memory needed for bks_tab: ',one * bks_size * b2Mb,' [Mb] <<< MEM' 1079 call wrtout(std_out, msg) 1080 1081 !ABI_MALLOC(wfd%bks_ranks, (wfd%mband, nkibz, nsppol)) 1082 1083 ABI_MALLOC(Wfd%bks_tab, (Wfd%mband, nkibz, nsppol, 0:Wfd%nproc-1)) 1084 Wfd%bks_tab = WFD_NOWAVE 1085 1086 ! Update the kbs table storing the distribution of the ug. 1087 call wfd%update_bkstab(show=-std_out) 1088 end select 1089 1090 call cwtime_report(" wfd_init", cpu, wall, gflops) 1091 1092 end subroutine wfd_init
m_wfd/wfd_mybands [ Functions ]
NAME
wfd_mybands
FUNCTION
Return the list of band indices of the ug owned by this node at given (k,s).
INPUTS
ik_ibz=Index of the k-point in the IBZ
spin=spin index
[how]=string defining which status is checked.
Possible mutually exclusive values: "Allocated", "Stored".
Only the first character is checked (no case-sensitive)
By default the list of bands whose status is either WFD_ALLOCATED or WFD_STORED is returned.
OUTPUT
how_manyb=The number of bands owned by this node my_band_list(Wfd%mband)=The first how_manyb values are the bands treated by this node.
SOURCE
2461 subroutine wfd_mybands(Wfd, ik_ibz, spin, how_manyb, my_band_list, how) 2462 2463 !Arguments ------------------------------------ 2464 !scalars 2465 integer,intent(in) :: ik_ibz,spin 2466 integer,intent(out) :: how_manyb 2467 character(len=*),optional,intent(in) :: how 2468 class(wfd_t),intent(in) :: Wfd 2469 !arrays 2470 integer,intent(out) :: my_band_list(Wfd%mband) 2471 2472 !Local variables ------------------------------ 2473 !scalars 2474 integer :: band 2475 logical :: do_have 2476 2477 !************************************************************************ 2478 2479 how_manyb=0; my_band_list=-1 2480 do band=1,Wfd%nband(ik_ibz,spin) 2481 if (present(how)) then 2482 do_have = wfd%ihave_ug(band, ik_ibz, spin, how=how) 2483 else 2484 do_have = wfd%ihave_ug(band, ik_ibz, spin) 2485 end if 2486 if (do_have) then 2487 how_manyb = how_manyb + 1 2488 my_band_list(how_manyb) = band 2489 end if 2490 end do 2491 2492 end subroutine wfd_mybands
m_wfd/wfd_norm2 [ Functions ]
NAME
wfd_norm2
FUNCTION
Compute <u_{bks}|u_{bks}> in G-space
INPUTS
Wfd<wfd_t>=the wavefunction descriptor. Cryst<crystal_t>=Structure describing the crystal structure and its symmetries. Pawtab(ntypat*usepaw)<type(pawtab_type)>=PAW tabulated starting data. band=Band index. ik_bz=Index of the k-point in the BZ. spin=Spin index
SOURCE
1293 function wfd_norm2(Wfd,Cryst,Pawtab,band,ik_ibz,spin) result(norm2) 1294 1295 !Arguments ------------------------------------ 1296 !scalars 1297 integer,intent(in) :: band,ik_ibz,spin 1298 real(dp) :: norm2 1299 type(crystal_t),intent(in) :: Cryst 1300 class(wfd_t),target,intent(inout) :: Wfd 1301 type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw) 1302 1303 !Local variables ------------------------------ 1304 !scalars 1305 integer :: npw_k,istwf_k 1306 complex(dpc) :: cdum 1307 type(wave_t),pointer :: wave 1308 character(len=500) :: msg 1309 !arrays 1310 real(dp) :: pawovlp(2) 1311 complex(gwpc),ABI_CONTIGUOUS pointer :: ug1(:) 1312 type(pawcprj_type),allocatable :: Cp1(:,:) 1313 1314 !************************************************************************ 1315 1316 ! Planewave part. 1317 npw_k = Wfd%npwarr(ik_ibz) 1318 istwf_k = Wfd%istwfk(ik_ibz) 1319 1320 ABI_CHECK(wfd%get_wave_ptr(band, ik_ibz, spin, wave, msg) == 0, msg) 1321 1322 ug1 => wave%ug 1323 cdum = xdotc(Wfd%nspinor*npw_k,ug1,1,ug1,1) 1324 1325 if (istwf_k>1) then 1326 cdum=two*DBLE(cdum) 1327 if (istwf_k==2) cdum=cdum-CONJG(ug1(1))*ug1(1) 1328 end if 1329 1330 ! Paw on-site term. 1331 if (Wfd%usepaw==1) then 1332 1333 ! Avoid the computation if Cprj are already in memory with the correct order. 1334 if (wave%has_cprj == WFD_STORED .and. wave%cprj_order == CPR_RANDOM) then 1335 1336 pawovlp = paw_overlap(wave%Cprj, wave%Cprj, Cryst%typat, Pawtab) 1337 cdum = cdum + CMPLX(pawovlp(1),pawovlp(2), kind=dpc) 1338 1339 else 1340 ! Compute Cproj 1341 ABI_MALLOC(Cp1,(Wfd%natom,Wfd%nspinor)) 1342 call pawcprj_alloc(Cp1,0,Wfd%nlmn_atm) 1343 1344 call wfd%get_cprj(band,ik_ibz,spin,Cryst,Cp1,sorted=.FALSE.) 1345 pawovlp = paw_overlap(Cp1,Cp1,Cryst%typat,Pawtab) 1346 cdum = cdum + CMPLX(pawovlp(1),pawovlp(2), kind=dpc) 1347 1348 call pawcprj_free(Cp1) 1349 ABI_FREE(Cp1) 1350 end if 1351 end if 1352 1353 norm2 = DBLE(cdum) 1354 1355 end function wfd_norm2
m_wfd/wfd_paw_get_aeur [ Functions ]
NAME
wfd_paw_get_aeur
FUNCTION
Compute the AE PAW wavefunction in real space.
INPUTS
band,ik_ibz,spin=indices specifying the band, the k-point and the spin. Psps<pseudopotential_type>=variables related to pseudopotentials Cryst<crystal_t>= data type gathering info on symmetries and unit cell. Wfd<wfd_t>=wavefunction descriptor. Pawtab(ntypat*usepaw)<type(pawtab_type)>=paw tabulated starting data. Pawfgrtab(natom)<pawfgrtab_type>= atomic data given on fine rectangular grid. NB: rpaw should be used in nhatgrid to initialize the datatype (optcut=1 option) instead of the radius for the shape functions (rpaw /= rshp). Paw_onsite(natom)<paw_pwaves_lmn_t>=3D PAW partial waves in real space for each FFT point in the PAW spheres.
OUTPUT
ur_ae(Wfd%nfft*Wfd%nspinor)=AE PAW wavefunction in real space. [ur_ae_onsite(Wfd%nfft*Wfd%nspinor)] [ur_ps_onsite(Wfd%nfft*Wfd%nspinor)]
NOTES
(1) The true wavefunction integrates in real space to the unit cell volume. The definition of the cprj matrix elements includes the term 1/SQRT(ucvol) that comes from the use of a normalized planewave e^(iG.r)/SQRT(omega) in the FFT transform G-->R (see e.g. opernla_ylm) On the contrary, the convention for the G-->R transform employed in the FFT routines used in abinit is u(r) = sum_G u(G) e^(iG.r); u(G) = one/omega \int u(r) e^(-iG.r)dr. Hence we have to multiply the onsite part by SQRT(uvol) before adding the smooth FFT part in real space. (2) Care has to be taken in the calculation of the onsite contribution when the FFT point belongs to the PAW sphere of a periodically repeated atom. In this case one evaluates the onsite term associated to the atom in the first unit cell then the contribution has to be multiplied by a k- dependent phase factor to account for the wrapping of the real-space point in the first unit cell.
SOURCE
4810 subroutine wfd_paw_get_aeur(Wfd,band,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae,ur_ae_onsite,ur_ps_onsite) 4811 4812 !Arguments ------------------------------------ 4813 !scalars 4814 integer,intent(in) :: band,ik_ibz,spin 4815 type(pseudopotential_type),intent(in) :: Psps 4816 type(crystal_t),intent(in) :: Cryst 4817 class(wfd_t),intent(inout) :: Wfd 4818 !arrays 4819 type(pawtab_type),intent(in) :: Pawtab(Cryst%ntypat) 4820 type(pawfgrtab_type),intent(in) :: Pawfgrtab(Cryst%natom) 4821 type(paw_pwaves_lmn_t),intent(in) :: Paw_onsite(Cryst%natom) 4822 complex(gwpc),intent(out) :: ur_ae(Wfd%nfft*Wfd%nspinor) 4823 complex(gwpc),optional,intent(out) :: ur_ae_onsite(Wfd%nfft*Wfd%nspinor) 4824 complex(gwpc),optional,intent(out) :: ur_ps_onsite(Wfd%nfft*Wfd%nspinor) 4825 4826 !Local variables------------------------------- 4827 !scalars 4828 integer :: itypat,ln_size,lmn_size,iatom,spinor 4829 integer :: nfgd,ifgd,jlmn,jl,jm,ifftsph 4830 real(dp) :: phj,tphj,arg,re_cp,im_cp 4831 complex(dpc) :: cp,cnorm 4832 !arrays 4833 real(dp) :: kpoint(3) 4834 complex(dpc),allocatable :: ceikr(:),phk_atm(:) 4835 type(pawcprj_type),allocatable :: Cp1(:,:) 4836 4837 ! ************************************************************************* 4838 4839 ! TODO ngfft should be included in pawfgrtab_type 4840 !% if (ANY(Wfd%ngfft(1:3)/=Pawfgrtab%ngfft(1:3)) then 4841 !! ABI_ERROR("Wfd%ngfft(1:3)/=Pawfgrtab%ngfft(1:3)") 4842 !% end if 4843 4844 call wfd%get_ur(band,ik_ibz,spin,ur_ae) 4845 4846 kpoint = Wfd%kibz(:,ik_ibz) 4847 ABI_MALLOC(ceikr, (Wfd%nfftot * wfd%nspinor)) 4848 4849 call calc_ceikr(kpoint, wfd%ngfft, Wfd%nfftot, wfd%nspinor, ceikr) 4850 ur_ae = ur_ae * ceikr 4851 4852 ABI_MALLOC(Cp1,(Wfd%natom,Wfd%nspinor)) 4853 call pawcprj_alloc(Cp1,0,Wfd%nlmn_atm) 4854 4855 call wfd%get_cprj(band,ik_ibz,spin,Cryst,Cp1,sorted=.FALSE.) 4856 4857 ! Add onsite term on the augmented FFT mesh. 4858 if (present(ur_ae_onsite)) ur_ae_onsite = czero 4859 if (present(ur_ps_onsite)) ur_ps_onsite = czero 4860 4861 ABI_CHECK(Wfd%nspinor==1,"nspinor==1 not coded") 4862 4863 do iatom=1,Cryst%natom 4864 itypat =Cryst%typat(iatom) 4865 lmn_size=Pawtab(itypat)%lmn_size 4866 ln_size =Pawtab(itypat)%basis_size ! no. of nl elements in PAW basis. 4867 nfgd =Pawfgrtab(iatom)%nfgd ! no. of points in the fine grid for this PAW sphere. 4868 4869 ABI_MALLOC(phk_atm,(nfgd)) 4870 do ifgd=1,nfgd 4871 arg = -two_pi* DOT_PRODUCT(Paw_onsite(iatom)%r0shift(:,ifgd),kpoint) 4872 phk_atm(ifgd) = DCMPLX(COS(arg),SIN(arg)) 4873 end do 4874 4875 do spinor=1,Wfd%nspinor 4876 do jlmn=1,lmn_size 4877 jl=Psps%indlmn(1,jlmn,itypat) 4878 jm=Psps%indlmn(2,jlmn,itypat) 4879 re_cp = Cp1(iatom,spinor)%cp(1,jlmn) 4880 im_cp = Cp1(iatom,spinor)%cp(2,jlmn) 4881 cp = DCMPLX(re_cp, im_cp) * SQRT(Cryst%ucvol) ! Pay attention here. see (1). 4882 4883 do ifgd=1,nfgd ! loop over fine grid points in current PAW sphere. 4884 ifftsph = Pawfgrtab(iatom)%ifftsph(ifgd) ! index of the point on the grid 4885 phj = Paw_onsite(iatom)% phi(ifgd,jlmn) 4886 tphj = Paw_onsite(iatom)%tphi(ifgd,jlmn) 4887 ur_ae(ifftsph) = ur_ae(ifftsph) + cp * (phj-tphj) * phk_atm(ifgd) 4888 if (present(ur_ae_onsite)) ur_ae_onsite(ifftsph) = ur_ae_onsite(ifftsph) + cp * phj * phk_atm(ifgd) 4889 if (present(ur_ps_onsite)) ur_ps_onsite(ifftsph) = ur_ps_onsite(ifftsph) + cp * tphj * phk_atm(ifgd) 4890 end do 4891 end do !jlmn 4892 end do !spinor 4893 4894 ABI_FREE(phk_atm) 4895 end do !iatom 4896 4897 ! Remove the phase e^{ikr}, u(r) is returned. 4898 ur_ae = ur_ae * CONJG(ceikr) 4899 cnorm = xdotc(Wfd%nfft*Wfd%nspinor,ur_ae,1,ur_ae,1)/Wfd%nfft 4900 !write(std_out,*)" AE PAW norm: (b,k,s)",band,ik_ibz,spin,REAL(cnorm) 4901 4902 if (present(ur_ae_onsite)) ur_ae_onsite = ur_ae_onsite * CONJG(ceikr) 4903 if (present(ur_ps_onsite)) ur_ps_onsite = ur_ps_onsite * CONJG(ceikr) 4904 4905 call pawcprj_free(Cp1) 4906 ABI_FREE(Cp1) 4907 ABI_FREE(ceikr) 4908 4909 end subroutine wfd_paw_get_aeur
m_wfd/wfd_print [ Functions ]
NAME
wfd_print
FUNCTION
Print the content of a wfd_t datatype
INPUTS
Wfd<wfd_t>=The datatype. [header]=String to be printed as header for additional info. [unit]=Unit number for output [prtvol]=Verbosity level [mode_paral]=Either "COLL" or "PERS". Defaults to "COLL".
OUTPUT
Only printing
SOURCE
1689 subroutine wfd_print(Wfd, header, unit, prtvol, mode_paral) 1690 1691 !Arguments ------------------------------------ 1692 integer,optional,intent(in) :: unit,prtvol 1693 character(len=4),optional,intent(in) :: mode_paral 1694 character(len=*),optional,intent(in) :: header 1695 class(wfd_t),intent(in) :: Wfd 1696 1697 !Local variables------------------------------- 1698 !scalars 1699 integer :: my_prtvol, my_unt, mpw, ib, ik, is, ug_cnt, ur_cnt, cprj_cnt, spin, ik_ibz, band 1700 real(dp) :: ug_size, ur_size, cprj_size !,kdata_bsize 1701 character(len=4) :: my_mode 1702 character(len=500) :: msg 1703 ! ************************************************************************* 1704 1705 my_unt =std_out; if (present(unit )) my_unt =unit 1706 my_prtvol=0 ; if (present(prtvol )) my_prtvol=prtvol 1707 my_mode ='COLL' ; if (present(mode_paral)) my_mode =mode_paral 1708 1709 msg=' ==== Info on the Wfd% object ==== ' 1710 if (present(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== ' 1711 call wrtout(my_unt,msg,my_mode) 1712 1713 write(msg,'(3(a,i0,a),a,i0,2a,f5.1)')& 1714 ' Number of irreducible k-points ........ ',Wfd%nkibz,ch10,& 1715 ' Number of spinorial components ........ ',Wfd%nspinor,ch10,& 1716 ' Number of spin-density components ..... ',Wfd%nspden,ch10,& 1717 ' Number of spin polarizations .......... ',Wfd%nsppol,ch10,& 1718 ' Plane wave cutoff energy .............. ',Wfd%ecut 1719 call wrtout(my_unt, msg, my_mode) 1720 1721 mpw = maxval(Wfd%npwarr) 1722 write(msg,'(3(a,i0,a))')& 1723 ' Max number of G-vectors ............... ',mpw,ch10,& 1724 ' Total number of FFT points ............ ',Wfd%nfftot,ch10,& 1725 ' Number of FFT points treated by me .... ',Wfd%nfft,ch10 1726 call wrtout(my_unt, msg, my_mode) 1727 1728 call print_ngfft(Wfd%ngfft, 'FFT mesh for wavefunctions', my_unt, my_mode, my_prtvol) 1729 1730 ug_cnt = 0; ur_cnt = 0; cprj_cnt = 0 1731 do spin=1,Wfd%nsppol 1732 do ik_ibz=1,Wfd%nkibz 1733 do band=1,Wfd%nband(ik_ibz,spin) 1734 ib = wfd%bks2wfd(1, band, ik_ibz, spin) 1735 if (ib == 0) cycle 1736 ik = wfd%bks2wfd(2, band, ik_ibz, spin) 1737 is = wfd%bks2wfd(3, band, ik_ibz, spin) 1738 if (wfd%s(is)%k(ik)%b(ib)%has_ug >= WFD_ALLOCATED) ug_cnt = ug_cnt + 1 1739 if (wfd%s(is)%k(ik)%b(ib)%has_ur >= WFD_ALLOCATED) ur_cnt = ur_cnt + 1 1740 if (wfd%s(is)%k(ik)%b(ib)%has_cprj >= WFD_ALLOCATED) cprj_cnt = cprj_cnt + 1 1741 end do 1742 end do 1743 end do 1744 1745 ! Info on memory needed for u(g), u(r) and PAW cprj 1746 write(msg, '(a,i0)')' Total number of (b,k,s) states stored by this rank: ', ug_cnt 1747 call wrtout(std_out, msg, pre_newlines=1) 1748 1749 ug_size = one * Wfd%nspinor * mpw * ug_cnt 1750 write(msg,'(a,f8.1,a)')' Memory allocated for Fourier components u(G): ',two*gwpc*ug_size*b2Mb,' [Mb] <<< MEM' 1751 call wrtout(std_out, msg) 1752 1753 if (any(wfd%keep_ur)) then 1754 ur_size = one * Wfd%nspinor * Wfd%nfft * ur_cnt 1755 write(msg,'(a,f8.1,a)')' Memory allocated for real-space u(r): ',two*gwpc*ur_size*b2Mb,' [Mb] <<< MEM' 1756 call wrtout(std_out, msg) 1757 end if 1758 1759 if (wfd%usepaw==1) then 1760 cprj_size = one * Wfd%nspinor * sum(Wfd%nlmn_atm) * cprj_cnt 1761 write(msg,'(a,f8.1,a)')' Memory allocated for PAW projections cprj: ',dp*cprj_size*b2Mb,' [Mb] <<< MEM' 1762 call wrtout(std_out, msg) 1763 end if 1764 1765 !TODO 1766 ! Add additionanl info 1767 !kdata_bsize = nkibz * (four * (3 * mpw) + dp * two * mpw * natom) 1768 !write(msg,'(a,f8.1,a)')' Memory allocated for Kdata = ',kdata_bsize * b2Mb,' [Mb] <<< MEM' 1769 1770 write(msg,'(a,f8.1,a)')' Memory needed for wfd%s datastructure: ',ABI_MEM_MB(wfd%s),' [Mb] <<< MEM' 1771 call wrtout(std_out, msg) 1772 write(msg,'(a,f8.1,a)')' Memory needed for wfd%s(0)%k datastructure: ',ABI_MEM_MB(wfd%s(1)%k),' [Mb] <<< MEM' 1773 call wrtout(std_out, msg) 1774 write(msg,'(a,f8.1,a)')' Memory allocated for Kdata array: ',ABI_MEM_MB(wfd%kdata),' [Mb] <<< MEM' 1775 call wrtout(std_out, msg, newlines=1) 1776 1777 end subroutine wfd_print
m_wfd/wfd_push_ug [ Functions ]
NAME
wfd_push_ug
FUNCTION
This routine changes the status of the object by saving the wavefunction in the correct slot inside Wfd%Wave. It also set the corresponding has_ug flag to WFD_STORED. If the status of the corresponding ur is (WFD_STORED|WFD_ALLOCATED), then an G->R FFT transform is done (see also update_ur)
INPUTS
band=Band index. ik_ibz=k-point index spin=Spin index. Cryst<crystal_t>=Object defining the unit cell and its symmetries. ug(npw_k*Wfd%nspinor)=The ug to be saved. [update_ur]=If .FALSE.: no G-->R transform is done even if ur is (WFD_STORED|WFD_ALLOCATED) so be careful. Defaults to .TRUE. [update_cprj]=If .FALSE.: <C|p_i> matrix elements are not recalculatedd even if cprj is (WFD_STORED|WFD_ALLOCATED) so be careful. Defaults to .TRUE.
SIDE EFFECTS
Wfd<wfd_t>=See above.
SOURCE
2176 subroutine wfd_push_ug(Wfd, band, ik_ibz, spin, Cryst, ug, update_ur, update_cprj) 2177 2178 !Arguments ------------------------------------ 2179 !scalars 2180 integer,intent(in) :: ik_ibz,spin,band 2181 logical,optional,intent(in) :: update_ur,update_cprj 2182 class(wfd_t),target,intent(inout) :: Wfd 2183 type(crystal_t),intent(in) :: Cryst 2184 !arrays 2185 complex(gwpc),intent(inout) :: ug(:) 2186 2187 !Local variables ------------------------------ 2188 !scalars 2189 integer,parameter :: choice1=1,idir0=0,tim_fourdp=5,ndat1=1 2190 integer :: npw_k, ib, ik, is 2191 logical :: do_update_ur,do_update_cprj,want_sorted 2192 character(len=500) :: msg 2193 type(wave_t),pointer :: wave 2194 2195 !************************************************************************ 2196 2197 if (size(ug) /= Wfd%npwarr(ik_ibz) * Wfd%nspinor) then 2198 ABI_ERROR("Wrong size in assumed shape array") 2199 end if 2200 2201 if (any(wfd%bks2wfd(:, band, ik_ibz, spin) == 0)) then 2202 write(msg,'(a,i0,a,3(i0,1x))')" MPI rank ",Wfd%my_rank," doesn't have ug for (band, ik_ibz, spin): ",band,ik_ibz,spin 2203 ABI_ERROR(msg) 2204 end if 2205 2206 ib = wfd%bks2wfd(1, band, ik_ibz, spin) 2207 ik = wfd%bks2wfd(2, band, ik_ibz, spin) 2208 is = wfd%bks2wfd(3, band, ik_ibz, spin) 2209 2210 wave => wfd%s(is)%k(ik)%b(ib) 2211 wave%ug = ug 2212 wave%has_ug = WFD_STORED 2213 2214 if (Wfd%debug_level>0) then 2215 if (wave%has_ug == WFD_NOWAVE) then 2216 write(msg,'(a,i0,a,3(i0,1x))')" MPI rank ",Wfd%my_rank," doesn't have ug for (band, ik_ibz, spin): ",band,ik_ibz,spin 2217 ABI_ERROR(msg) 2218 end if 2219 end if 2220 2221 if (Wfd%usepaw==1) then 2222 ! Update the corresponding cprj if required. 2223 do_update_cprj=.TRUE.; if (present(update_cprj)) do_update_cprj=update_cprj 2224 if (do_update_cprj) then 2225 want_sorted = (wave%cprj_order == CPR_SORTED) 2226 call wfd%ug2cprj(band, ik_ibz, spin, choice1, idir0, wfd%natom, cryst, wave%cprj, sorted=want_sorted) 2227 wave%has_cprj = WFD_STORED 2228 else 2229 wave%has_cprj = WFD_ALLOCATED 2230 end if 2231 end if 2232 2233 if (any(wave%has_ur == [WFD_STORED, WFD_ALLOCATED])) then 2234 ! Update the corresponding ur if required. 2235 do_update_ur=.TRUE.; if (present(update_ur)) do_update_ur=update_ur 2236 2237 if (do_update_ur) then 2238 npw_k = Wfd%npwarr(ik_ibz) 2239 call fft_ug(npw_k,Wfd%nfft,Wfd%nspinor,ndat1,Wfd%mgfft,Wfd%ngfft,Wfd%istwfk(ik_ibz),& 2240 Wfd%Kdata(ik_ibz)%kg_k,Wfd%Kdata(ik_ibz)%gbound,ug,wave%ur) 2241 wave%has_ur = WFD_STORED 2242 else 2243 wave%has_ur = WFD_ALLOCATED 2244 end if 2245 end if 2246 2247 end subroutine wfd_push_ug
m_wfd/wfd_read_wfk [ Functions ]
NAME
wfd_read_wfk
FUNCTION
This routine reads the WFK file completing the initialization of the wavefunction descriptor
INPUTS
wfk_fname=Name of the WFK file. iomode=Option specifying the fileformat as well as the IO mode to be used.
OUTPUT
[out_hdr]=Header of the WFK file.
SIDE EFFECTS
Wfd<wfd_t>=All the states owned by this node whose status is (STORED|ALLOCATED) read.
SOURCE
4410 subroutine wfd_read_wfk(Wfd, wfk_fname, iomode, out_hdr) 4411 4412 !Arguments ------------------------------------ 4413 !scalars 4414 integer,intent(in) :: iomode 4415 character(len=*),intent(in) :: wfk_fname 4416 class(wfd_t),target,intent(inout) :: Wfd 4417 type(Hdr_type),optional,intent(inout) :: out_hdr ! ifort and others are buggy for optional intent(out) structured types 4418 4419 !Local variables ------------------------------ 4420 !scalars 4421 integer,parameter :: formeig0 = 0, optkg1 = 1, method = 2 4422 integer :: wfk_unt,npw_disk,nmiss,ig,sc_mode,ii 4423 integer :: io_comm,master,my_rank,spin,ik_ibz,fform,ierr ! ,igp 4424 integer :: mcg,nband_wfd,nband_disk,band,mband_disk,bcount,istwfk_disk 4425 integer :: spinor,cg_spad,gw_spad,icg,igw,cg_bpad, allcg_bpad, ib, ik, is 4426 integer :: my_bmin, my_bmax, bmin, bmax 4427 logical :: change_gsphere, master_only, iread 4428 real(dp) :: cpu, wall, gflops, cpu_ks, wall_ks, gflops_ks 4429 character(len=500) :: msg 4430 type(Wfk_t) :: Wfk 4431 type(Hdr_type) :: Hdr 4432 type(wave_t),pointer :: wave 4433 !arrays 4434 integer,allocatable :: gf2wfd(:), kg_k(:,:), all_countks(:,:) 4435 integer :: work_ngfft(18),gmax_wfd(3),gmax_disk(3),gmax(3) 4436 real(dp) :: tsec(2) 4437 real(dp),allocatable :: eig_k(:), cg_k(:,:), out_cg(:,:), work(:,:,:,:), allcg_k(:,:) 4438 logical,allocatable :: my_readmask(:,:,:) 4439 character(len=6) :: tag_spin(2) 4440 4441 !************************************************************************ 4442 4443 ! Keep track of time spent in wfd_read_wfk 4444 call timab(300, 1, tsec) 4445 4446 if (any(iomode == [IO_MODE_NETCDF, IO_MODE_FORTRAN_MASTER])) then 4447 ABI_ERROR(sjoin("Unsupported value for iomode: ", itoa(iomode))) 4448 end if 4449 4450 ! IO_MODE_FORTRAN --> only master reads and broadcasts data. 4451 ! IO_MODE_MPI --> all procs read with collective operations. 4452 my_rank = Wfd%my_rank; master = Wfd%master 4453 io_comm = wfd%comm; sc_mode = xmpio_collective; master_only = .False.; iread = .True. 4454 !if (iomode == IO_MODE_FORTRAN) then 4455 io_comm = xmpi_comm_self; sc_mode = xmpio_single; master_only = .True.; iread = my_rank == wfd%master 4456 !end if 4457 4458 call wrtout(std_out, sjoin( & 4459 " wfd_read_wfk: Reading file:", wfk_fname, & 4460 " with iomode:", iomode2str(iomode),", master_only:", yesno(master_only)), pre_newlines=2) 4461 if (iomode == IO_MODE_MPI) then 4462 call wrtout(std_out, sjoin( & 4463 " If MPI-IO is too slow, use the command line option `abinit --enforce-fortran-io ...`", ch10, & 4464 " to make the master proc read data with Fortran-IO and then broadcast (requires more memory)"), do_flush=.True.) 4465 end if 4466 4467 if (iread) then 4468 wfk_unt = get_unit() 4469 call wfk_open_read(Wfk, wfk_fname, formeig0, iomode, wfk_unt, io_comm, Hdr_out=Hdr) 4470 end if 4471 4472 if (master_only) call hdr%bcast(wfd%master, wfd%my_rank, wfd%comm) 4473 if (present(out_hdr)) call hdr_copy(hdr, out_hdr) 4474 4475 ! TODO: Perform more consistency check btw Hdr and Wfd. 4476 ! Output the header of the GS wavefunction file. 4477 fform = 0 4478 if (wfd%prtvol /= 0 .and. wfd%my_rank == 0) call hdr%echo(fform, 4, unit=std_out) 4479 4480 mband_disk = MAXVAL(Hdr%nband) 4481 ABI_CHECK_ILEQ(Wfd%mband, mband_disk, "Not enough bands stored on WFK file") 4482 4483 ! Each node will read the waves whose status if (WFD_ALLOCATED|WFD_STORED). 4484 ! all_countks is a global array used to skip (ik_ibz, spin) if all MPI procs do not need bands for this (k, s) 4485 ABI_MALLOC(my_readmask, (mband_disk, Wfd%nkibz, Wfd%nsppol)) 4486 my_readmask = .False. 4487 my_bmin = huge(1); my_bmax = -huge(1) 4488 ABI_ICALLOC(all_countks, (wfd%nkibz, wfd%nsppol)) 4489 4490 do spin=1,Wfd%nsppol 4491 do ik_ibz=1,Wfd%nkibz 4492 do band=1,Wfd%nband(ik_ibz,spin) 4493 if (wfd%ihave_ug(band, ik_ibz, spin)) then 4494 my_bmin = min(my_bmin, band) 4495 my_bmax = max(my_bmax, band) 4496 my_readmask(band, ik_ibz, spin) = .True. 4497 all_countks(ik_ibz, spin) = 1 4498 if (wfd%ihave_ug(band, ik_ibz, spin, how="Stored")) then 4499 ABI_ERROR("Wavefunction is already stored!") 4500 end if 4501 end if 4502 end do 4503 end do 4504 end do 4505 4506 ! All procs must agree when skipping (k, s) states 4507 ! We also need bmin/bmax for master_only option. 4508 call xmpi_sum(all_countks, wfd%comm, ierr) 4509 call xmpi_min(my_bmin, bmin, wfd%comm, ierr) 4510 call xmpi_max(my_bmax, bmax, wfd%comm, ierr) 4511 4512 call wrtout(std_out, sjoin(" About to read: ",itoa(count(my_readmask)), " (b, k, s) states in total.")) 4513 do spin=1,wfd%nsppol 4514 call wrtout(std_out, sjoin(" For spin:", itoa(spin), & 4515 ", will read:", itoa(count(any(my_readmask(:,:,spin), dim=1))), " k-points out of:", itoa(wfd%nkibz))) 4516 end do 4517 tag_spin(: )= [' ',' ']; if (Wfd%nsppol==2) tag_spin(:)= [' UP ',' DOWN '] 4518 if (wfd%prtvol > 0) call wrtout(std_out,' k eigenvalues [eV]') 4519 4520 call cwtime(cpu, wall, gflops, "start") 4521 4522 if (method == 1) then 4523 do spin=1,wfd%nsppol 4524 do ik_ibz=1,Wfd%nkibz 4525 if (all_countks(ik_ibz, spin) == 0) cycle 4526 npw_disk = Hdr%npwarr(ik_ibz) 4527 nband_disk = Hdr%nband(ik_ibz+(spin-1)*Hdr%nkpt) 4528 istwfk_disk = hdr%istwfk(ik_ibz) 4529 change_gsphere = istwfk_disk /= wfd%istwfk(ik_ibz) 4530 ABI_CHECK(.not. change_gsphere, "different istwfk values are not coded") 4531 4532 nband_wfd = Wfd%nband(ik_ibz,spin) 4533 if (nband_wfd > nband_disk) then 4534 write(msg,'(a,2(i0,1x))')" nband_wfd to be read cannot be greater than nband_disk while: ",nband_wfd, nband_disk 4535 ABI_ERROR(msg) 4536 end if 4537 4538 mcg = npw_disk*Wfd%nspinor*nband_wfd 4539 4540 ABI_MALLOC(eig_k,((2*Wfk%mband)**formeig0*Wfk%mband)) 4541 4542 ABI_MALLOC(kg_k,(3,optkg1*npw_disk)) 4543 ABI_MALLOC_OR_DIE(cg_k,(2,mcg), ierr) 4544 4545 call wfk%read_band_block([1, nband_wfd], ik_ibz, spin, sc_mode, kg_k=kg_k, cg_k=cg_k, eig_k=eig_k) 4546 4547 if (wfd%prtvol > 0 .and. Wfd%my_rank==Wfd%master) then 4548 if (Wfd%nsppol==2) then 4549 write(std_out,'(i3,a,10f7.2/50(10x,10f7.2/))') ik_ibz,tag_spin(spin),(eig_k(ib)*Ha_eV,ib=1,nband_wfd) 4550 else 4551 write(std_out,'(i3,7x,10f7.2/50(10x,10f7.2/))')ik_ibz,(eig_k(ib)*Ha_eV,ib=1,nband_wfd) 4552 end if 4553 end if 4554 4555 ! Table with the correspondence btw the k-centered sphere of the WFK file 4556 ! and the one used in Wfd (possibly smaller due to ecutwfn). 4557 ABI_MALLOC(gf2wfd,(npw_disk)) 4558 if (any(my_readmask(:,ik_ibz,spin))) then 4559 call kg_map(wfd%npwarr(ik_ibz), wfd%kdata(ik_ibz)%kg_k, npw_disk, kg_k, gf2wfd, nmiss) 4560 end if 4561 !if (nmiss/=0) then 4562 ! write(msg,'(a,2(1x,i0),a,i0)')" For (k,s) ",ik_ibz,spin," the number of missing G is ",nmiss 4563 ! ABI_WARNING(msg) 4564 !end if 4565 4566 ! Conversion of the basis set. 4567 do band=1,Wfd%nband(ik_ibz,spin) 4568 if (my_readmask(band, ik_ibz, spin)) then 4569 4570 ABI_CHECK(all(wfd%bks2wfd(:, band, ik_ibz, spin) /= 0), "state in not allocated") 4571 ib = wfd%bks2wfd(1, band, ik_ibz, spin) 4572 ik = wfd%bks2wfd(2, band, ik_ibz, spin) 4573 is = wfd%bks2wfd(3, band, ik_ibz, spin) 4574 wave => wfd%s(is)%k(ik)%b(ib) 4575 wave%ug = czero 4576 4577 cg_bpad=npw_disk*Wfd%nspinor*(band-1) 4578 do spinor=1,Wfd%nspinor 4579 cg_spad=(spinor-1)*npw_disk 4580 gw_spad=(spinor-1)*Wfd%npwarr(ik_ibz) 4581 do ig=1,npw_disk 4582 icg = ig+cg_spad+cg_bpad 4583 igw = gf2wfd(ig)+gw_spad 4584 if (gf2wfd(ig) /= 0) then 4585 wave%ug(igw) = CMPLX(cg_k(1,icg), cg_k(2,icg), kind=gwpc) 4586 end if 4587 end do 4588 end do 4589 wave%has_ug = WFD_STORED 4590 4591 end if 4592 end do 4593 4594 ABI_FREE(eig_k) 4595 ABI_FREE(kg_k) 4596 ABI_FREE(cg_k) 4597 ABI_FREE(gf2wfd) 4598 end do !ik_ibz 4599 end do !spin 4600 4601 else if (method==2) then 4602 ! DEFAULT ALGO: This seems to be the most efficient one. 4603 4604 do spin=1,Wfd%nsppol 4605 do ik_ibz=1,Wfd%nkibz 4606 if (all_countks(ik_ibz, spin) == 0) cycle 4607 call cwtime(cpu_ks, wall_ks, gflops_ks, "start") 4608 !write(std_out,*)" About to read ik_ibz: ",ik_ibz,", spin: ",spin 4609 4610 npw_disk = Hdr%npwarr(ik_ibz) 4611 nband_disk = Hdr%nband(ik_ibz+(spin-1)*Hdr%nkpt) 4612 istwfk_disk = hdr%istwfk(ik_ibz) 4613 change_gsphere = istwfk_disk /= wfd%istwfk(ik_ibz) 4614 nband_wfd = Wfd%nband(ik_ibz, spin) 4615 4616 if (nband_wfd > nband_disk) then 4617 write(msg,'(2(a, i0))')& 4618 "nband_wfd to be read: ", nband_wfd,", cannot be greater than nband_disk: ",nband_disk 4619 ABI_ERROR(msg) 4620 end if 4621 4622 ! Allocate full array for eigenvalues and G-vectors. 4623 ABI_MALLOC(eig_k, ((2*nband_disk)**formeig0*nband_disk)) 4624 ABI_MALLOC(kg_k, (3,optkg1*npw_disk)) 4625 4626 ! Allocate array to store my wavefunctions and read data 4627 mcg = npw_disk * wfd%nspinor * count(my_readmask(:, ik_ibz, spin)) 4628 ABI_MALLOC_OR_DIE(cg_k, (2, mcg), ierr) 4629 4630 if (.not. master_only) then 4631 ! All procs read. 4632 call wfk%read_bmask(my_readmask(:,ik_ibz, spin), ik_ibz, spin, sc_mode, kg_k=kg_k, cg_k=cg_k, eig_k=eig_k) 4633 4634 else 4635 ! Master reads full set of bands and broadcasts data, then each proc extract its own set of wavefunctions. 4636 ! TODO: Should read in blocks to reduce memory footprint 4637 ABI_MALLOC_OR_DIE(allcg_k, (2, npw_disk*wfd%nspinor*(bmax-bmin+1)), ierr) 4638 if (my_rank == master) then 4639 call wfk%read_band_block([bmin, bmax], ik_ibz, spin, xmpio_single, kg_k=kg_k, cg_k=allcg_k, eig_k=eig_k) 4640 end if 4641 call xmpi_bcast(kg_k, wfd%master, wfd%comm, ierr) 4642 call xmpi_bcast(eig_k, wfd%master, wfd%comm, ierr) 4643 call xmpi_bcast(allcg_k, wfd%master, wfd%comm, ierr) 4644 4645 bcount = 0 4646 do band=bmin,bmax 4647 if (my_readmask(band, ik_ibz, spin)) then 4648 bcount = bcount + 1 4649 cg_bpad = npw_disk * wfd%nspinor * (bcount-1) 4650 allcg_bpad = npw_disk * wfd%nspinor * (band - bmin) 4651 call cg_zcopy(npw_disk * wfd%nspinor, allcg_k(:, allcg_bpad+1), cg_k(:, cg_bpad+1)) 4652 end if 4653 end do 4654 ABI_FREE(allcg_k) 4655 end if 4656 4657 if (Wfd%my_rank == Wfd%master .and. wfd%prtvol > 0) then 4658 if (Wfd%nsppol==2) then 4659 write(std_out,'(i3,a,10f7.2/50(10x,10f7.2/))') ik_ibz,tag_spin(spin),(eig_k(ib)*Ha_eV,ib=1,nband_wfd) 4660 else 4661 write(std_out,'(i3,7x,10f7.2/50(10x,10f7.2/))')ik_ibz,(eig_k(ib)*Ha_eV,ib=1,nband_wfd) 4662 end if 4663 end if 4664 4665 ! Table with the correspondence btw the k-centered sphere of the WFK file 4666 ! and the one used in Wfd (possibly smaller due to ecutwfn). 4667 ! TODO: Here I should treat the case in which istwfk in wfd differs from the one on disk. 4668 ABI_MALLOC(gf2wfd, (npw_disk)) 4669 if (any(my_readmask(:,ik_ibz,spin))) then 4670 call kg_map(wfd%npwarr(ik_ibz), wfd%kdata(ik_ibz)%kg_k, npw_disk, kg_k, gf2wfd, nmiss) 4671 end if 4672 !if (nmiss/=0) then 4673 ! write(msg,'(a,2(1x,i0),a,i0)')" For (k,s) ",ik_ibz,spin," the number of missing G is ",nmiss 4674 ! ABI_WARNING(msg) 4675 !end if 4676 4677 if (change_gsphere .and. any(my_readmask(:,ik_ibz,spin))) then 4678 ! Prepare call to ctgk_change_sphere 4679 ! FFT box must enclose the two spheres (wfd(k), wfk(k)) 4680 gmax_wfd = maxval(abs(wfd%kdata(ik_ibz)%kg_k), dim=2) 4681 gmax_disk = maxval(abs(kg_k), dim=2) 4682 do ii=1,3 4683 gmax(ii) = max(gmax_wfd(ii), gmax_disk(ii)) 4684 end do 4685 gmax = 2*gmax + 1 4686 call ngfft_seq(work_ngfft, gmax) 4687 ABI_MALLOC(work, (2, work_ngfft(4),work_ngfft(5),work_ngfft(6))) 4688 ABI_MALLOC(out_cg, (2, wfd%npwarr(ik_ibz) * wfd%nspinor)) 4689 end if 4690 4691 ! Conversion of the basis set. 4692 bcount = 0 4693 do band=1,Wfd%nband(ik_ibz,spin) 4694 if (my_readmask(band, ik_ibz, spin)) then 4695 4696 ib = wfd%bks2wfd(1, band, ik_ibz, spin) 4697 ik = wfd%bks2wfd(2, band, ik_ibz, spin) 4698 is = wfd%bks2wfd(3, band, ik_ibz, spin) 4699 ABI_CHECK(all(wfd%bks2wfd(:, band, ik_ibz, spin) /= 0), "state in not allocated") 4700 4701 wave => wfd%s(is)%k(ik)%b(ib) 4702 wave%ug = czero 4703 4704 bcount = bcount + 1 4705 cg_bpad = npw_disk*Wfd%nspinor*(bcount-1) 4706 4707 if (change_gsphere) then 4708 ! Different istwfk storage. 4709 call cgtk_change_gsphere(wfd%nspinor, & 4710 npw_disk, istwfk_disk, kg_k, cg_k(:, cg_bpad+1:), & 4711 wfd%npwarr(ik_ibz), wfd%istwfk(ik_ibz), wfd%kdata(ik_ibz)%kg_k, out_cg, work_ngfft, work) 4712 4713 wave%ug(:) = CMPLX(out_cg(1, :), out_cg(2, :), kind=gwpc) 4714 !call wfd%push_ug(band, ik_ibz, spin, cryst, out_cg) 4715 else 4716 do spinor=1,Wfd%nspinor 4717 cg_spad=(spinor-1)*npw_disk 4718 gw_spad=(spinor-1)*Wfd%npwarr(ik_ibz) 4719 do ig=1,npw_disk 4720 icg = ig+cg_spad+cg_bpad 4721 igw = gf2wfd(ig)+gw_spad 4722 if (gf2wfd(ig) /= 0) then 4723 wave%ug(igw) = CMPLX(cg_k(1,icg),cg_k(2,icg), kind=gwpc) 4724 end if 4725 end do 4726 !call wfd%push_ug(band, ik_ibz, spin, cryst, out_cg) 4727 end do 4728 end if 4729 4730 wave%has_ug = WFD_STORED 4731 end if 4732 end do 4733 4734 ABI_FREE(eig_k) 4735 ABI_FREE(kg_k) 4736 ABI_FREE(cg_k) 4737 ABI_FREE(gf2wfd) 4738 ABI_SFREE(work) 4739 ABI_SFREE(out_cg) 4740 4741 if (ik_ibz <= 10 .or. mod(ik_ibz, 200) == 0) then 4742 write(msg,'(4x,4(a,i0),a)') "Reading kpt [", ik_ibz, "/", wfd%nkibz, "] spin [", spin, "/", wfd%nsppol, "]" 4743 call cwtime_report(msg, cpu_ks, wall_ks, gflops_ks) 4744 end if 4745 end do !ik_ibz 4746 end do !spin 4747 4748 else 4749 ABI_ERROR(sjoin("Wrong method: ", itoa(method))) 4750 end if 4751 4752 call wfk%close() 4753 call Hdr%free() 4754 4755 ABI_FREE(my_readmask) 4756 ABI_FREE(all_countks) 4757 4758 ! Update the kbs table storing the distribution of the ug and set the MPI communicators. 4759 select type (wfd) 4760 class is (wfdgw_t) 4761 call wfd%update_bkstab() 4762 end select 4763 4764 call cwtime_report(" WFK IO", cpu, wall, gflops, end_str=ch10) 4765 call timab(300, 2, tsec) 4766 4767 end subroutine wfd_read_wfk
m_wfd/wfd_reset_ur_cprj [ Functions ]
NAME
wfd_reset_ur_cprj
FUNCTION
Reinitialize the storage mode of the ur treated by this node.
SOURCE
1459 subroutine wfd_reset_ur_cprj(Wfd) 1460 1461 !Arguments ------------------------------------ 1462 class(wfd_t),intent(inout) :: Wfd 1463 1464 !Local variables ------------------------------ 1465 !scalars 1466 integer :: ib, ik, is 1467 !************************************************************************ 1468 1469 do is=1,size(wfd%s) 1470 do ik=1,size(wfd%s(is)%k) 1471 do ib=1,size(wfd%s(is)%k(ik)%b) 1472 if (wfd%s(is)%k(ik)%b(ib)%has_ur == WFD_STORED) wfd%s(is)%k(ik)%b(ib)%has_ur = WFD_ALLOCATED 1473 if (wfd%usepaw == 1) then 1474 if (wfd%s(is)%k(ik)%b(ib)%has_cprj == WFD_STORED) wfd%s(is)%k(ik)%b(ib)%has_cprj = WFD_ALLOCATED 1475 end if 1476 end do 1477 end do 1478 end do 1479 1480 end subroutine wfd_reset_ur_cprj
m_wfd/wfd_sym_ug_kg [ Functions ]
NAME
wfd_sym_ug_kg
FUNCTION
Use crystalline symmetries and time reversal to reconstruct wavefunctions at kk_bz from the IBZ image kk_ibz. Return periodic part in G-space as well as list of G-vectors belonging to the G-sphere centered on kk_bz
INPUTS
ecut: Cutoff energy for planewave basis set. kk_bz: k-point in the BZ for output wavefunctions and G-vectors. kk_ibz: Symmetrical image of kk_bz in the IBZ. bstart: Initial band nband: Number of bands to symmetrize. spin: Spin index mpw: Maximum number of planewaves used to dimension arrays. indkk: Symmetry map kk_bz -> kk_ibz as computed by listkk. cryst: Crystalline structure and symmetries work_ngfft: Define the size of the workspace array work work: Workspace array used to symmetrize wavefunctions
OUTPUT
istwf_kbz: Time-reversal flag associated to output wavefunctions npw_kbz: Number of G-vectors in kk_bz G-sphere kg_kbz: G-vectors in reduced coordinates. cgs_kbz: Periodic part of wavefunctions at kk_bz
SOURCE
4136 subroutine wfd_sym_ug_kg(self, ecut, kk_bz, kk_ibz, bstart, nband, spin, mpw, indkk, cryst, & 4137 work_ngfft, work, istwf_kbz, npw_kbz, kg_kbz, cgs_kbz) 4138 4139 !Arguments ------------------------------------ 4140 !scalars 4141 integer,intent(in) :: bstart, nband, spin, mpw 4142 type(crystal_t),intent(in) :: cryst 4143 class(wfd_t),intent(in) :: self 4144 integer,intent(out) :: istwf_kbz, npw_kbz 4145 real(dp),intent(in) :: ecut 4146 !arrays 4147 integer :: work_ngfft(18) 4148 integer,intent(in) :: indkk(6) 4149 integer,intent(out) :: kg_kbz(3, mpw) 4150 real(dp),intent(in) :: kk_bz(3), kk_ibz(3) 4151 real(dp),intent(out) :: cgs_kbz(2, mpw*self%nspinor, nband) 4152 real(dp),intent(out) :: work(2, work_ngfft(4), work_ngfft(5), work_ngfft(6)) 4153 4154 !Local variables ------------------------------ 4155 !scalars 4156 integer,parameter :: ndat1 = 1 4157 integer :: ik_ibz, isym_k, trev_k, ib, band, istwf_kirr, npw_kirr 4158 logical :: isirr_k 4159 !arrays 4160 integer :: g0_k(3) 4161 integer,allocatable :: gtmp(:,:) 4162 real(dp),allocatable :: cg_kirr(:,:) 4163 4164 !************************************************************************ 4165 4166 ! As reported by listkk via symrel 4167 ik_ibz = indkk(1); isym_k = indkk(2); trev_k = indkk(6); g0_k = indkk(3:5) 4168 isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0)) 4169 4170 ! Get npw_kbz, kg_kbz and symmetrize wavefunctions from IBZ (if needed). 4171 ! Be careful with time-reversal symmetry. 4172 if (isirr_k) then 4173 ! Copy u_k(G) 4174 istwf_kbz = self%istwfk(ik_ibz); npw_kbz = self%npwarr(ik_ibz) 4175 ABI_CHECK(mpw >= npw_kbz, "mpw < npw_kbz") 4176 kg_kbz(:,1:npw_kbz) = self%kdata(ik_ibz)%kg_k 4177 do ib=1,nband 4178 band = ib + bstart - 1 4179 call self%copy_cg(band, ik_ibz, spin, cgs_kbz(1,1,ib)) 4180 end do 4181 else 4182 ! Reconstruct u_k(G) from the IBZ image. 4183 istwf_kbz = 1 4184 call get_kg(kk_bz, istwf_kbz, ecut, cryst%gmet, npw_kbz, gtmp) 4185 ABI_CHECK(mpw >= npw_kbz, "mpw < npw_kbz") 4186 kg_kbz(:,1:npw_kbz) = gtmp(:,:npw_kbz) 4187 ABI_FREE(gtmp) 4188 4189 ! Use cg_kirr as workspace array, results stored in cgs_kbz. 4190 istwf_kirr = self%istwfk(ik_ibz); npw_kirr = self%npwarr(ik_ibz) 4191 ABI_MALLOC(cg_kirr, (2, npw_kirr*self%nspinor)) 4192 do ib=1,nband 4193 band = ib + bstart - 1 4194 call self%copy_cg(band, ik_ibz, spin, cg_kirr) 4195 call cgtk_rotate(cryst, kk_ibz, isym_k, trev_k, g0_k, self%nspinor, ndat1, & 4196 npw_kirr, self%kdata(ik_ibz)%kg_k, & 4197 npw_kbz, kg_kbz, istwf_kirr, istwf_kbz, cg_kirr, cgs_kbz(:,:,ib), work_ngfft, work) 4198 end do 4199 ABI_FREE(cg_kirr) 4200 end if 4201 4202 end subroutine wfd_sym_ug_kg
m_wfd/wfd_sym_ur [ Functions ]
NAME
wfd_sym_ur
FUNCTION
Symmetrize a wave function in real space
INPUTS
Wfd<wfd_t>=the wavefunction descriptor. Cryst<crystal_t>=Structure describing the crystal structure and its symmetries. Kmesh<kmesh_t>=Structure describing the BZ sampling band=Band index. ik_bz=Index of the k-point in the BZ. spin=Spin index [trans] = "N" if only the symmetried wavefunction is needed, "C" if the complex conjugate is required. Default is "N" [with_umklp] = Optional flag. If .True. (Default) the umklapp G0 vector in the relation kbz = Sk + G0 is taken into account when constructing u_kbz.
NOTES
This method is deprecated. See wfd_sym_ug_kg for symmetrization in G-space
OUTPUT
ur_kbz(Wfd%nfft*Wfd%nspinor)=The symmetrized wavefunction in real space. [ur_kibz(Wfd%nfft*Wfd%nspinor)]= Optional output: u(r) in the IBZ.
SOURCE
3959 subroutine wfd_sym_ur(Wfd,Cryst,Kmesh,band,ik_bz,spin,ur_kbz,trans,with_umklp,ur_kibz) 3960 3961 !Arguments ------------------------------------ 3962 !scalars 3963 integer,intent(in) :: band,ik_bz,spin 3964 character(len=*),optional,intent(in) :: trans 3965 logical,optional,intent(in) :: with_umklp 3966 type(crystal_t),intent(in) :: Cryst 3967 type(kmesh_t),intent(in) :: Kmesh 3968 class(wfd_t),intent(inout) :: Wfd 3969 !arrays 3970 complex(gwpc),intent(out) :: ur_kbz(Wfd%nfft*Wfd%nspinor) 3971 complex(gwpc),optional,intent(out) :: ur_kibz(Wfd%nfft*Wfd%nspinor) 3972 3973 !Local variables ------------------------------ 3974 !scalars 3975 integer :: ik_ibz,isym_k,itim_k,nr,ispinor,spad,ir,ir2 3976 integer :: fft_idx,ix,iy,iz,nx,ny,nz,irot 3977 real(dp) :: gdotr 3978 complex(dpc) :: ph_mkt,u2b,u2a 3979 complex(gwpc) :: gwpc_ph_mkt 3980 logical :: isirred,my_with_umklp 3981 character(len=1) :: my_trans 3982 !character(len=500) :: msg 3983 !arrays 3984 integer :: umklp(3) 3985 real(dp) :: kbz(3),spinrot_k(4) 3986 complex(dpc) :: spinrot_mat(2,2) 3987 complex(gwpc),allocatable :: ur(:) 3988 3989 !************************************************************************ 3990 3991 my_trans = "N"; if (present(trans)) my_trans = toupper(trans(1:1)) 3992 my_with_umklp = .TRUE.; if (present(with_umklp)) my_with_umklp = with_umklp 3993 3994 ! k_bz = S k - G0 ==> u_{k_bz} = e^{iG0.r} u_{Sk} 3995 ! k_bz = -S k - G0 ==> u_{k_bz} = e^{iG0.r} u_{Sk}^* 3996 3997 ! u(r,b,kbz)=e^{-2i\pi kibz.(R^{-1}t} u (R{^-1}(r-t),b,kibz) 3998 ! =e^{+2i\pi kibz.(R^{-1}t} u*({R^-1}(r-t),b,kibz) for time-reversal 3999 ! 4000 ! Get ik_ibz, non-symmorphic phase, ph_mkt, and symmetries from ik_bz. 4001 call Kmesh%get_BZ_item(ik_bz,kbz,ik_ibz,isym_k,itim_k,ph_mkt,umklp,isirred) 4002 gwpc_ph_mkt = ph_mkt 4003 4004 if (isirred) then 4005 ! Avoid symmetrization if this point is irreducible. 4006 call wfd%get_ur(band,ik_ibz,spin,ur_kbz) 4007 if (present(ur_kibz)) call xcopy(Wfd%nfft*Wfd%nspinor,ur_kbz,1,ur_kibz,1) 4008 if (my_trans=="C") ur_kbz = GWPC_CONJG(ur_kbz) 4009 RETURN 4010 end if 4011 4012 ! Reconstruct ur in the BZ from the corresponding wavefunction in IBZ. 4013 ABI_MALLOC(ur, (Wfd%nfft*Wfd%nspinor)) 4014 4015 call wfd%get_ur(band,ik_ibz,spin,ur) 4016 if (present(ur_kibz)) call xcopy(Wfd%nfft*Wfd%nspinor,ur,1,ur_kibz,1) 4017 4018 ! Wfd%irottb(:,isym_k) is the table for rotated FFT points 4019 SELECT CASE (Wfd%nspinor) 4020 4021 CASE (1) 4022 ! Rotation in real space 4023 do ir=1,Wfd%nfft 4024 irot = Wfd%irottb(ir,isym_k) 4025 ur_kbz(ir) = ur(irot) * gwpc_ph_mkt 4026 end do 4027 4028 ! Apply time-reversal symmetry if needed. 4029 if (itim_k==2) ur_kbz = GWPC_CONJG(ur_kbz) 4030 4031 ! Take into account a possible umklapp. 4032 if (ANY(umklp/=0).and. my_with_umklp) then 4033 ! Compute ur_kbz = ur_kbz*eig0r 4034 nx = Wfd%ngfft(1); ny = Wfd%ngfft(2); nz = Wfd%ngfft(3) 4035 fft_idx=0 4036 do iz=0,nz-1 4037 do iy=0,ny-1 4038 do ix=0,nx-1 4039 gdotr= two_pi*( umklp(1)*(ix/DBLE(nx)) & 4040 +umklp(2)*(iy/DBLE(ny)) & 4041 +umklp(3)*(iz/DBLE(nz)) ) 4042 fft_idx = fft_idx+1 4043 ur_kbz(fft_idx) = ur_kbz(fft_idx) * DCMPLX(DCOS(gdotr),DSIN(gdotr)) 4044 end do 4045 end do 4046 end do 4047 end if 4048 4049 if (my_trans=="C") ur_kbz = GWPC_CONJG(ur_kbz) 4050 4051 CASE (2) 4052 ABI_ERROR("Implementation has to be tested") 4053 4054 nr = Wfd%nfft 4055 spinrot_k = Cryst%spinrot(:,isym_k) 4056 ! 4057 ! ==== Apply Time-reversal if required ==== 4058 ! \psi_{-k}^1 = (\psi_k^2)^* 4059 ! \psi_{-k}^2 = -(\psi_k^1)^* 4060 if (itim_k==1) then 4061 ur_kbz = ur 4062 else if (itim_k==2) then 4063 ur_kbz(1:nr) = GWPC_CONJG(ur(nr+1:2*nr)) 4064 ur_kbz(nr+1:2*nr)=-GWPC_CONJG(ur(1:nr)) 4065 else 4066 ABI_ERROR('Wrong i2 in spinor') 4067 end if 4068 ! 4069 ! Rotate wavefunctions in real space. 4070 do ispinor=1,Wfd%nspinor 4071 spad=(ispinor-1)*nr 4072 do ir=1,nr 4073 ir2 = Wfd%irottb(ir,isym_k) 4074 ur_kbz(ir+spad) = ur_kbz(ir2+spad) * gwpc_ph_mkt 4075 end do 4076 end do 4077 ! 4078 ! Rotation in spinor space. 4079 spinrot_mat(1,1)= spinrot_k(1) + j_dpc*spinrot_k(4) 4080 spinrot_mat(1,2)= spinrot_k(3) + j_dpc*spinrot_k(2) 4081 spinrot_mat(2,1)=-spinrot_k(3) + j_dpc*spinrot_k(2) 4082 spinrot_mat(2,2)= spinrot_k(1) - j_dpc*spinrot_k(4) 4083 4084 do ir=1,nr 4085 u2a=ur_kbz(ir) 4086 u2b=ur_kbz(ir+nr) 4087 ur_kbz(ir) =spinrot_mat(1,1)*u2a+spinrot_mat(1,2)*u2b 4088 ur_kbz(ir+nr)=spinrot_mat(2,1)*u2a+spinrot_mat(2,2)*u2b 4089 end do 4090 4091 if (ANY(umklp /=0)) then 4092 !ur_kbz(1:Wfd%nfft) = ur_kbz(1:Wfd%nfft) *eig0r 4093 !ur_kbz(Wfd%nfft+1:) = ur_kbz(Wfd%nfft+1:)*eig0r 4094 end if 4095 4096 CASE DEFAULT 4097 ABI_ERROR(sjoin("Wrong value for nspinor: ", itoa(Wfd%nspinor))) 4098 END SELECT 4099 4100 ABI_FREE(ur) 4101 4102 end subroutine wfd_sym_ur
m_wfd/wfd_t [ Types ]
NAME
wfd_t
FUNCTION
Container gathering information on the set of wavefunctions treated by this node This is a base class without the bks_tab that is used in the GW/BSE code.
SOURCE
270 type, public :: wfd_t 271 272 integer :: debug_level = 0 ! Internal flag defining the debug level. 273 integer :: lmnmax 274 integer :: mband ! MAX(nband) 275 integer :: mgfft ! Maximum size of 1D FFTs i.e. MAXVAL(ngfft(1:3)), used to dimension some arrays. 276 integer :: natom 277 integer :: nfft ! Number of FFT points treated by this processor 278 integer :: nfftot ! Total number of points in the FFT grid 279 integer :: nkibz ! Number of irreducible k-points 280 integer :: nspden ! Number of independent spin-density components 281 integer :: nspinor ! Number of spinor components 282 integer :: nsppol ! Number of independent spin polarizations 283 integer :: ntypat 284 integer :: paral_kgb ! Option for kgb parallelism 285 integer :: usepaw ! 1 if PAW is used, 0 otherwise. 286 integer :: prtvol ! Verbosity level. 287 integer :: pawprtvol ! Verbosity level for PAW. 288 integer :: usewvl ! 1 if BigDFT is used, 0 otherwise. 289 integer :: comm ! The MPI communicator for this pool of processors. 290 integer :: master ! The rank of master node in comm. 291 integer :: my_rank ! The rank of my processor inside the MPI communicator comm. 292 integer :: nproc ! The number of processors in MPI comm. 293 294 integer :: my_nspins ! Number of spins treated by this MPI proc 295 296 integer,allocatable :: my_nkspin(:) 297 ! (%nsppol)) 298 ! Number of k-points treated by this MPI proc. 299 300 logical :: rfft_is_symok ! .TRUE. if the real space FFT mesh is compatible with the rotational 301 ! part of the space group. 302 303 logical :: use_fnl_dir0der0 ! .TRUE. if the the Wfd must store fnl dir0der0. 304 305 real(dp) :: dilatmx 306 307 real(dp) :: ecut 308 ! Cutoff for plane wave basis set. 309 310 real(dp) :: ecutsm 311 ! smearing energy for plane wave kinetic energy (Ha) 312 ! Cutoff for plane wave basis set. 313 314 integer :: ngfft(18) 315 ! Information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft 316 317 integer :: nloalg(3) 318 ! Governs the choice of the algorithm for nonlocal operator. See doc. 319 320 integer,allocatable :: irottb(:,:) 321 ! (nfftot, nsym) 322 ! Index of $R^{-1}(r-\tau)$ in the FFT box. 323 324 integer,allocatable :: istwfk(:) 325 ! (nkibz) 326 ! Storage mode for this k-point. 327 328 integer,allocatable :: nband(:,:) 329 ! (nkibz,nsppol) 330 ! Number of bands at each k-point and spin. 331 332 integer,allocatable :: indlmn(:,:,:) 333 ! (6, lmnmax, ntypat) 334 ! array giving l,m,n,lm,ln,spin for i=ln (if useylm=0) 335 ! or i=lmn (if useylm=1) 336 337 integer,allocatable :: nlmn_atm(:) 338 ! (natom) 339 ! Number of (n,l,m) channels for each atom. Only for PAW 340 341 integer,allocatable :: nlmn_sort(:) 342 ! (natom) 343 ! Number of (n,l,m) channels for each atom (sorted by atom type). Only for PAW 344 345 integer,allocatable :: nlmn_type(:) 346 ! (ntypat) 347 ! Number of (n,l,m) channels for each type of atom. Only for PAW. 348 349 integer,allocatable :: npwarr(:) 350 ! (nkibz) 351 ! Number of plane waves for this k-point. 352 353 integer, allocatable :: bks2wfd(:,:,:,:) 354 ! (3, mband, nkibz, nsppol) 355 ! Maps global (band, ik_ibz, spin) to index in the wave store. 356 ! Set to 0 if the (b, k, s) state is not in the store. 357 358 real(dp),allocatable :: kibz(:,:) 359 ! (3, nkibz) 360 ! Reduced coordinates of the k-points in the IBZ. 361 362 real(dp),allocatable :: ph1d(:,:) 363 ! (2,3*(2*mgfft+1)*natom) 364 ! 1-dim structure factor phase information. 365 366 logical,private, allocatable :: keep_ur(:,:,:) 367 ! TODO: To be removed 368 ! keep(mband,nkibz,nsppol) 369 ! Storage strategy: keep or not keep calculated u(r) in memory. 370 371 type(kdata_t),allocatable :: kdata(:) 372 ! (nkibz) 373 ! datatype storing k-dependent quantities. 374 375 type(spin_store_t),allocatable :: s(:) 376 ! (my_nsppol) 377 ! wfd%s(is)%k(ik)%b(ib) 378 379 type(MPI_type) :: MPI_enreg 380 ! The MPI_type structured datatype gather different information about the MPI parallelisation: 381 ! number of processors, the index of my processor, the different groups of processors, etc ... 382 383 !type(pseudopotential_type), pointer :: psps 384 !type(pawtab_type), pointer :: pawtab(:) 385 386 contains 387 388 procedure :: free => wfd_free 389 ! Free memory. 390 391 procedure :: norm2 => wfd_norm2 392 ! Compute <u(g)|u(g)> for the same k-point and spin. 393 394 procedure :: xdotc => wfd_xdotc 395 ! Compute <u_{b1ks}|u_{b2ks}> in G-space 396 397 procedure :: reset_ur_cprj => wfd_reset_ur_cprj 398 ! Reinitialize memory storage of u(r) and <p_i|psi> 399 400 procedure :: get_many_ur => wfd_get_many_ur 401 ! Get many wavefunctions in real space from its (bands(:),k,s) indices. 402 403 procedure :: copy_cg => wfd_copy_cg 404 ! Return a copy of u(g) in a real(2,npw_k)) array (Abinit convention) 405 406 procedure :: get_ur => wfd_get_ur 407 ! Get one wavefunction in real space from its (b,k,s) indices. 408 409 procedure :: get_cprj => wfd_get_cprj 410 ! Get one PAW projection <Proj_i|Cnk> with all NL projectors from its (b,k,s) indices. 411 412 procedure :: change_ngfft => wfd_change_ngfft 413 ! Reinitialize internal FFT tables. 414 415 procedure :: print => wfd_print 416 ! Printout of basic info. 417 418 procedure :: ug2cprj => wfd_ug2cprj 419 ! Get PAW cprj from its (b,k,s) indices. 420 421 procedure :: wave_free => wfd_wave_free 422 ! Free internal buffers used to store the wavefunctions. 423 424 procedure :: get_wave_ptr => wfd_get_wave_ptr 425 ! Return pointer to wave_t from its (b,k,s) indices 426 427 procedure :: push_ug => wfd_push_ug 428 ! Modify the value of u(g)_ks stored in the object. 429 430 procedure :: extract_cgblock => wfd_extract_cgblock 431 ! Extract a block of wavefunctions for a given spin and k-points (uses the cg storage mode) 432 433 procedure :: ihave_ug => wfd_ihave_ug 434 ! True if the node has this ug with the specified status. 435 436 procedure :: mybands => wfd_mybands 437 ! Returns the list of band indices of the u(g) owned by this node at given (k,s). 438 439 procedure :: test_ortho => wfd_test_ortho 440 ! Test the orthonormalization of the wavefunctions. 441 442 procedure :: sym_ur => wfd_sym_ur 443 ! Symmetrize a wave function in real space 444 ! This routine is deprecated, see wfd_sym_ug_kg for algo in G-space. 445 446 procedure :: sym_ug_kg => wfd_sym_ug_kg 447 ! Symmetrize a wave function in G-space 448 ! Used in phgamma only, see sigmaph for a more efficient version. 449 450 procedure :: paw_get_aeur => wfd_paw_get_aeur 451 ! Compute the AE PAW wavefunction in real space. 452 453 procedure :: read_wfk => wfd_read_wfk 454 ! Read u(g) from the WFK file completing the initialization of the object. 455 456 procedure :: dump_errinfo => wfd_dump_errinfo 457 458 end type wfd_t 459 460 public :: wfd_init ! Main creation method. 461 462 !public :: wfd_get_socpert 463 public :: test_charge
m_wfd/wfd_test_ortho [ Functions ]
NAME
wfd_test_ortho
FUNCTION
Test the orthonormalization of the wavefunctions stored in Wfd.
INPUTS
Wfd<wfd_t>=wavefunction descriptor. Cryst<crystal_t>=Object defining the unit cell and its symmetries. Pawtab(ntypat*usepaw)<type(pawtab_type)>=PAW tabulated starting data.
OUTPUT
Only writing.
SOURCE
3780 subroutine wfd_test_ortho(Wfd,Cryst,Pawtab,unit,mode_paral) 3781 3782 !Arguments ------------------------------------ 3783 !scalars 3784 integer,intent(in),optional :: unit 3785 character(len=4),optional,intent(in) :: mode_paral 3786 type(crystal_t),intent(in) :: Cryst 3787 class(wfd_t),target,intent(inout) :: Wfd 3788 !array 3789 type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw) 3790 3791 !Local variables ------------------------------ 3792 !scalars 3793 integer :: ik_ibz,spin,band,band1,band2,ib,ib1,ib2,ierr,how_manyb,my_unt,npw_k,istwf_k 3794 real(dp) :: glob_cinf,my_cinf,glob_csup,my_csup,glob_einf,min_norm2,glob_esup,max_norm2 3795 complex(dpc) :: cdum 3796 logical :: bands_are_spread 3797 character(len=4) :: my_mode 3798 character(len=500) :: msg 3799 type(wave_t),pointer :: wave1, wave2 3800 !arrays 3801 integer :: my_bandlist(Wfd%mband) 3802 real(dp) :: pawovlp(2) 3803 complex(gwpc),ABI_CONTIGUOUS pointer :: ug1(:),ug2(:) 3804 !complex(gwpc) :: ur(Wfd%nfft*Wfd%nspinor) 3805 character(len=6) :: tag_spin(2) 3806 type(pawcprj_type),allocatable :: Cp1(:,:),Cp2(:,:) 3807 3808 !************************************************************************ 3809 3810 tag_spin(:)=(/' ',' '/); if (Wfd%nsppol==2) tag_spin(:)=(/' UP ',' DOWN '/) 3811 3812 my_unt =std_out; if (present(unit )) my_unt =unit 3813 my_mode ='COLL' ; if (present(mode_paral)) my_mode =mode_paral 3814 3815 if (Wfd%usepaw==1) then 3816 ABI_MALLOC(Cp1,(Wfd%natom,Wfd%nspinor)) 3817 call pawcprj_alloc(Cp1,0,Wfd%nlmn_atm) 3818 ABI_MALLOC(Cp2,(Wfd%natom,Wfd%nspinor)) 3819 call pawcprj_alloc(Cp2,0,Wfd%nlmn_atm) 3820 end if 3821 3822 bands_are_spread = .FALSE. 3823 3824 do spin=1,Wfd%nsppol 3825 min_norm2=greatest_real; max_norm2=-greatest_real 3826 my_cinf=greatest_real; my_csup=-greatest_real 3827 do ik_ibz=1,Wfd%nkibz 3828 istwf_k = Wfd%istwfk(ik_ibz) 3829 npw_k = Wfd%npwarr(ik_ibz) 3830 3831 ! Select my band indices. 3832 call wfd%mybands(ik_ibz,spin,how_manyb,my_bandlist, how="Stored") 3833 if (how_manyb/=Wfd%nband(ik_ibz,spin)) bands_are_spread = .TRUE. 3834 3835 ! 1) Normalization. 3836 do ib=1,how_manyb 3837 band = my_bandlist(ib) 3838 ABI_CHECK(wfd%get_wave_ptr(band, ik_ibz, spin, wave1, msg) == 0, msg) 3839 ug1 => wave1%ug 3840 cdum = xdotc(npw_k*Wfd%nspinor,ug1,1,ug1,1) 3841 if (istwf_k>1) then 3842 cdum=two*DBLE(cdum) 3843 if (istwf_k==2) cdum=cdum-CONJG(ug1(1))*ug1(1) 3844 end if 3845 if (Wfd%usepaw==1) then 3846 call wfd%get_cprj(band,ik_ibz,spin,Cryst,Cp1,sorted=.FALSE.) 3847 pawovlp = paw_overlap(Cp1,Cp1,Cryst%typat,Pawtab,spinor_comm=Wfd%MPI_enreg%comm_spinor) 3848 cdum = cdum + CMPLX(pawovlp(1),pawovlp(2), kind=dpc) 3849 end if 3850 !write(std_out,*)"ik_ibz, band, spin, cdum: ",ik_ibz,band,spin,cdum 3851 if (REAL(cdum)<min_norm2) min_norm2=REAL(cdum) 3852 if (REAL(cdum)>max_norm2) max_norm2=REAL(cdum) 3853 end do 3854 3855 ! TODO should use the communicator for this spin 3856 call xmpi_min(min_norm2,glob_einf,Wfd%comm,ierr) 3857 call xmpi_max(max_norm2,glob_esup,Wfd%comm,ierr) 3858 3859 ! 2) Orthogonality of wavefunctions. 3860 do ib1=1,how_manyb 3861 band1 = my_bandlist(ib1) 3862 ABI_CHECK(wfd%get_wave_ptr(band1, ik_ibz, spin, wave1, msg) == 0, msg) 3863 ug1 => wave1%ug 3864 if (Wfd%usepaw==1) call wfd%get_cprj(band1,ik_ibz,spin,Cryst,Cp1,sorted=.FALSE.) 3865 3866 do ib2=ib1+1,how_manyb 3867 band2 = my_bandlist(ib2) 3868 ABI_CHECK(wfd%get_wave_ptr(band2, ik_ibz, spin, wave2, msg) == 0, msg) 3869 ug2 => wave2%ug 3870 if (Wfd%usepaw==1) call wfd%get_cprj(band2,ik_ibz,spin,Cryst,Cp2,sorted=.FALSE.) 3871 3872 cdum = xdotc(npw_k*Wfd%nspinor,ug1,1,ug2,1) 3873 if (istwf_k>1) then 3874 cdum=two*DBLE(cdum) 3875 if (istwf_k==2) cdum=cdum-CONJG(ug1(1))*ug2(1) 3876 end if 3877 if (Wfd%usepaw==1) then 3878 pawovlp = paw_overlap(Cp1,Cp2,Cryst%typat,Pawtab,spinor_comm=Wfd%MPI_enreg%comm_spinor) 3879 cdum = cdum + CMPLX(pawovlp(1),pawovlp(2), kind=dpc) 3880 end if 3881 3882 if (ABS(cdum)<my_cinf) my_cinf=ABS(cdum) 3883 if (ABS(cdum)>my_csup) my_csup=ABS(cdum) 3884 !if (ABS(cdum) > 0.1) write(std_out,*)" ib1,ib2,ABS_dotprod: ",ib1,ib2,ABS(cdum) 3885 end do !ib2 3886 end do !ib 3887 3888 ! TODO should use the communicator for this spin 3889 call xmpi_min(my_cinf,glob_cinf,Wfd%comm,ierr) 3890 call xmpi_max(my_csup,glob_csup,Wfd%comm,ierr) 3891 end do ! ik_ibz 3892 3893 ! Output results for this spin 3894 write(msg,'(2a)')ch10,' test on the normalization of the wavefunctions' 3895 if (Wfd%nsppol==2) write(msg,'(3a)')ch10,' test on the normalization of the wavefunctions with spin ',tag_spin(spin) 3896 call wrtout(my_unt,msg,mode_paral) 3897 write(msg,'(a,f9.6,a,a,f9.6)')& 3898 ' min sum_G |a(n,k,G)| = ',glob_einf,ch10,& 3899 ' max sum_G |a(n,k,G)| = ',glob_esup 3900 call wrtout(my_unt,msg,mode_paral) 3901 3902 write(msg,'(a)')' test on the orthogonalization of the wavefunctions (NB: this is not invariant for degenerate states)' 3903 if (Wfd%nsppol==2) write(msg,'(2a)')' test on the orthogonalization of the wavefunctions with spin ',tag_spin(spin) 3904 call wrtout(my_unt,msg,mode_paral) 3905 write(msg,'(a,f9.6,a,a,f9.6,a)')& 3906 '- min sum_G a(n,k,G)a(n",k,G) = ',glob_cinf,ch10,& 3907 '- max sum_G a(n,k,G)a(n",k,G) = ',glob_csup,ch10 3908 call wrtout(my_unt,msg,mode_paral) 3909 3910 end do ! spin 3911 3912 if (bands_are_spread) then 3913 write(msg,'(3a)')& 3914 'Note that the test on the orthogonalization is not complete ',ch10,& 3915 'since bands are spread among different processors' 3916 call wrtout(my_unt,msg,mode_paral) 3917 end if 3918 3919 if (Wfd%usepaw==1) then 3920 call pawcprj_free(Cp1) 3921 ABI_FREE(Cp1) 3922 call pawcprj_free(Cp2) 3923 ABI_FREE(Cp2) 3924 end if 3925 3926 end subroutine wfd_test_ortho
m_wfd/wfd_ug2cprj [ Functions ]
NAME
wfd_ug2cprj
FUNCTION
Calculate the projected wave function <Proj_i|Cnk> with all NL projectors for a single k-point, band and spin.
INPUTS
Wfd<wfd_t>=Structure containing the wave functions for the GW. ik_ibz=Index of the required k-point spin=Required spin index. choice=chooses possible output: In addition to projected wave function: choice=1 => nothing else =2 => 1st gradients with respect to atomic position(s) =3 => 1st gradients with respect to strain(s) =23=> 1st gradients with respect to atm. pos. and strain(s) =4 => 2nd derivatives with respect to atomic pos. =24=> 1st and 2nd derivatives with respect to atomic pos. =5 => 1st gradients with respect to k wavevector =6 => 2nd derivatives with respect to strain and atm. pos. idir=direction of the derivative, i.e. dir. of - atom to be moved in the case choice=2 - strain component in the case choice=3 - k point direction in the case choice=5 Compatible only with choice=2,3,5; if idir=0, all derivatives are computed natom Cryst [sorted]=Logical flags defining if the output Cprj has to be sorted by atom type or not. By default, Cprj matrix elements are unsorted.
OUTPUT
cwaveprj
SOURCE
1818 subroutine wfd_ug2cprj(Wfd,band,ik_ibz,spin,choice,idir,natom,Cryst,cwaveprj,sorted) 1819 1820 !Arguments ------------------------------- 1821 !scalars 1822 integer,intent(in) :: choice,idir,natom,band,ik_ibz,spin 1823 logical,optional,intent(in) :: sorted 1824 class(wfd_t),target,intent(inout) :: Wfd 1825 type(crystal_t),intent(in) :: Cryst 1826 !arrays 1827 type(pawcprj_type),intent(inout) :: cwaveprj(natom,Wfd%nspinor) 1828 1829 !Local variables------------------------------- 1830 !scalars 1831 integer :: cpopt,istwf_k,npw_k,nkpg 1832 integer :: ia,iatm,dimffnl,itypat,iatom,isp 1833 character(len=500) :: msg 1834 type(wave_t),pointer :: wave 1835 logical :: want_sorted 1836 !arrays 1837 integer,ABI_CONTIGUOUS pointer :: kg_k(:,:) 1838 integer,allocatable :: dimcprj_srt(:) 1839 real(dp) :: kpoint(3) 1840 real(dp),ABI_CONTIGUOUS pointer :: phkxred(:,:) 1841 real(dp),allocatable :: cwavef(:,:), kpg(:,:) 1842 !real(dp),allocatable :: ph1d(2,3*(2*mgfft+1)*natom) 1843 real(dp),ABI_CONTIGUOUS pointer :: ph3d(:,:,:) ! ph3d(2,npw_k,matblk) 1844 real(dp),ABI_CONTIGUOUS pointer :: ffnl(:,:,:,:) ! ffnl(npw_k,dimffnl,lmnmax,ntypat) 1845 type(pawcprj_type),allocatable :: Cprj_srt(:,:) 1846 1847 ! ********************************************************************* 1848 1849 ! Different form factors have to be calculated and stored in Kdata. 1850 ABI_CHECK_IEQ(choice, 1, "choice/=1 not coded") 1851 1852 dimffnl = 1 1853 npw_k = Wfd%npwarr(ik_ibz) 1854 istwf_k = Wfd%istwfk(ik_ibz) 1855 kpoint = Wfd%kibz(:,ik_ibz) 1856 1857 kg_k => Wfd%Kdata(ik_ibz)%kg_k 1858 ph3d => Wfd%Kdata(ik_ibz)%ph3d 1859 ffnl => Wfd%Kdata(ik_ibz)%fnl_dir0der0 1860 phkxred => Wfd%Kdata(ik_ibz)%phkxred 1861 1862 ! Compute (k+G) vectors 1863 nkpg=0 1864 !% if (choice==3.or.choice==2.or.choice==23) nkpg=3*Wfd%nloalg(3) 1865 !% if (choice==4.or.choice==24) nkpg=9*Wfd%nloalg(3) 1866 ABI_MALLOC(kpg,(npw_k,nkpg)) 1867 if (nkpg>0) call mkkpg(kg_k,kpg,kpoint,nkpg,npw_k) 1868 1869 ! Copy wavefunction in G-space 1870 ABI_MALLOC(cwavef, (2,npw_k*Wfd%nspinor)) 1871 ABI_CHECK(wfd%get_wave_ptr(band, ik_ibz, spin, wave, msg) == 0, msg) 1872 cwavef(1,:) = DBLE (wave%ug) 1873 cwavef(2,:) = AIMAG(wave%ug) 1874 1875 cpopt = 0 ! Nothing is already calculated. 1876 1877 want_sorted=.FALSE.; if (present(sorted)) want_sorted=sorted 1878 1879 if (want_sorted) then 1880 ! Output cprj are sorted. 1881 call getcprj(choice,cpopt,cwavef,cwaveprj,ffnl,& 1882 idir,Wfd%indlmn,istwf_k,kg_k,kpg,kpoint,Wfd%lmnmax,Wfd%mgfft,Wfd%MPI_enreg,& 1883 Cryst%natom,Cryst%nattyp,Wfd%ngfft,Wfd%nloalg,npw_k,Wfd%nspinor,Cryst%ntypat,& 1884 phkxred,Wfd%ph1d,ph3d,Cryst%ucvol,1) 1885 1886 else 1887 ! Output cprj are unsorted. 1888 ABI_MALLOC(dimcprj_srt,(Cryst%natom)) 1889 ia=0 1890 do itypat=1,Cryst%ntypat 1891 dimcprj_srt(ia+1:ia+Cryst%nattyp(itypat))=Wfd%nlmn_type(itypat) 1892 ia=ia+Cryst%nattyp(itypat) 1893 end do 1894 1895 ABI_MALLOC(Cprj_srt,(natom,Wfd%nspinor)) 1896 call pawcprj_alloc(Cprj_srt,0,dimcprj_srt) 1897 ABI_FREE(dimcprj_srt) 1898 1899 ! Calculate sorted cprj. 1900 call getcprj(choice,cpopt,cwavef,Cprj_srt,ffnl,& 1901 idir,Wfd%indlmn,istwf_k,kg_k,kpg,kpoint,Wfd%lmnmax,Wfd%mgfft,Wfd%MPI_enreg,& 1902 Cryst%natom,Cryst%nattyp,Wfd%ngfft,Wfd%nloalg,npw_k,Wfd%nspinor,Cryst%ntypat,& 1903 phkxred,Wfd%ph1d,ph3d,Cryst%ucvol,1) 1904 1905 ! Reorder cprj (sorted --> unsorted) 1906 do iatom=1,Cryst%natom 1907 iatm=Cryst%atindx(iatom) 1908 do isp=1,Wfd%nspinor 1909 cwaveprj(iatom,isp)%cp=Cprj_srt(iatm,isp)%cp 1910 end do 1911 end do 1912 1913 call pawcprj_free(Cprj_srt) 1914 ABI_FREE(Cprj_srt) 1915 end if 1916 1917 ABI_FREE(cwavef) 1918 ABI_FREE(kpg) 1919 1920 end subroutine wfd_ug2cprj
m_wfd/wfd_wave_free [ Functions ]
NAME
wfd_wave_free
FUNCTION
Collection procedure that frees the set of waves specified by mask.
INPUTS
mask(Wfd%mband,Wfd%nkibz,Wfd%nsppol)=.TRUE. if the memory allocated for this state has to be freed [what]=String specifying which array have to be deallocated. Possible values (no case-sensitive). "All"= To free both ug and ur and PAW Cprj, if any. Default "G" = Only ug "R" = Only ur. "C" = Only PAW Cprj.
SIDE EFFECTS
Wfd<wfd_t>=See above.
SOURCE
2671 subroutine wfd_wave_free(Wfd, what, bks_mask) 2672 2673 !Arguments ------------------------------------ 2674 !scalars 2675 class(wfd_t),target,intent(inout) :: Wfd 2676 character(len=*),optional,intent(in) :: what 2677 !arrays 2678 logical,optional,intent(in) :: bks_mask(Wfd%mband,Wfd%nkibz,Wfd%nsppol) 2679 2680 !Local variables ------------------------------ 2681 !scalars 2682 integer :: ik_ibz, spin, band, ib, ik, is 2683 logical :: do_free 2684 type(wave_t),pointer :: wave 2685 !character(len=500) :: msg 2686 character(len=10) :: my_what 2687 !************************************************************************ 2688 2689 my_what="ALL"; if (present(what)) my_what=toupper(what) 2690 2691 do spin=1,Wfd%nsppol 2692 do ik_ibz=1,Wfd%nkibz 2693 do band=1,Wfd%nband(ik_ibz,spin) 2694 do_free=.TRUE.; if (present(bks_mask)) do_free=bks_mask(band,ik_ibz,spin) 2695 if (do_free) then 2696 ib = wfd%bks2wfd(1, band, ik_ibz, spin) 2697 ik = wfd%bks2wfd(2, band, ik_ibz, spin) 2698 is = wfd%bks2wfd(3, band, ik_ibz, spin) 2699 if (ib /= 0) then 2700 wave => wfd%s(is)%k(ik)%b(ib) 2701 call wave%free(what=my_what) 2702 end if 2703 select type (wfd) 2704 class is (wfdgw_t) 2705 ! Update the associated flags if we release the G-space. 2706 if ( firstchar(my_what, ["A", "G"])) Wfd%bks_tab(band, ik_ibz, spin, Wfd%my_rank) = WFD_NOWAVE 2707 end select 2708 end if 2709 end do 2710 end do 2711 end do 2712 2713 end subroutine wfd_wave_free
m_wfd/wfd_xdotc [ Functions ]
NAME
wfd_xdotc
FUNCTION
Compute <u_{b1ks}|u_{b2ks}> in G-space
INPUTS
Wfd<wfd_t>=the wavefunction descriptor. Cryst<crystal_t>=Structure describing the crystal structure and its symmetries. Pawtab(ntypat*usepaw)<type(pawtab_type)>=PAW tabulated starting data. band1, band2=Band indices. ik_bz=Index of the k-point in the BZ. spin=Spin index
SOURCE
1377 function wfd_xdotc(Wfd,Cryst,Pawtab,band1,band2,ik_ibz,spin) 1378 1379 !Arguments ------------------------------------ 1380 !scalars 1381 integer,intent(in) :: band1,band2,ik_ibz,spin 1382 complex(gwpc) :: wfd_xdotc 1383 class(wfd_t),target,intent(inout) :: Wfd 1384 type(crystal_t),intent(in) :: Cryst 1385 !arrays 1386 type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw) 1387 1388 !Local variables ------------------------------ 1389 !scalars 1390 integer :: npw_k,istwf_k 1391 type(wave_t),pointer :: wave1, wave2 1392 character(len=500) :: msg 1393 !arrays 1394 real(dp) :: pawovlp(2) 1395 complex(gwpc),ABI_CONTIGUOUS pointer :: ug1(:),ug2(:) 1396 type(pawcprj_type),allocatable :: Cp1(:,:),Cp2(:,:) 1397 1398 !************************************************************************ 1399 1400 ! Planewave part. 1401 npw_k = Wfd%npwarr(ik_ibz) 1402 istwf_k = Wfd%istwfk(ik_ibz) 1403 1404 ABI_CHECK(wfd%get_wave_ptr(band1, ik_ibz, spin, wave1, msg) == 0, msg) 1405 ABI_CHECK(wfd%get_wave_ptr(band2, ik_ibz, spin, wave2, msg) == 0, msg) 1406 ug1 => wave1%ug 1407 ug2 => wave2%ug 1408 1409 wfd_xdotc = xdotc(npw_k*Wfd%nspinor,ug1,1,ug2,1) 1410 if (istwf_k>1) then 1411 wfd_xdotc=two*DBLE(wfd_xdotc) 1412 if (istwf_k==2) wfd_xdotc = wfd_xdotc-CONJG(ug1(1))*ug2(1) 1413 end if 1414 1415 ! Paw on-site term. 1416 if (Wfd%usepaw==1) then 1417 ! Avoid the computation if Cprj are already in memory with the correct order. 1418 if (wave1%has_cprj == WFD_STORED .and. wave1%cprj_order == CPR_RANDOM .and. & 1419 wave2%has_cprj == WFD_STORED .and. wave2%cprj_order == CPR_RANDOM) then 1420 1421 pawovlp = paw_overlap(wave1%Cprj, wave2%Cprj,& 1422 Cryst%typat,Pawtab,spinor_comm=Wfd%MPI_enreg%comm_spinor) 1423 wfd_xdotc = wfd_xdotc + CMPLX(pawovlp(1),pawovlp(2), kind=gwpc) 1424 1425 else 1426 ! Compute Cprj 1427 ABI_MALLOC(Cp1,(Wfd%natom,Wfd%nspinor)) 1428 call pawcprj_alloc(Cp1,0,Wfd%nlmn_atm) 1429 ABI_MALLOC(Cp2,(Wfd%natom,Wfd%nspinor)) 1430 call pawcprj_alloc(Cp2,0,Wfd%nlmn_atm) 1431 1432 call wfd%get_cprj(band1,ik_ibz,spin,Cryst,Cp1,sorted=.FALSE.) 1433 call wfd%get_cprj(band2,ik_ibz,spin,Cryst,Cp2,sorted=.FALSE.) 1434 1435 pawovlp = paw_overlap(Cp1,Cp2,Cryst%typat,Pawtab,spinor_comm=Wfd%MPI_enreg%comm_spinor) 1436 wfd_xdotc = wfd_xdotc + CMPLX(pawovlp(1),pawovlp(2), kind=gwpc) 1437 1438 call pawcprj_free(Cp1) 1439 ABI_FREE(Cp1) 1440 call pawcprj_free(Cp2) 1441 ABI_FREE(Cp2) 1442 end if 1443 end if 1444 1445 end function wfd_xdotc
m_wfd/wfdgw_bands_of_rank [ Functions ]
NAME
wfdgw_bands_of_rank
FUNCTION
Return the list of band index of the ug owned by a given processor at given (k,s).
INPUTS
Wfd rank=The MPI rank of the processor. ik_ibz=Index of the k-point in the IBZ spin=spin index
OUTPUT
how_manyb=The number of bands owned by this node rank_band_list(Wfd%mband)=The first how_manyb values are the bands treated by the node.
SOURCE
2569 subroutine wfdgw_bands_of_rank(Wfd,rank,ik_ibz,spin,how_manyb,rank_band_list) 2570 2571 !Arguments ------------------------------------ 2572 !scalars 2573 integer,intent(in) :: ik_ibz,spin,rank 2574 integer,intent(out) :: how_manyb 2575 class(wfdgw_t),intent(in) :: Wfd 2576 !arrays 2577 integer,intent(out) :: rank_band_list(Wfd%mband) 2578 2579 !Local variables ------------------------------ 2580 !scalars 2581 integer :: band 2582 logical :: it_has 2583 2584 !************************************************************************ 2585 2586 how_manyb=0; rank_band_list=-1 2587 do band=1,Wfd%nband(ik_ibz,spin) 2588 it_has = Wfd%rank_has_ug(rank, band, ik_ibz, spin) 2589 if (it_has) then 2590 how_manyb = how_manyb +1 2591 rank_band_list(how_manyb)=band 2592 end if 2593 end do 2594 2595 end subroutine wfdgw_bands_of_rank
m_wfd/wfdgw_bks_distrb [ Functions ]
NAME
wfdgw_bks_distrb
FUNCTION
Build a local logical table indexed by bands, k-points and spin that defines the distribution of the load inside the loops according to the availability of the ug.
INPUTS
Wfd<wfd_t>= [bks_mask(Wfd%mband,Wfd%nkibz,Wfd%nsppol)]=Mask used to skip selecter (b,k,s) entries. [got(Wfd%nproc)]=The number of tasks already assigned to the nodes.
OUTPUT
bks_distrbk(Wfd%mband,Wfd%nkibz,Wfd%nsppol)=Global table with the rank of the node treating (b,k,s)
SOURCE
3167 subroutine wfdgw_bks_distrb(Wfd, bks_distrb, got, bks_mask) 3168 3169 !Arguments ------------------------------------ 3170 !scalars 3171 class(wfdgw_t),intent(in) :: Wfd 3172 !arrays 3173 integer,intent(out) :: bks_distrb(Wfd%mband,Wfd%nkibz,Wfd%nsppol) 3174 integer,optional,intent(inout) :: got(Wfd%nproc) 3175 logical,optional,intent(in) :: bks_mask(Wfd%mband,Wfd%nkibz,Wfd%nsppol) 3176 3177 !Local variables ------------------------------ 3178 !scalars 3179 integer :: ik_ibz,spin,band,how_many,idle 3180 character(len=500) :: msg 3181 !arrays 3182 integer :: get_more(Wfd%nproc),proc_ranks(Wfd%nproc) 3183 logical :: rank_mask(Wfd%nproc) 3184 3185 !************************************************************************ 3186 3187 get_more=0; if (present(got)) get_more=got 3188 3189 ! Initialize the table here to avoid problem with the cycle instruction below. 3190 bks_distrb = xmpi_undefined_rank 3191 3192 do spin=1,Wfd%nsppol 3193 do ik_ibz=1,Wfd%nkibz 3194 do band=1,Wfd%nband(ik_ibz,spin) 3195 if (present(bks_mask)) then 3196 if (.not.bks_mask(band, ik_ibz, spin)) CYCLE 3197 end if 3198 3199 call wfdgw_who_has_ug(Wfd, band, ik_ibz, spin, how_many, proc_ranks) 3200 3201 if (how_many == 1) then 3202 ! I am the only one owing this band. Add it to list. 3203 bks_distrb(band, ik_ibz, spin) = proc_ranks(1) 3204 3205 else if (how_many>1) then 3206 ! This band is duplicated. Assign it trying to obtain a good load distribution. 3207 rank_mask=.FALSE.; rank_mask(proc_ranks(1:how_many)+1)=.TRUE. 3208 idle = imin_loc(get_more,mask=rank_mask) 3209 get_more(idle) = get_more(idle) + 1 3210 bks_distrb(band,ik_ibz,spin) = proc_ranks(idle) 3211 3212 else 3213 call wfd%dump_errinfo() 3214 write(msg,'(a,3(i0,1x))')" Nobody has (band, ik_ibz, spin): ",band,ik_ibz,spin 3215 ABI_ERROR(msg) 3216 end if 3217 3218 end do 3219 end do 3220 end do 3221 3222 if (present(got)) got=get_more 3223 3224 end subroutine wfdgw_bks_distrb
m_wfd/wfdgw_copy [ Functions ]
NAME
wfdgw_copy
FUNCTION
Copy a wfd_t data type.
SOURCE
1178 subroutine wfdgw_copy(Wfd_in, Wfd_out) 1179 1180 !Arguments ------------------------------------ 1181 class(wfdgw_t),intent(inout) :: Wfd_in,Wfd_out 1182 1183 !Local variables ------------------------------ 1184 !scalars 1185 integer :: band, ik_ibz, spin, cnt_s, cnt_b, ib, ik, is 1186 1187 !************************************************************************ 1188 1189 !@wfd_t 1190 call deep_copy(Wfd_in%debug_level ,Wfd_out%debug_level) 1191 call deep_copy(Wfd_in%lmnmax ,Wfd_out%lmnmax) 1192 call deep_copy(Wfd_in%mband ,Wfd_out%mband) 1193 call deep_copy(Wfd_in%mgfft ,Wfd_out%mgfft) 1194 call deep_copy(Wfd_in%natom ,Wfd_out%natom) 1195 call deep_copy(Wfd_in%nfft ,Wfd_out%nfft) 1196 call deep_copy(Wfd_in%nfftot ,Wfd_out%nfftot) 1197 call deep_copy(Wfd_in%nkibz ,Wfd_out%nkibz) 1198 call deep_copy(Wfd_in%nspden ,Wfd_out%nspden) 1199 call deep_copy(Wfd_in%nspinor ,Wfd_out%nspinor) 1200 call deep_copy(Wfd_in%nsppol ,Wfd_out%nsppol) 1201 call deep_copy(Wfd_in%ntypat ,Wfd_out%ntypat) 1202 call deep_copy(Wfd_in%paral_kgb ,Wfd_out%paral_kgb) 1203 call deep_copy(Wfd_in%usepaw ,Wfd_out%usepaw) 1204 call deep_copy(Wfd_in%prtvol ,Wfd_out%prtvol) 1205 call deep_copy(Wfd_in%pawprtvol ,Wfd_out%pawprtvol) 1206 call deep_copy(Wfd_in%usewvl ,Wfd_out%usewvl) 1207 call deep_copy(Wfd_in%comm ,Wfd_out%comm) 1208 call deep_copy(Wfd_in%master ,Wfd_out%master) 1209 call deep_copy(Wfd_in%my_rank ,Wfd_out%my_rank) 1210 call deep_copy(Wfd_in%nproc ,Wfd_out%nproc) 1211 call deep_copy(Wfd_in%my_nspins ,Wfd_out%my_nspins) 1212 call deep_copy(Wfd_in%rfft_is_symok ,Wfd_out%rfft_is_symok) 1213 call deep_copy(Wfd_in%dilatmx ,Wfd_out%dilatmx) 1214 call deep_copy(Wfd_in%ecut ,Wfd_out%ecut) 1215 call deep_copy(Wfd_in%ecutsm ,Wfd_out%ecutsm) 1216 1217 !arrays 1218 Wfd_out%ngfft =Wfd_in%ngfft 1219 Wfd_out%nloalg=Wfd_in%nloalg 1220 1221 call alloc_copy(Wfd_in%my_nkspin ,Wfd_out%my_nkspin) 1222 call alloc_copy(Wfd_in%irottb ,Wfd_out%irottb) 1223 call alloc_copy(Wfd_in%istwfk ,Wfd_out%istwfk) 1224 call alloc_copy(Wfd_in%nband ,Wfd_out%nband) 1225 call alloc_copy(Wfd_in%indlmn ,Wfd_out%indlmn) 1226 call alloc_copy(Wfd_in%nlmn_atm ,Wfd_out%nlmn_atm) 1227 call alloc_copy(Wfd_in%nlmn_sort ,Wfd_out%nlmn_sort) 1228 call alloc_copy(Wfd_in%nlmn_type ,Wfd_out%nlmn_type) 1229 call alloc_copy(Wfd_in%npwarr ,Wfd_out%npwarr) 1230 call alloc_copy(Wfd_in%kibz ,Wfd_out%kibz) 1231 call alloc_copy(Wfd_in%bks2wfd ,Wfd_out%bks2wfd) 1232 call alloc_copy(Wfd_in%bks_tab ,Wfd_out%bks_tab) 1233 call alloc_copy(Wfd_in%ph1d ,Wfd_out%ph1d) 1234 call alloc_copy(Wfd_in%keep_ur ,Wfd_out%keep_ur) 1235 1236 ! types 1237 if (size(Wfd_in%Kdata,DIM=1) /= size(Wfd_out%Kdata,DIM=1)) then 1238 ABI_REMALLOC(Wfd_out%Kdata, (Wfd_out%nkibz)) 1239 end if 1240 1241 call kdata_copy(Wfd_in%Kdata, Wfd_out%Kdata) 1242 1243 ! Allocate ragged array. 1244 ABI_MALLOC(wfd_out%s, (wfd_out%my_nspins)) 1245 cnt_s = 0 1246 do spin=1,wfd_out%nsppol 1247 if (wfd_out%my_nkspin(spin) > 0) then 1248 cnt_s = cnt_s + 1 1249 ABI_MALLOC(wfd_out%s(cnt_s)%k, (wfd_out%my_nkspin(spin))) 1250 do ik=1,wfd_out%my_nkspin(spin) 1251 cnt_b = size(wfd_out%s(cnt_s)%k(ik)%b) 1252 ABI_MALLOC(wfd_out%s(cnt_s)%k(ik)%b, (cnt_b)) 1253 end do 1254 end if 1255 end do 1256 1257 ! Copy waves 1258 do spin=1,wfd_in%nsppol 1259 do ik_ibz=1,wfd_in%nkibz 1260 do band=1,wfd_in%nband(ik_ibz, spin) 1261 ib = wfd_in%bks2wfd(1, band, ik_ibz, spin) 1262 ik = wfd_in%bks2wfd(2, band, ik_ibz, spin) 1263 is = wfd_in%bks2wfd(3, band, ik_ibz, spin) 1264 if (ib /= 0) wfd_out%s(is)%k(ik)%b(ib) = wfd_in%s(is)%k(ik)%b(ib)%copy() 1265 end do 1266 end do 1267 end do 1268 1269 call copy_mpi_enreg(Wfd_in%MPI_enreg, Wfd_out%MPI_enreg) 1270 1271 end subroutine wfdgw_copy
m_wfd/wfdgw_distribute_bands [ Functions ]
NAME
wfdgw_distribute_bands
FUNCTION
Distribute a set of bands taking into account the distribution of the ug.
INPUTS
band=the index of the band.
ik_ibz=Index of the k-point in the IBZ
spin=spin index
[got(Wfd%nproc)]=The number of tasks already assigned to the nodes.
[bmask(Wfd%mband)]=The routine will raise an error if one band index
is not treated by any processor. bmask can be used to select the subset of
indices that are expected to be available.
OUTPUT
my_nband=The number of bands that will be treated by this node. my_band_list(1:my_nband)=The band indices for this node
SOURCE
2909 subroutine wfdgw_distribute_bands(Wfd,ik_ibz,spin,my_nband,my_band_list,got,bmask) 2910 2911 !Arguments ------------------------------------ 2912 !scalars 2913 integer,intent(in) :: ik_ibz,spin 2914 integer,intent(out) :: my_nband 2915 class(wfdgw_t),intent(in) :: Wfd 2916 !arrays 2917 integer,intent(out) :: my_band_list(Wfd%mband) 2918 integer,optional,intent(inout) :: got(Wfd%nproc) 2919 logical,optional,intent(in) :: bmask(Wfd%mband) 2920 2921 !Local variables ------------------------------ 2922 !scalars 2923 integer :: band,how_many,idle 2924 character(len=500) :: msg 2925 !arrays 2926 integer :: proc_ranks(Wfd%nproc),get_more(Wfd%nproc) 2927 logical :: rank_mask(Wfd%nproc) 2928 2929 !************************************************************************ 2930 2931 my_nband=0; my_band_list=0 2932 get_more=0; if (present(got)) get_more = got 2933 2934 do band=1,Wfd%nband(ik_ibz,spin) 2935 if (present(bmask)) then 2936 if (.not.bmask(band)) CYCLE 2937 end if 2938 2939 call wfdgw_who_has_ug(Wfd, band, ik_ibz, spin, how_many, proc_ranks) 2940 2941 if (how_many == 1) then 2942 ! I am the only one owing this band. Add it to list. 2943 if (proc_ranks(1) == Wfd%my_rank) then 2944 my_nband=my_nband + 1 2945 my_band_list(my_nband) = band 2946 end if 2947 else if (how_many > 1) then 2948 ! This band is duplicated. Assign it trying to obtain a good load distribution. 2949 rank_mask=.FALSE.; rank_mask(proc_ranks(1:how_many)+1)=.TRUE. 2950 idle = imin_loc(get_more, mask=rank_mask) 2951 get_more(idle) = get_more(idle) + 1 2952 if (Wfd%my_rank==idle-1) then 2953 my_nband=my_nband + 1 2954 my_band_list(my_nband) = band 2955 end if 2956 else 2957 write(msg,'(a,3(i0,1x))')" No processor has (band, ik_ibz, spin): ",band,ik_ibz,spin 2958 ABI_ERROR(msg) 2959 end if 2960 end do 2961 2962 if (present(got)) got = get_more 2963 2964 end subroutine wfdgw_distribute_bands
m_wfd/wfdgw_distribute_bbp [ Functions ]
NAME
wfdgw_distribute_bbp
FUNCTION
Distribute a set of (b,b') indices taking into account the MPI distribution of the ug. It is used to calculate matrix elements of the form <b,k,s|O|b',k,s>
INPUTS
Wfd<wfd_t>= ik_ibz=The index of the k-point in the IBZ. spin=Spin index. allup=String used to select or not the upper triangle. Possible values: "All" =Entire (b,b') matrix will be distributed. "Upper"=Only the upper triangle is distributed. [got(%nproc)]=The number of tasks already assigned to the nodes. Used to optimize the work load. Be careful when this routine is called inside several loops since each node should call the routine at each iteration with the same (local) copy of got so that bbp_distrb will assume the same value on each node. [bbp_mask(%mband,%mband)]= mask used to select a subset of (b,b') indices.
OUTPUT
my_nbbp=The number of (b,b') indices treated by this node. bbp_distrb(%mband%mband)=The rank of the node that will treat (b,b').
SOURCE
3388 subroutine wfdgw_distribute_bbp(Wfd,ik_ibz,spin,allup,my_nbbp,bbp_distrb,got,bbp_mask) 3389 3390 !Arguments ------------------------------------ 3391 !scalars 3392 integer,intent(in) :: ik_ibz,spin 3393 integer,intent(out) :: my_nbbp 3394 class(wfdgw_t),intent(in) :: Wfd 3395 character(len=*),intent(in) :: allup 3396 !arrays 3397 integer,intent(out) :: bbp_distrb(Wfd%mband,Wfd%mband) 3398 integer,optional,intent(inout) :: got(Wfd%nproc) 3399 logical,optional,intent(in) :: bbp_mask(Wfd%mband,Wfd%mband) 3400 3401 !Local variables ------------------------------ 3402 !arrays 3403 integer :: loc_got(Wfd%nproc) 3404 3405 !************************************************************************ 3406 3407 ! Just a wrapper around wfdgw_distribute_kb_kpbp. 3408 loc_got=0; if (present(got)) loc_got = got 3409 3410 if (present(bbp_mask)) then 3411 call wfd%distribute_kb_kpbp(ik_ibz,ik_ibz,spin,allup,my_nbbp,bbp_distrb,loc_got,bbp_mask) 3412 else 3413 call wfd%distribute_kb_kpbp(ik_ibz,ik_ibz,spin,allup,my_nbbp,bbp_distrb,loc_got) 3414 end if 3415 3416 end subroutine wfdgw_distribute_bbp
m_wfd/wfdgw_distribute_kb_kpbp [ Functions ]
NAME
wfdgw_distribute_kb_kpbp
FUNCTION
This routines distributes as set of (b,b') indices taking into account the MPI distribution of the ug. It is used to calculate matrix elements of the form <b,k,s|O|b',k',s>
INPUTS
Wfd<wfd_t>= ik_ibz =The index of the k-point k in the IBZ. ikp_ibz=The index of the k-point k' in the IBZ. spin=Spin index. allup=String used to select the upper triangle of the (b,b') matrix. Possible values: "All" =Entire (b,b') matrix will be distributed. "Upper"=Only the upper triangle is distributed. [got(%nproc)]=The number of tasks already assigned to the nodes. Used to optimize the distribution of the tasks. Be careful when this routine is called inside several loops since each node should call the routine at each iteration with the same (local) copy of got so that bbp_distrb will assume the same value on each node. [bbp_mask(%mband,%mband)]= mask used to select a subset of (b,b') indices.
OUTPUT
my_nbbp=The number of (b,b') indices treated by this node. bbp_distrb(%mband%mband)=The rank of the node that will treat (b,b').
SOURCE
3448 subroutine wfdgw_distribute_kb_kpbp(Wfd, ik_ibz, ikp_ibz, spin, allup, my_nbbp, bbp_distrb, & 3449 got, bbp_mask) ! optional 3450 3451 !Arguments ------------------------------------ 3452 !scalars 3453 integer,intent(in) :: ik_ibz,ikp_ibz,spin 3454 integer,intent(out) :: my_nbbp 3455 class(wfdgw_t),intent(in) :: Wfd 3456 character(len=*),intent(in) :: allup 3457 !arrays 3458 integer,intent(out) :: bbp_distrb(Wfd%mband,Wfd%mband) 3459 integer,optional,intent(inout) :: got(Wfd%nproc) 3460 logical,optional,intent(in) :: bbp_mask(Wfd%mband,Wfd%mband) 3461 3462 !Local variables ------------------------------ 3463 !scalars 3464 integer :: my_nband,ib1,ib2,pcb2,pcb1,howmany_b,howmany_bp,workload_min 3465 integer :: rank,ncpus,idle,b1_stop,ierr 3466 character(len=500) :: msg 3467 !arrays 3468 integer :: rank_bandlist_k(Wfd%mband),rank_bandlist_kp(Wfd%mband) 3469 integer :: get_more(Wfd%nproc),my_band_list_k(Wfd%mband) 3470 integer,allocatable :: whocan_k(:,:),whocan_kp(:,:) 3471 logical :: b_mask(Wfd%mband) 3472 3473 !************************************************************************ 3474 3475 ABI_MALLOC_OR_DIE(whocan_k ,(Wfd%mband,Wfd%nproc), ierr) 3476 ABI_MALLOC_OR_DIE(whocan_kp,(Wfd%mband,Wfd%nproc), ierr) 3477 whocan_k =0 ! Will be set to 1 if this node can calculate something containing (k,b) 3478 whocan_kp=0 ! Will be set to 1 if this node can calculate something containing (kp,bp) 3479 3480 do rank=0,Wfd%nproc-1 3481 3482 call wfd%bands_of_rank(rank,ik_ibz ,spin,howmany_b, rank_bandlist_k ) 3483 do pcb1=1,howmany_b 3484 ib1 = rank_bandlist_k(pcb1) 3485 whocan_k(ib1,rank+1) = 1 3486 end do 3487 3488 call wfd%bands_of_rank(rank,ikp_ibz,spin,howmany_bp,rank_bandlist_kp) 3489 do pcb2=1,howmany_bp 3490 ib2 = rank_bandlist_kp(pcb2) 3491 whocan_kp(ib2,rank+1) = 1 3492 end do 3493 3494 end do 3495 3496 get_more=0; if (present(got)) get_more=got 3497 b1_stop=Wfd%nband(ik_ibz,spin) 3498 3499 bbp_distrb = xmpi_undefined_rank 3500 3501 do ib2=1,Wfd%nband(ikp_ibz,spin) 3502 b_mask = .TRUE.; if (present(bbp_mask)) b_mask = bbp_mask(:,ib2) 3503 if (ANY(b_mask)) then 3504 my_nband=0; my_band_list_k=0 3505 ! Only the upper triangle of the (b1,b2) matrix. 3506 if (firstchar(allup, ["U","u"])) b1_stop = MIN(ib2,Wfd%nband(ik_ibz,spin)) 3507 3508 do ib1=1,b1_stop 3509 if (b_mask(ib1)) then 3510 ! 3511 ! find which CPUs can do the calculation (k,b)->(kp,bp) 3512 ! find the one which is less busy 3513 ncpus=0 3514 workload_min=HUGE(0) 3515 do rank=0,Wfd%nproc-1 3516 if( whocan_k(ib1,rank+1)==1 .AND. whocan_kp(ib2,rank+1)==1 ) then 3517 ncpus=ncpus+1 3518 if( get_more(rank+1) < workload_min ) then 3519 idle=rank+1 3520 workload_min=get_more(idle) 3521 end if 3522 3523 end if 3524 end do 3525 3526 if(ncpus>0) then 3527 bbp_distrb(ib1,ib2)=idle-1 3528 get_more(idle) = get_more(idle) + 1 3529 3530 else 3531 call wfd%dump_errinfo() 3532 write(msg,'(a,5(i0,1x))')" Nobody has (band1, ik_ibz) (band2, ikp_ibz) spin: ",ib1,ik_ibz,ib2,ikp_ibz,spin 3533 ABI_ERROR(msg) 3534 end if 3535 3536 end if 3537 end do ! ib1 3538 end if 3539 end do ! ib2 3540 3541 ABI_FREE(whocan_k) 3542 ABI_FREE(whocan_kp) 3543 3544 my_nbbp = COUNT(bbp_distrb==Wfd%my_rank) 3545 if (present(got)) got=get_more 3546 3547 end subroutine wfdgw_distribute_kb_kpbp
m_wfd/wfdgw_get_nl_me [ Functions ]
NAME
wfdgw_get_nl_me
FUNCTION
Compute the non-local potential using the Wfd information.
INPUTS
cryst<crystal_t>= data type gathering info on symmetries and unit cell psps<pseudopotential_type>=variables related to pseudopotentials pawtab(psps%ntypat) <type(pawtab_type)>=paw tabulated starting data paw_ij(natom)<type(paw_ij_type)>=data structure containing PAW arrays given on (i,j) channels.
OUTPUT
SOURCE
5109 subroutine wfdgw_get_nl_me(Wfd, cryst, psps, pawtab, bks_mask, nl_bks) 5110 5111 !Arguments ------------------------------------ 5112 !scalars 5113 class(wfdgw_t),target,intent(inout) :: Wfd 5114 type(crystal_t),intent(in) :: cryst 5115 type(pseudopotential_type),intent(in) :: psps 5116 ! arrays 5117 logical,intent(in) :: bks_mask(Wfd%mband, Wfd%nkibz, Wfd%nsppol) 5118 real(dp),allocatable,intent(out) :: nl_bks(:, :, :) 5119 type(Pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 5120 5121 !Local variables ------------------------------ 5122 !scalars 5123 integer,parameter :: nspinor1=1,nspden1=1,nsppol1=1,spin1=1 5124 integer,parameter :: ndat1=1,nnlout1=1,tim_nonlop0=0,idir0=0 5125 integer :: natom,ib,ispin,ik_ibz,npw_k,istwf_k,nkpg 5126 integer :: choice,cpopt,paw_opt,signs,ierr 5127 character(len=500) :: msg 5128 type(gs_hamiltonian_type) :: ham_k 5129 type(wave_t),pointer :: wave 5130 !arrays 5131 integer :: bks_distrb(Wfd%mband, Wfd%nkibz, Wfd%nsppol) 5132 integer, ABI_CONTIGUOUS pointer :: kg_k(:,:) 5133 real(dp) :: kpoint(3),enlout(1) 5134 real(dp),allocatable :: kpg_k(:,:),vnl_psi(:,:),vectin(:,:) 5135 real(dp) :: opaw_psi(1,1) 5136 real(dp),ABI_CONTIGUOUS pointer :: ffnl_k(:,:,:,:),ph3d_k(:,:,:) 5137 complex(gwpc),ABI_CONTIGUOUS pointer :: ug1(:) 5138 type(pawcprj_type),allocatable :: cprj(:,:) 5139 5140 !************************************************************************ 5141 5142 ABI_CHECK(Wfd%paral_kgb == 0, "paral_kgb not coded") 5143 5144 natom = cryst%natom 5145 5146 signs = 1 ! => get non-local energy in G-space. 5147 choice = 1 ! => <G|V_nonlocal|vectin>. 5148 cpopt =-1; paw_opt= 0 5149 if (Wfd%usepaw==1) then 5150 ABI_ERROR("The construction of the non-local contribution is not tested/implemented for usepaw==1!") 5151 end if 5152 ! Initialize the Hamiltonian on the coarse FFT mesh. 5153 call init_hamiltonian(ham_k, psps, pawtab, nspinor1, nsppol1, nspden1, natom, cryst%typat, cryst%xred, & 5154 Wfd%nfft, Wfd%mgfft, Wfd%ngfft, cryst%rprimd, Wfd%nloalg) 5155 5156 ! Continue to prepare the GS Hamiltonian (note spin1) 5157 call ham_k%load_spin(spin1, with_nonlocal=.true.) 5158 5159 ! Distribute (b, k, s) states. 5160 call Wfd%bks_distrb(bks_distrb, bks_mask=bks_mask) 5161 5162 ABI_CALLOC(nl_bks, (Wfd%mband, Wfd%nkibz, Wfd%nsppol)) 5163 nl_bks(:,:,:) = czero 5164 5165 write(std_out,'(a)') " " 5166 call wrtout(std_out,sjoin(" Will calculate ",itoa(count(bks_mask))," <b,k,s|Vnl|b,k,s> matrix elements in wfdgw_get_nl_me.")) 5167 do ispin=1,Wfd%nsppol 5168 if (ispin/=1) then 5169 ABI_WARNING("In the construction of the non-local contribution, the case nsppol/=1 is not tested.") 5170 end if 5171 do ik_ibz=1,Wfd%nkibz 5172 if (all(bks_distrb(:, ik_ibz, ispin) /= Wfd%my_rank)) cycle 5173 5174 kpoint = Wfd%kibz(:, ik_ibz) 5175 istwf_k = Wfd%istwfk(ik_ibz) 5176 npw_k = Wfd%Kdata(ik_ibz)%npw 5177 kg_k => Wfd%kdata(ik_ibz)%kg_k 5178 ffnl_k => Wfd%Kdata(ik_ibz)%fnl_dir0der0 5179 ph3d_k => Wfd%Kdata(ik_ibz)%ph3d 5180 5181 ABI_MALLOC(vectin, (2, npw_k * nspinor1)) 5182 ABI_MALLOC(vnl_psi, (2, npw_k * nspinor1)) 5183 vectin=zero 5184 vnl_psi=zero 5185 ! Compute (k+G) vectors (only if psps%useylm=1) 5186 nkpg = 3 * Wfd%nloalg(3) 5187 ABI_MALLOC(kpg_k, (npw_k, nkpg)) 5188 if (nkpg > 0) then 5189 call mkkpg(kg_k, kpg_k, kpoint, nkpg, npw_k) 5190 end if 5191 5192 ! Load k-dependent part in the Hamiltonian datastructure 5193 call ham_k%load_k(kpt_k=kpoint,istwf_k=istwf_k,npw_k=npw_k,kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl_k,& 5194 & ph3d_k=ph3d_k,compute_ph3d=(Wfd%paral_kgb/=1),compute_gbound=(Wfd%paral_kgb/=1)) 5195 5196 ! ======================================================== 5197 ! ==== Compute nonlocal form factors ffnl at all (k+G) ==== 5198 ! ======================================================== 5199 ! Calculate <G|Vnl|psi> for this k-point 5200 do ib=1,Wfd%nband(ik_ibz, ispin) 5201 if (bks_distrb(ib, ik_ibz, ispin) /= Wfd%my_rank) cycle 5202 5203 ABI_CHECK(Wfd%get_wave_ptr(ib, ik_ibz, ispin, wave, msg) == 0, msg) 5204 ug1 => wave%ug 5205 ! Input wavefunction coefficients <G|Cnk>. 5206 ! vectin, (2, npw_k * nspinor1)) 5207 vectin(:,:) = zero 5208 if (ispin == 1) then 5209 vectin(1, 1:npw_k) = real(ug1) 5210 vectin(2, 1:npw_k) = aimag(ug1) 5211 else 5212 vectin(1, npw_k+1:) = real(ug1) 5213 vectin(2, npw_k+1:) = aimag(ug1) 5214 end if 5215 5216 if (Wfd%usepaw == 1) then 5217 call Wfd%get_cprj(ib, ik_ibz, ispin, cryst, cprj, sorted=.True.) 5218 end if 5219 5220 ! Compute nonlocal band energy. We could use nonlop to get matrix elements but we only need diagonal contributions. 5221 ! enlout(1) is the band energy 5222 call nonlop(choice, cpopt, cprj, enlout, ham_k, idir0, (/zero/), Wfd%mpi_enreg, ndat1, nnlout1, & 5223 paw_opt, signs, opaw_psi, tim_nonlop0, vectin, vnl_psi) 5224 nl_bks(ib, ik_ibz, ispin) = enlout(1) 5225 end do ! ib 5226 5227 ABI_FREE(vectin) 5228 ABI_FREE(vnl_psi) 5229 ABI_FREE(kpg_k) 5230 end do ! ik_ibz 5231 end do ! ispin 5232 5233 call xmpi_sum(nl_bks, Wfd%comm, ierr) 5234 5235 call ham_k%free() 5236 5237 end subroutine wfdgw_get_nl_me
m_wfd/wfdgw_iterator_bks [ Functions ]
NAME
wfdgw_iterator_bks
FUNCTION
Iterator used to loop over bands, k-points and spin indices taking into account the distribution of the ug.
INPUTS
Wfd<wfd_t>= bks_mask(Wfd%mband.Wfd%nkibz,Wfd%nsppol)= mask used to select the (b,k,s) indices.
OUTPUT
iter_bks<iter2_t>=Iterator over the bands treated by this node for each k-point and spin.
SOURCE
3115 type(iter2_t) function wfdgw_iterator_bks(Wfd, bks_mask) result(iter_bks) 3116 3117 !Arguments ------------------------------------ 3118 !scalars 3119 class(wfdgw_t),intent(in) :: Wfd 3120 !arrays 3121 logical,optional,intent(in) :: bks_mask(Wfd%mband,Wfd%nkibz,Wfd%nsppol) 3122 3123 !Local variables ------------------------------ 3124 !scalars 3125 integer :: ik_ibz,spin,my_nband 3126 !arrays 3127 integer :: my_band_list(Wfd%mband) 3128 3129 !************************************************************************ 3130 3131 call iter_alloc(iter_bks,(/Wfd%nkibz,Wfd%nsppol/)) 3132 3133 do spin=1,Wfd%nsppol 3134 do ik_ibz=1,Wfd%nkibz 3135 if (present(bks_mask)) then 3136 call wfd%distribute_bands(ik_ibz,spin,my_nband,my_band_list,bmask=bks_mask(:,ik_ibz,spin)) 3137 else 3138 call wfd%distribute_bands(ik_ibz,spin,my_nband,my_band_list) 3139 end if 3140 call iter_push(iter_bks,ik_ibz,spin,my_band_list(1:my_nband)) 3141 end do 3142 end do 3143 3144 end function wfdgw_iterator_bks
m_wfd/wfdgw_mkrho [ Functions ]
NAME
wfdgw_mkrho
FUNCTION
Calculate the charge density on the fine FFT grid in real space.
INPUTS
Wfd<wfd_t)=datatype gathering info on the wavefunctions. ngfftf(18)=array containing all the information for the "fine" FFT. Cryst<crystal_t> Info on the crystalline structure optcalc=option for calculation. If =0 (default value) perform calculation of electronic density. If =1, perform calculation of kinetic energy density. In both cases, the result is returned in rhor. Psps<type(pseudopotential_type)>=variables related to pseudopotentials nfftf=Total number of points on the fine FFT grid (for this processor) [optcalc]=Optional option used to calculate the kinetic energy density. Defaults to 0.
OUTPUT
rhor(nfftf,nspden)=The density in the real space on the fine FFT grid. If nsppol==2, total charge in first half, spin-up component in second half. (for non-collinear magnetism, first element: total density, 3 next ones: mx,my,mz in units of hbar/2) If optcalc==1 (optional argument, default value is 0), then rhor will actually contains kinetic energy density (taur) instead of electronic density.
NOTES
In the case of PAW calculations:
All computations are done on the fine FFT grid.
All variables (nfftf,ngfftf,mgfftf) refer to this fine FFT grid.
All arrays (densities/potentials...) are computed on this fine FFT grid.
Developers have to be careful when introducing others arrays:
they have to be stored on the fine FFT grid.
In the case of norm-conserving calculations:
The mesh is the usual augmented FFT grid to treat correctly the convolution.
SOURCE
5467 subroutine wfdgw_mkrho(wfd, cryst, psps, ebands, ngfftf, nfftf, rhor, & 5468 optcalc) ! optional arguments 5469 5470 !Arguments ------------------------------------ 5471 !scalars 5472 integer,intent(in) :: nfftf 5473 integer,intent(in),optional :: optcalc 5474 type(ebands_t),intent(in) :: ebands 5475 type(crystal_t),intent(in) :: cryst 5476 type(Pseudopotential_type),intent(in) :: psps 5477 class(wfdgw_t),intent(inout) :: wfd 5478 !arrays 5479 integer,intent(in) :: ngfftf(18) 5480 real(dp),intent(out) :: rhor(nfftf, Wfd%nspden) 5481 5482 !Local variables ------------------------------ 5483 !scalars 5484 integer,parameter :: ndat1=1 5485 integer :: cplex,ib,ib_iter,ierr,ik,ir,is,n1,n2,n3,nfftotf 5486 integer :: alpha,nalpha,ipw,myoptcalc 5487 real(dp) :: kpt_cart,kg_k_cart,gp2pi1,gp2pi2,gp2pi3,cwftmp,bks_weight 5488 character(len=500) :: msg 5489 type(wave_t),pointer :: wave 5490 !arrays 5491 integer,allocatable :: irrzon(:,:,:) 5492 real(dp),allocatable :: phnons(:,:,:),rhog(:,:),rhor_down(:),rhor_mx(:),rhor_my(:),cwavef(:,:) 5493 complex(dpc),allocatable :: wfr_x(:),wfr_y(:) 5494 complex(gwpc),allocatable :: gradug(:),work(:) 5495 complex(gwpc),allocatable,target :: wfr(:) 5496 complex(gwpc), ABI_CONTIGUOUS pointer :: cwavef1(:),cwavef2(:) 5497 type(iter2_t) :: Iter_bks 5498 5499 !************************************************************************* 5500 5501 ! Consistency check. 5502 ABI_CHECK(Wfd%nsppol == ebands%nsppol, "Mismatch in nsppol") 5503 5504 if (ANY(ngfftf(1:3) /= Wfd%ngfft(1:3))) call wfd%change_ngfft(Cryst,Psps,ngfftf) 5505 5506 ! Calculate IBZ contribution to the charge density. 5507 ABI_MALLOC(wfr, (nfftf*Wfd%nspinor)) 5508 5509 if (wfd%nspden == 4) then 5510 ABI_MALLOC(wfr_x, (nfftf)) 5511 ABI_MALLOC(wfr_y, (nfftf)) 5512 ABI_MALLOC(rhor_down, (nfftf)) 5513 ABI_MALLOC(rhor_mx, (nfftf)) 5514 ABI_MALLOC(rhor_my, (nfftf)) 5515 rhor_down = zero; rhor_mx = zero; rhor_my = zero 5516 end if 5517 5518 ! Update the (b,k,s) distribution table. 5519 call wfd%update_bkstab() 5520 5521 ! Calculate the unsymmetrized density. 5522 rhor=zero 5523 myoptcalc=0; if (present(optcalc)) myoptcalc=optcalc 5524 nalpha=1; if (myoptcalc==1) nalpha=3 5525 if (myoptcalc == 1 .and. wfd%nspinor == 2) then 5526 ABI_ERROR("kinetic energy density with nspinor == 2 not implemented") 5527 end if 5528 5529 ! Build the iterator that will distribute the states in an automated way. 5530 Iter_bks = wfd%iterator_bks(bks_mask=ABS(ebands%occ)>=tol8) 5531 5532 do alpha=1,nalpha 5533 do is=1,Wfd%nsppol 5534 do ik=1,Wfd%nkibz 5535 do ib_iter=1,iter_len(Iter_bks,ik,is) 5536 ib = iter_yield(Iter_bks,ib_iter,ik,is) 5537 bks_weight = ebands%occ(ib,ik,is) * ebands%wtk(ik) / Cryst%ucvol 5538 5539 call wfd%get_ur(ib,ik,is,wfr) 5540 5541 cwavef1 => wfr(1:nfftf) 5542 if (myoptcalc == 1) then 5543 ABI_MALLOC(gradug,(Wfd%Kdata(ik)%npw)) 5544 ABI_MALLOC(cwavef,(2,Wfd%Kdata(ik)%npw)) 5545 ABI_MALLOC(work,(nfftf)) 5546 5547 ABI_CHECK(wfd%get_wave_ptr(ib, ik, is, wave, msg) == 0, msg) 5548 cwavef(1,:)= REAL(wave%ug(:)) 5549 cwavef(2,:)=AIMAG(wave%ug(:)) 5550 5551 ! Multiplication by 2pi i (k+G)_alpha 5552 gp2pi1=Cryst%gprimd(alpha,1)*two_pi 5553 gp2pi2=Cryst%gprimd(alpha,2)*two_pi 5554 gp2pi3=Cryst%gprimd(alpha,3)*two_pi 5555 kpt_cart=gp2pi1*Wfd%kibz(1,ik)+gp2pi2*Wfd%kibz(2,ik)+gp2pi3*Wfd%kibz(3,ik) 5556 do ipw=1,Wfd%Kdata(ik)%npw 5557 kg_k_cart= gp2pi1*Wfd%Kdata(ik)%kg_k(1,ipw) + & 5558 gp2pi2*Wfd%Kdata(ik)%kg_k(2,ipw) + & 5559 gp2pi3*Wfd%Kdata(ik)%kg_k(3,ipw)+kpt_cart 5560 ! ipwsp=ipw!+(ispinor-1)*Wfd%Kdata(ik)%npw 5561 cwftmp=-cwavef(2,ipw)*kg_k_cart 5562 cwavef(2,ipw)=cwavef(1,ipw)*kg_k_cart 5563 cwavef(1,ipw)=cwftmp 5564 end do 5565 gradug(:)=CMPLX(cwavef(1,:),cwavef(2,:),gwpc) 5566 call fft_ug(Wfd%npwarr(ik),nfftf,Wfd%nspinor,ndat1,Wfd%mgfft,Wfd%ngfft,& 5567 Wfd%istwfk(ik),Wfd%Kdata(ik)%kg_k,Wfd%Kdata(ik)%gbound,gradug,work) 5568 cwavef1(:)=work(:) 5569 ABI_FREE(work) 5570 ABI_FREE(cwavef) 5571 ABI_FREE(gradug) 5572 end if 5573 5574 !$OMP PARALLEL DO 5575 do ir=1,nfftf 5576 rhor(ir,is) = rhor(ir,is) + CONJG(cwavef1(ir)) * cwavef1(ir) * bks_weight 5577 end do 5578 !call cplx_addtorho(n1,n2,n3,n4,n5,n6,ndat,weight_r,ur,rho) 5579 5580 if (wfd%nspinor == 2 .and. wfd%nspden == 1) then 5581 cwavef2 => wfr(1+nfftf:2*nfftf) 5582 do ir=1,nfftf 5583 rhor(ir, 1) = rhor(ir, 1) + CONJG(cwavef2(ir)) * cwavef2(ir) * bks_weight 5584 end do 5585 end if 5586 5587 if (wfd%nspinor == 2 .and. wfd%nspden == 4) then 5588 cwavef2 => wfr(1+nfftf:2*nfftf) 5589 wfr_x(:) = cwavef1(:) + cwavef2(:) ! $(\Psi^{1}+\Psi^{2})$ 5590 wfr_y(:) = cwavef1(:) -j_dpc*cwavef2(:) ! $(\Psi^{1}-i\Psi^{2})$ 5591 !$OMP PARALLEL DO 5592 do ir=1,nfftf 5593 rhor_down(ir) = rhor_down(ir) + CONJG(cwavef2(ir)) * cwavef2(ir) * bks_weight 5594 rhor_mx(ir) = rhor_mx(ir) + CONJG(wfr_x(ir)) * wfr_x(ir) * bks_weight 5595 rhor_my(ir) = rhor_my(ir) + CONJG(wfr_y(ir)) * wfr_y(ir) * bks_weight 5596 end do 5597 end if 5598 5599 end do 5600 end do 5601 end do 5602 5603 end do ! enddo alpha 5604 5605 call iter_free(Iter_bks) 5606 5607 select case (myoptcalc) 5608 case (0) 5609 ! density 5610 if (wfd%nspden == 4) then 5611 rhor(:, 2) = rhor_mx 5612 rhor(:, 3) = rhor_my 5613 rhor(:, 4) = rhor_down 5614 end if 5615 case (1) 5616 ! convention for taur = 1/2 Sum_i |grad phi_i|^2 5617 rhor(:,:)=half*rhor(:,:) 5618 5619 case default 5620 ABI_ERROR(sjoin("Wrong myoptcalc:", itoa(myoptcalc))) 5621 end select 5622 5623 call xmpi_sum(rhor,Wfd%comm,ierr) 5624 5625 ! Symmetrization in G-space implementing also the AFM case 5626 n1=ngfftf(1); n2=ngfftf(2); n3=ngfftf(3); nfftotf=n1*n2*n3 5627 5628 ABI_MALLOC(irrzon,(nfftotf**(1-1/Cryst%nsym),2,(Wfd%nspden/Wfd%nsppol)-3*(Wfd%nspden/4))) 5629 ABI_MALLOC(phnons,(2,nfftotf,(Wfd%nspden/Wfd%nsppol)-3*(Wfd%nspden/4))) 5630 5631 if (Cryst%nsym/=1) then 5632 call irrzg(irrzon,Wfd%nspden,Wfd%nsppol,Cryst%nsym,n1,n2,n3,phnons,Cryst%symafm,Cryst%symrel,Cryst%tnons) 5633 end if 5634 5635 ! Symmetrize rho(r), and pack nspden components following abinit conventions. 5636 cplex=1 5637 ABI_MALLOC(rhog,(2,cplex*nfftf)) 5638 5639 call symrhg(cplex,Cryst%gprimd,irrzon,Wfd%MPI_enreg,nfftf,nfftotf,ngfftf,Wfd%nspden,Wfd%nsppol,& 5640 Cryst%nsym,phnons,rhog,rhor,Cryst%rprimd,Cryst%symafm,Cryst%symrel,Cryst%tnons) 5641 5642 ABI_FREE(rhog) 5643 ABI_FREE(phnons) 5644 ABI_FREE(irrzon) 5645 5646 ! Find and print minimum and maximum total electron density 5647 ! (or total kinetic energy density, or total element of kinetic energy density tensor) and locations 5648 !call wrtout(std_out,'mkrho: echo density (plane-wave part only)','COLL') 5649 !call prtrhomxmn(std_out,wfd%mpi_enreg,nfftf,ngfftf,wfd%nspden,1,rhor,optrhor=optcalc,ucvol=crystl%ucvol) 5650 5651 write(msg,'(a,f9.4)')' planewave contribution to nelect: ',SUM(rhor(:,1))*Cryst%ucvol/nfftf 5652 call wrtout(std_out, msg) 5653 5654 if (Wfd%nspden==4) then 5655 write(msg,'(a,3f9.4)')& 5656 ' mx, my, mz: ',SUM(rhor(:,2))*Cryst%ucvol/nfftf,SUM(rhor(:,3))*Cryst%ucvol/nfftf,SUM(rhor(:,4))*Cryst%ucvol/nfftf 5657 call wrtout(std_out, msg) 5658 end if 5659 5660 ABI_FREE(wfr) 5661 5662 if (Wfd%nspden == 4) then 5663 ABI_FREE(wfr_x) 5664 ABI_FREE(wfr_y) 5665 ABI_FREE(rhor_down) 5666 ABI_FREE(rhor_mx) 5667 ABI_FREE(rhor_my) 5668 end if 5669 5670 end subroutine wfdgw_mkrho
m_wfd/wfdgw_pawrhoij [ Functions ]
NAME
wfdgw_pawrhoij
FUNCTION
Calculate the PAW quantities rhoij (augmentation occupancies)
Remember:for each atom, rho_ij=Sum_{n,k} {occ(n,k)*<Cnk|p_i><p_j|Cnk>}
INPUTS
atindx1(natom)=index table for atoms, inverse of atindx
cprj(natom,nspinor*mband*mkmem*nsppol)= wave functions projected with non-local projectors:
cprj_nk(i)=<p_i|Cnk> where p_i is a non-local projector.
istwfk(nkpt)=parameter that describes the storage of wfs
kptopt=option for the generation of k points
mband=maximum number of bands
natom=number of atoms in cell
nkpt=number of k points
nspinor=number of spinorial components of the wavefunctions
nsppol=1 for unpolarized, 2 for spin-polarized
occ(mband*nkpt*nsppol)=occupation number for each band for each k
pawprtvol=control print volume and debugging output for PAW
SIDE EFFECTS
pawrhoij(natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data On input: arrays dimensions On output: pawrhoij(:)%rhoij_(lmn2_size,nspden)= Sum_{n,k} {occ(n,k)*conjugate[cprj_nk(ii)].cprj_nk(jj)} (non symetrized)
SOURCE
5799 subroutine wfdgw_pawrhoij(Wfd,Cryst,Bst,kptopt,pawrhoij,pawprtvol) 5800 5801 !Arguments --------------------------------------------- 5802 !scalars 5803 integer,intent(in) :: kptopt,pawprtvol 5804 type(crystal_t),intent(in) :: Cryst 5805 class(wfdgw_t),intent(inout) :: Wfd 5806 type(ebands_t),intent(in) :: Bst 5807 !arrays 5808 type(pawrhoij_type),intent(inout) :: pawrhoij(Wfd%natom) 5809 5810 !Local variables --------------------------------------- 5811 !scalars 5812 integer :: cplex,cplex_rhoij,qphase,iatom,band,ik_ibz 5813 integer :: spin,natinc,nband_k,option,lmn2_size,nspden 5814 logical :: use_timerev,use_zeromag 5815 real(dp) :: occup,wtk_k 5816 character(len=500) :: msg 5817 !arrays 5818 !real(dp) :: tsec(2) 5819 character(len=8),parameter :: dspin(6)=(/"up ","down ","dens (n)","magn (x)","magn (y)","magn (z)"/) 5820 type(pawcprj_type),allocatable :: cwaveprj(:,:) 5821 integer :: bks_distrb(Wfd%mband,Wfd%nkibz,Wfd%nsppol) 5822 integer :: got(Wfd%nproc) 5823 logical :: bks_mask(Wfd%mband,Wfd%nkibz,Wfd%nsppol) 5824 5825 !************************************************************************ 5826 5827 ! Allocate temporary cwaveprj storage (sorted by atom type) 5828 ABI_MALLOC(cwaveprj,(Wfd%natom,Wfd%nspinor)) 5829 call pawcprj_alloc(cwaveprj,0,Wfd%nlmn_sort) 5830 5831 ! Initialize output quantities if not already done. 5832 do iatom=1,Wfd%natom 5833 if (pawrhoij(iatom)%use_rhoij_==0) then 5834 cplex_rhoij= pawrhoij(iatom)%cplex_rhoij 5835 qphase = pawrhoij(iatom)%qphase 5836 lmn2_size = pawrhoij(iatom)%lmn2_size 5837 nspden = pawrhoij(iatom)%nspden 5838 ABI_MALLOC(pawrhoij(iatom)%rhoij_,(cplex_rhoij*qphase*lmn2_size,nspden)) 5839 pawrhoij(iatom)%use_rhoij_=1 5840 end if 5841 pawrhoij(iatom)%rhoij_=zero 5842 end do 5843 5844 option=1 5845 use_timerev=(kptopt>0.and.kptopt<3) 5846 use_zeromag=.false. ; if (Wfd%natom>0) use_zeromag=(pawrhoij(1)%nspden==4.and.Wfd%nspden==1) 5847 5848 ! Distribute (b,k,s). 5849 where (ABS(Bst%occ)>tol8) 5850 bks_mask=.TRUE. 5851 else where 5852 bks_mask=.FALSE. 5853 end where 5854 got = 0 5855 5856 call wfd%bks_distrb(bks_distrb,got,bks_mask) 5857 5858 do spin=1,Wfd%nsppol 5859 do ik_ibz=1,Wfd%nkibz 5860 5861 nband_k=Wfd%nband(ik_ibz,spin) 5862 wtk_k=Bst%wtk(ik_ibz) 5863 5864 cplex=2; if (Wfd%istwfk(ik_ibz)>1) cplex=1 5865 5866 do band=1,nband_k 5867 5868 if (bks_distrb(band,ik_ibz,spin) == Wfd%my_rank) then 5869 !locc_test = (abs(Bst%occ(band,ik_ibz,spin))>tol8) 5870 occup = Bst%occ(band,ik_ibz,spin) 5871 5872 ! Extract cprj for current band cwaveprj are sorted by atom type. 5873 call wfd%get_cprj(band,ik_ibz,spin,Cryst,cwaveprj,sorted=.TRUE.) 5874 5875 ! Accumulate contribution from (occupied) current band 5876 !if (locc_test) then 5877 call pawaccrhoij(Cryst%atindx,cplex,cwaveprj,cwaveprj ,0,spin,Wfd%natom,Wfd%natom,& 5878 & Wfd%nspinor,occup,option,pawrhoij,use_timerev,use_zeromag,wtk_k) 5879 !end if 5880 end if 5881 end do !band 5882 5883 end do !ik_ibz 5884 end do !spin 5885 ! 5886 ! Free temporary cwaveprj storage. 5887 call pawcprj_free(cwaveprj) 5888 ABI_FREE(cwaveprj) 5889 ! 5890 !========================================== 5891 ! MPI: need to exchange arrays between procs 5892 ! TODO it should be tested. 5893 call pawrhoij_mpisum_unpacked(pawrhoij,Wfd%comm) 5894 5895 ! Print info. 5896 if (abs(pawprtvol)>=1) then 5897 natinc=1; if(Wfd%natom>1.and.pawprtvol>=0) natinc=Wfd%natom-1 5898 write(msg, '(7a)') ch10," PAW TEST:",ch10,' ========= Values of RHOIJ in wfdgw_pawrhoij =========',ch10 5899 call wrtout(std_out, msg) 5900 do iatom=1,Cryst%natom,natinc 5901 call pawrhoij_print_rhoij(pawrhoij(iatom)%rhoij_,pawrhoij(iatom)%cplex_rhoij,& 5902 pawrhoij(iatom)%qphase,iatom,Cryst%natom,& 5903 unit=std_out,opt_prtvol=pawprtvol) 5904 end do 5905 end if 5906 5907 end subroutine wfdgw_pawrhoij
m_wfd/wfdgw_plot_ur [ Functions ]
NAME
wfdgw_plot_ur
FUNCTION
This routine writes the squared modulus of the wavefunctions in real space to an external files, one for each (k,b,s). File are written in the XSF format (Xcrysden). A subset of (b,k,s) states can be specified via the bks_mask. The routine is MPI parallelized.
INPUTS
Wfd<wfd_t>=Wavefunction descriptor. Cryst<crystal_t>= Information on symmetries and unit cell. Psps<pseudopotential_type>=Pseudopotential info. Pawtab(ntypat*usepaw)<type(pawtab_type)>=PAW tabulated starting data. Pawrad(ntypat*usepaw)<type(pawrad_type)>=paw radial mesh and related data. ngfftf(18)=The FFT mesh used for plotting |wfr|**2, it can differ from the one internally used in Wfd. For example, PAW wavefunctions should be plotted on a much finer FFT mesh. bks_mask(mband,nkibz,nsppol)=logical mask used to select the states to be plotted.
OUTPUT
Output is written on file.
SOURCE
4938 subroutine wfdgw_plot_ur(Wfd,Cryst,Psps,Pawtab,Pawrad,ngfftf,bks_mask) 4939 4940 !Arguments ------------------------------------ 4941 !scalars 4942 type(crystal_t),intent(in) :: Cryst 4943 type(Pseudopotential_type),intent(in) :: Psps 4944 class(wfdgw_t),intent(inout) :: Wfd 4945 !arrays 4946 integer,intent(in) :: ngfftf(18) 4947 logical,target,intent(in) :: bks_mask(Wfd%mband,Wfd%nkibz,Wfd%nsppol) 4948 type(Pawtab_type),intent(in) :: Pawtab(Cryst%ntypat*Wfd%usepaw) 4949 type(Pawrad_type),intent(in) :: Pawrad(Cryst%ntypat*Wfd%usepaw) 4950 4951 !Local variables ------------------------------ 4952 !scalars 4953 integer :: spin,band,ik_ibz,optcut,optgr0,optgr1,optgr2,optrad 4954 integer :: n1,n2,n3,my_nplots,plot,funt,my_nband,cplex 4955 character(len=500) :: msg 4956 character(len=fnlen) :: xsf_fname 4957 !arrays 4958 integer :: got(Wfd%nproc) 4959 integer,allocatable :: l_size_atm(:),my_plot_list(:,:) 4960 integer :: my_band_list(Wfd%mband) 4961 real(dp),allocatable :: data_plot(:) 4962 logical,ABI_CONTIGUOUS pointer :: bmask(:) 4963 complex(gwpc),allocatable :: ur_ae(:),nc_ur(:) 4964 type(Pawfgrtab_type),allocatable :: Pawfgrtab(:) 4965 type(paw_pwaves_lmn_t),allocatable :: Paw_onsite(:) 4966 4967 !************************************************************************ 4968 4969 if (ALL(.not.bks_mask)) RETURN 4970 4971 call wrtout(std_out," Plotting |wfs|^2 ...") 4972 ! 4973 ! Change the FFT mesh if needed because we want u(r) on the ngfftf mesh (pawecutd for PAW). 4974 call wfd%change_ngfft(Cryst,Psps,ngfftf) 4975 n1 = ngfftf(1); n2 = ngfftf(2); n3 = ngfftf(3) 4976 4977 ! Distribute the plots among the nodes taking into account the distribution of the waves. 4978 ! my_plot_list gives the list of (b,k,s) states plotted by this node. 4979 ABI_MALLOC(my_plot_list,(3,Wfd%mband*Wfd%nkibz*Wfd%nsppol)) 4980 4981 my_nplots=0; got=0 4982 do spin=1,Wfd%nsppol 4983 do ik_ibz=1,Wfd%nkibz 4984 bmask => bks_mask(:,ik_ibz,spin) 4985 call wfd%distribute_bands(ik_ibz,spin,my_nband,my_band_list,got,bmask) 4986 4987 if (my_nband>0) then 4988 my_plot_list(1,my_nplots+1:my_nplots+my_nband) = my_band_list(1:my_nband) 4989 my_plot_list(2,my_nplots+1:my_nplots+my_nband) = ik_ibz 4990 my_plot_list(3,my_nplots+1:my_nplots+my_nband) = spin 4991 my_nplots = my_nplots + my_nband 4992 end if 4993 end do 4994 end do 4995 4996 if (Wfd%usepaw==1) then 4997 ABI_WARNING("Testing the calculation of AE PAW wavefunctions.") 4998 ! Use a local pawfgrtab to make sure we use the correction in the paw spheres 4999 ! the usual pawfgrtab uses r_shape which may not be the same as r_paw. 5000 cplex=1 5001 call pawtab_get_lsize(Pawtab,l_size_atm,Cryst%natom,Cryst%typat) 5002 ABI_MALLOC(Pawfgrtab,(Cryst%natom)) 5003 call pawfgrtab_init(Pawfgrtab,cplex,l_size_atm,Wfd%nspden,Cryst%typat) 5004 ABI_FREE(l_size_atm) 5005 5006 optcut=1 ! use rpaw to construct Pawfgrtab. 5007 optgr0=0; optgr1=0; optgr2=0 ! dont need gY terms locally. 5008 optrad=1 ! do store r-R. 5009 5010 call nhatgrid(Cryst%atindx1,Cryst%gmet,Cryst%natom,Cryst%natom,Cryst%nattyp,ngfftf,Cryst%ntypat,& 5011 optcut,optgr0,optgr1,optgr2,optrad,Pawfgrtab,Pawtab,Cryst%rprimd,Cryst%typat,Cryst%ucvol,Cryst%xred) 5012 5013 !Pawfgrtab is ready to use 5014 5015 if (Wfd%pawprtvol>0) then 5016 call pawfgrtab_print(Pawfgrtab,natom=Cryst%natom,unit=std_out,& 5017 prtvol=Wfd%pawprtvol,mode_paral="COLL") 5018 end if 5019 5020 ABI_MALLOC(Paw_onsite,(Cryst%natom)) 5021 call paw_pwaves_lmn_init(Paw_onsite,Cryst%natom,Cryst%natom,Cryst%ntypat,& 5022 Cryst%rprimd,Cryst%xcart,Pawtab,Pawrad,Pawfgrtab) 5023 5024 ABI_MALLOC(ur_ae,(Wfd%nfft*Wfd%nspinor)) 5025 ABI_MALLOC(data_plot,(Wfd%nfft)) 5026 5027 do plot=1,my_nplots 5028 band =my_plot_list(1,plot) 5029 ik_ibz=my_plot_list(2,plot) 5030 spin =my_plot_list(3,plot) 5031 5032 call wfd%paw_get_aeur(band,ik_ibz,spin,Cryst,Paw_onsite,Psps,Pawtab,Pawfgrtab,ur_ae) 5033 5034 data_plot = DBLE(ur_ae(1:Wfd%nfft)*CONJG(ur_ae(1:Wfd%nfft)))/Cryst%ucvol 5035 if (Wfd%nspinor==2) data_plot = data_plot + DBLE(ur_ae(Wfd%nfft+1:)*CONJG(ur_ae(Wfd%nfft+1:)))/Cryst%ucvol 5036 5037 write(xsf_fname,'(3(a,i0),a)') 'PAW_AE_wfk2_sp',spin,'_kpt',ik_ibz,'_bd',band,'.xsf' 5038 if (open_file(xsf_fname,msg,newunit=funt,status='unknown',form='formatted') /= 0) then 5039 ABI_ERROR(msg) 5040 end if 5041 5042 call printxsf(n1,n2,n3,data_plot,Cryst%rprimd,(/zero,zero,zero/),& 5043 Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xcart,Cryst%znucl,funt,0) 5044 5045 close(funt) 5046 end do 5047 5048 ABI_FREE(ur_ae) 5049 ABI_FREE(data_plot) 5050 5051 call pawfgrtab_free(Pawfgrtab) 5052 ABI_FREE(Pawfgrtab) 5053 call paw_pwaves_lmn_free(Paw_onsite) 5054 ABI_FREE(Paw_onsite) 5055 5056 else 5057 ! NC case. Just a simple FFT G-->R and then dump the results. 5058 ABI_MALLOC(nc_ur,(Wfd%nfft*Wfd%nspinor)) 5059 ABI_MALLOC(data_plot,(Wfd%nfft)) 5060 5061 do plot=1,my_nplots 5062 band =my_plot_list(1,plot) 5063 ik_ibz=my_plot_list(2,plot) 5064 spin =my_plot_list(3,plot) 5065 5066 call wfd%get_ur(band,ik_ibz,spin,nc_ur) 5067 5068 data_plot = DBLE(nc_ur(1:Wfd%nfft)*CONJG(nc_ur(1:Wfd%nfft)))/Cryst%ucvol 5069 if (Wfd%nspinor==2) data_plot = data_plot + DBLE(nc_ur(Wfd%nfft+1:)*CONJG(nc_ur(Wfd%nfft+1:)))/Cryst%ucvol 5070 5071 write(xsf_fname,'(3(a,i0),a)') 'NC_wfk2_sp',spin,'_kpt',ik_ibz,'_bd',band,'.xsf' 5072 if (open_file(xsf_fname,msg,newunit=funt,status='unknown',form='formatted') /= 0) then 5073 ABI_ERROR(msg) 5074 end if 5075 call printxsf(n1,n2,n3,data_plot,Cryst%rprimd,(/zero,zero,zero/),& 5076 Cryst%natom,Cryst%ntypat,Cryst%typat,Cryst%xcart,Cryst%znucl,funt,0) 5077 5078 close(funt) 5079 end do 5080 5081 ABI_FREE(nc_ur) 5082 ABI_FREE(data_plot) 5083 end if 5084 5085 ABI_FREE(my_plot_list) 5086 5087 end subroutine wfdgw_plot_ur
m_wfd/wfdgw_rank_has_ug [ Functions ]
NAME
wfdgw_rank_has_ug
FUNCTION
This function is used to ask a particular processor whether it has a particular ug and with which status.
INPUTS
rank=The MPI rank of the processor. band=Band index. ik_ibz=k-point index spin=Spin index.
NOTES
A zero index can be used to inquire the status of a bunch of states. Thus (band,ik_ibz,spin) = (0,1,1) means: Do you have at least one band for the first k-point and the first spin.
SOURCE
2338 function wfdgw_rank_has_ug(Wfd,rank,band,ik_ibz,spin) 2339 2340 !Arguments ------------------------------------ 2341 !scalars 2342 integer,intent(in) :: band,ik_ibz,spin,rank 2343 logical :: wfdgw_rank_has_ug 2344 class(wfdgw_t),intent(in) :: Wfd 2345 2346 !Local variables ------------------------------ 2347 !scalars 2348 integer :: nzeros 2349 integer(c_int8_t) :: bks_flag 2350 !arrays 2351 integer :: indices(3) 2352 2353 !************************************************************************ 2354 2355 indices = [band,ik_ibz,spin] 2356 bks_flag = WFD_STORED 2357 2358 if (ALL(indices/= [0,0,0])) then 2359 wfdgw_rank_has_ug = (Wfd%bks_tab(band,ik_ibz,spin,rank) == bks_flag); RETURN 2360 else 2361 nzeros = COUNT(indices==0) 2362 if (nzeros==3) ABI_ERROR("All indices are zero!") 2363 2364 if (band==0) then 2365 if (nzeros==1) wfdgw_rank_has_ug = ANY( Wfd%bks_tab(:,ik_ibz,spin,rank)==bks_flag) 2366 if (ik_ibz==0) wfdgw_rank_has_ug = ANY( Wfd%bks_tab(:,:,spin,rank) ==bks_flag) 2367 if (spin ==0) wfdgw_rank_has_ug = ANY( Wfd%bks_tab(:,ik_ibz,:,rank) ==bks_flag) 2368 2369 else if (ik_ibz==0) then 2370 if (nzeros==1) wfdgw_rank_has_ug = ANY( Wfd%bks_tab(band,:,spin,rank)==bks_flag) 2371 if (band ==0) wfdgw_rank_has_ug = ANY( Wfd%bks_tab(:,:,spin,rank) ==bks_flag) 2372 if (spin ==0) wfdgw_rank_has_ug = ANY( Wfd%bks_tab(band,:,:,rank) ==bks_flag) 2373 2374 else 2375 if (nzeros==1) wfdgw_rank_has_ug = ANY( Wfd%bks_tab(band,ik_ibz,:,rank)==bks_flag) 2376 if (ik_ibz==0) wfdgw_rank_has_ug = ANY( Wfd%bks_tab(band,:,:,rank) ==bks_flag) 2377 if (band ==0) wfdgw_rank_has_ug = ANY( Wfd%bks_tab(:,ik_ibz,:,rank) ==bks_flag) 2378 end if 2379 end if 2380 2381 end function wfdgw_rank_has_ug
m_wfd/wfdgw_rotate [ Functions ]
NAME
wfdgw_rotate
FUNCTION
This routine performs a linear transformation of the wavefunctions stored in Wfd taking into account memory distribution. The transformation is done in G-space therefore all the ug should be available. Wavefunctions in real space are then obtained via FFT. The implementation assumes that the matrix associated to the linear transformation is sparse (No BLAS-3 calls here).
INPUTS
Cryst<crystal_t>=Object defining the unit cell and its symmetries. m_ks_to_qp(mband,mband,nkibz,nsppol)= expansion of the QP amplitudes in terms of KS wavefunctions. [bmask(mband,nkibz,nsppol)]=The routine will raise an error if one band index is not treated by any processor. bmask can be used to select the subset of indices that are expected to be available.
SIDE EFFECTS
Wfd<wfd_t>=See above.
SOURCE
2993 subroutine wfdgw_rotate(Wfd, Cryst, m_ks_to_qp, bmask) 2994 2995 !Arguments ------------------------------------ 2996 !scalars 2997 class(wfdgw_t),intent(inout) :: Wfd 2998 type(crystal_t),intent(in) :: Cryst 2999 !arrays 3000 complex(dpc),target,intent(in) :: m_ks_to_qp(Wfd%mband,Wfd%mband,Wfd%nkibz,Wfd%nsppol) 3001 logical,optional,intent(in) :: bmask(Wfd%mband,Wfd%nkibz,Wfd%nsppol) 3002 3003 !Local variables------------------------------- 3004 !scalars 3005 integer :: band,ik_ibz,spin,ierr,icol,nnew,inew,my_nband,ib,npw_k,istwf_k 3006 character(len=500) :: msg 3007 type(wave_t),pointer :: wave 3008 !arrays 3009 integer :: new_list(Wfd%mband),my_band_list(Wfd%mband) 3010 complex(dpc),ABI_CONTIGUOUS pointer :: umat_sk(:,:) 3011 complex(gwpc) :: mcol(Wfd%mband) 3012 complex(gwpc),allocatable :: new_ug(:,:) !, new_ur(:) 3013 3014 !************************************************************************ 3015 3016 ! Update the distribution table, first. 3017 call wfd%update_bkstab() 3018 3019 ! Calculate: $\Psi^{QP}_{r,b} = \sum_n \Psi^{KS}_{r,n} M_{n,b}$ 3020 do spin=1,Wfd%nsppol 3021 do ik_ibz=1,Wfd%nkibz 3022 npw_k = Wfd%npwarr(ik_ibz) 3023 istwf_k = Wfd%istwfk(ik_ibz) 3024 if (istwf_k /= 1) then 3025 ABI_WARNING("wfdgw_rotate with istwfk /= 1") 3026 end if 3027 umat_sk => m_ks_to_qp(:,:,ik_ibz,spin) 3028 3029 ! Select only those states that are mixed by the (sparse) m_ks_to_qp. 3030 nnew=0; new_list=0 3031 do icol=1,Wfd%nband(ik_ibz,spin) 3032 mcol = m_ks_to_qp(:,icol,ik_ibz,spin) 3033 mcol(icol) = mcol(icol) - cone 3034 if (ANY(ABS(mcol)>tol12)) then ! Avoid a simple copy. 3035 nnew=nnew+1 3036 new_list(nnew)=icol 3037 end if 3038 end do 3039 if (nnew==0) CYCLE ! Nothing to do. 3040 3041 ! Retrieve the set of band indices that have to be treated by 3042 ! this node taking into account a possible duplication. 3043 if (present(bmask)) then 3044 call wfd%distribute_bands(ik_ibz,spin,my_nband,my_band_list,bmask=bmask(:,ik_ibz,spin)) 3045 else 3046 call wfd%distribute_bands(ik_ibz,spin,my_nband,my_band_list) 3047 end if 3048 3049 !if (my_nband>0) then 3050 ! write(std_out,*)" At (ik_ibz,spin) ",ik_ibz,spin,& 3051 ! & ", rank ",Wfd%my_rank," will sum ",my_nband," bands, my_band_list: ",my_band_list(1:my_nband) 3052 !end if 3053 ABI_MALLOC(new_ug,(npw_k*Wfd%nspinor,nnew)) 3054 new_ug=czero 3055 do inew=1,nnew 3056 icol = new_list(inew) 3057 do ib=1,my_nband 3058 band = my_band_list(ib) 3059 if (ABS(umat_sk(band,icol))>tol12) then 3060 ABI_CHECK(wfd%get_wave_ptr(band, ik_ibz, spin, wave, msg) == 0, msg) 3061 new_ug(:,inew) = new_ug(:,inew) + umat_sk(band, icol) * wave%ug 3062 end if 3063 end do 3064 end do 3065 3066 !if (istwf_k /= 1) then 3067 ! ABI_MALLOC(new_ur, (wfd%nfft * wfd%nspinor * nnew)) 3068 ! call fft_ug_dpc(npw_k, wfd%nfft, wfd%nspinor, nnew, wfd%mgfft, wfd%ngfft, istwf_k, & 3069 ! wfd%kdata(ik_ibz)%kg_k, wfd%kdata(ik_ibz)%gbound, new_ug, new_ur) 3070 ! new_ur = real(new_ur) 3071 ! call fft_ur_dpc(npw_k, wfd%nfft, wfd%nspinor, nnew, wfd%mgfft, wfd%ngfft, istwf_k, & 3072 ! wfd%kdata(ik_ibz)%kg_k, wfd%kdata(ik_ibz)%gbound, new_ur, new_ug) 3073 ! ABI_FREE(new_ur) 3074 !end if 3075 3076 call xmpi_sum(new_ug,Wfd%comm,ierr) 3077 3078 ! Update the input wave functions 3079 do inew=1,nnew 3080 band = new_list(inew) 3081 if (wfd%ihave_ug(band, ik_ibz, spin)) call wfd%push_ug(band, ik_ibz, spin, Cryst, new_ug(:,inew)) 3082 end do 3083 3084 ABI_FREE(new_ug) 3085 end do !ik_ibz 3086 end do !spin 3087 3088 ! Reinit the storage mode of Wfd as ug have been changed. 3089 ! This is needed only if FFTs are not done in wfd_push_ug. Do not know which one is faster. 3090 !call wfd%reset_ur_cprj() 3091 call xmpi_barrier(Wfd%comm) 3092 3093 end subroutine wfdgw_rotate
m_wfd/wfdgw_sanity_check [ Functions ]
NAME
wfdgw_sanity_check
FUNCTION
Debugging tool
INPUTS
Wfd<wfd_t>=
OUTPUT
SOURCE
3243 subroutine wfdgw_sanity_check(Wfd) 3244 3245 !Arguments ------------------------------------ 3246 !scalars 3247 class(wfdgw_t),intent(inout) :: Wfd 3248 3249 !Local variables ------------------------------ 3250 !scalars 3251 integer :: ik_ibz,spin,band,mpi_ierr,ierr,how_manyb,unt_dbg,irank 3252 character(len=500) :: msg 3253 !arrays 3254 integer :: my_band_list(Wfd%mband) 3255 3256 !************************************************************************ 3257 3258 call wfd%update_bkstab() 3259 ierr=0 3260 3261 do spin=1,Wfd%nsppol 3262 do ik_ibz=1,Wfd%nkibz 3263 do band=1,Wfd%nband(ik_ibz,spin) 3264 if (Wfd%bks_tab(band, ik_ibz, spin, Wfd%my_rank) == WFD_STORED .and. & 3265 .not. wfd%ihave_ug(band, ik_ibz, spin, how="Stored") ) then 3266 write(msg,'(a,3(i0,1x))')" Found inconsistency in bks_tab for (band, ik_ibz, spin): ",band,ik_ibz,spin 3267 call wrtout(std_out, msg) 3268 ierr = ierr+1 3269 end if 3270 end do 3271 end do 3272 end do 3273 3274 call xmpi_sum(ierr,Wfd%comm,mpi_ierr) 3275 3276 if (ierr/=0) then 3277 if (open_file("__WFD_DEBUG__",msg,newunit=unt_dbg,form="formatted") /= 0) then 3278 ABI_ERROR(msg) 3279 end if 3280 3281 do irank=0,Wfd%nproc-1 3282 if (irank==Wfd%my_rank) then 3283 write(unt_dbg,*)" (k,b,s) states owned by rank: ",Wfd%my_rank 3284 3285 do spin=1,Wfd%nsppol 3286 do ik_ibz=1,Wfd%nkibz 3287 write(unt_dbg,*)" (spin,ik_ibz) ",spin,ik_ibz 3288 call wfd%mybands(ik_ibz, spin, how_manyb, my_band_list, how="Stored") 3289 write(unt_dbg,*) (my_band_list(band),band=1,how_manyb) 3290 end do 3291 end do 3292 3293 end if 3294 end do 3295 close(unt_dbg) 3296 call xmpi_barrier(Wfd%comm) 3297 ABI_ERROR("Sanity check failed. Check WFD_DEBUG") 3298 end if 3299 3300 end subroutine wfdgw_sanity_check
m_wfd/wfdgw_show_bkstab [ Functions ]
NAME
wfdgw_show_bkstab
FUNCTION
Print a table showing the distribution of the wavefunctions.
SOURCE
2506 subroutine wfdgw_show_bkstab(Wfd, unit) 2507 2508 !Arguments ------------------------------------ 2509 !scalars 2510 integer,intent(in) :: unit 2511 class(wfdgw_t),intent(in) :: Wfd 2512 2513 !Local variables ------------------------------ 2514 !scalars 2515 integer :: ik_ibz,spin,band,nband_k,width 2516 character(len=1) :: chlist(0:Wfd%nproc-1) 2517 character(len=500) :: fmt 2518 2519 !************************************************************************ 2520 2521 width = max(80, Wfd%nproc) 2522 2523 write(fmt,"(a,i0,a)")"(i5,3x,",Wfd%nproc,"(a))" 2524 2525 do spin=1,Wfd%nsppol 2526 do ik_ibz=1,Wfd%nkibz 2527 write(unit,"(a)")repeat("=",width) 2528 write(unit,"(2(a,i0))")"Spin: ",spin,", ik_ibz: ",ik_ibz 2529 write(unit,"(a)")"MPI rank ----> (A=allocated, S=Stored, N=NoWave)." 2530 nband_k = Wfd%nband(ik_ibz, spin) 2531 do band=1,nband_k 2532 where (Wfd%bks_tab(band, ik_ibz, spin,:) == WFD_NOWAVE) 2533 chlist = "N" 2534 elsewhere (Wfd%bks_tab(band, ik_ibz, spin,:) == WFD_ALLOCATED) 2535 chlist = "A" 2536 elsewhere (Wfd%bks_tab(band, ik_ibz, spin,:) == WFD_STORED) 2537 chlist = "S" 2538 end where 2539 write(unit,fmt)band,chlist(:) 2540 end do 2541 write(unit,"(a)")repeat("=",width) 2542 end do 2543 end do 2544 2545 end subroutine wfdgw_show_bkstab
m_wfd/wfdgw_t [ Types ]
NAME
wfdgw_t
FUNCTION
This is a subclass class with the bks_tab used to parallelize the GW/BSE code. Unfortunately, the size of the bks_tab increases with the number of procs This design facilated the implementation of the MPI-algorithms but it leads to a big scalability issue when nprocs/nband/nkibz is large. New algorithms implemented outside of the GW/BSE code should use wfd_t. wfdgw_t is kept to avoid breaking the GW/BSE code.
SOURCE
482 type, extends (wfd_t), public :: wfdgw_t 483 484 integer(c_int8_t), private, allocatable :: bks_tab(:,:,:,:) 485 ! bks_tab(mband,nkibz,nsppol,0:nproc-1) 486 ! Global table used to keep trace of the distribution of the (b,k,s) states on each node inside Wfd%comm. 487 ! 1 if the node has this state. 0 otherwise. 488 ! A node owns a wavefunction if the corresponding ug is allocated AND computed. 489 ! If a node owns ur but not ug, or ug is just allocated then its entry in the table is zero. 490 491 contains 492 493 procedure :: show_bkstab => wfdgw_show_bkstab 494 ! Print a table showing the distribution of the wavefunctions. 495 496 procedure :: distribute_bands => wfdgw_distribute_bands 497 ! Distribute a set of bands taking into account the distribution of the ug. 498 499 procedure :: iterator_bks => wfdgw_iterator_bks 500 ! Iterator used to loop over bands, k-points and spin indices 501 502 procedure :: bks_distrb => wfdgw_bks_distrb 503 ! Distribute bands, k-points and spins 504 505 procedure :: update_bkstab => wfdgw_update_bkstab 506 ! Update the internal table with info on the distribution of the ugs. 507 508 procedure :: rotate => wfdgw_rotate 509 ! Linear transformation of the wavefunctions stored in Wfd 510 511 procedure :: sanity_check => wfdgw_sanity_check 512 ! Debugging tool 513 514 procedure :: distribute_bbp => wfdgw_distribute_bbp 515 ! Distribute a set of (b,b') indices 516 517 procedure :: distribute_kb_kpbp => wfdgw_distribute_kb_kpbp 518 519 procedure :: plot_ur => wfdgw_plot_ur 520 ! Write u(r) to an external file in XSF format. 521 522 procedure :: mkrho => wfdgw_mkrho 523 ! Calculate the charge density on the fine FFT grid in real space. 524 525 procedure :: pawrhoij => wfdgw_pawrhoij 526 527 procedure :: get_nl_me => wfdgw_get_nl_me 528 529 procedure :: rank_has_ug => wfdgw_rank_has_ug 530 531 procedure :: bands_of_rank => wfdgw_bands_of_rank 532 533 procedure :: write_wfk => wfdgw_write_wfk 534 ! Write u(g) to a WFK file. 535 536 end type wfdgw_t 537 538 public :: wfdgw_copy
m_wfd/wfdgw_update_bkstab [ Functions ]
NAME
wfdgw_update_bkstab
FUNCTION
This routine should be called by all the nodes before any MPI operation involving the object. It updates the bks_tab storing information on the distribution of ug. INPUT [show]=If present and > 0, print tabs to unit show.
SIDE EFFECTS
Wfd%bks_tab
SOURCE
2814 subroutine wfdgw_update_bkstab(Wfd, show) 2815 2816 !Arguments ------------------------------------ 2817 !scalars 2818 class(wfdgw_t),intent(inout) :: Wfd 2819 integer,optional,intent(in) :: show 2820 2821 !Local variables ------------------------------ 2822 !scalars 2823 integer :: ierr, nelem, spin, ik_ibz, band, is, ik, ib 2824 integer(c_int8_t),allocatable :: my_vtab(:),gather_vtabs(:) 2825 !logical,allocatable :: tab_ranks(:) 2826 2827 !************************************************************************ 2828 2829 ! Fill my slice of the global table. 2830 do spin=1,wfd%nsppol 2831 do ik_ibz=1,wfd%nkibz 2832 do band=1,Wfd%nband(ik_ibz, spin) 2833 ib = wfd%bks2wfd(1, band, ik_ibz, spin) 2834 ik = wfd%bks2wfd(2, band, ik_ibz, spin) 2835 is = wfd%bks2wfd(3, band, ik_ibz, spin) 2836 if (ib /= 0) then 2837 Wfd%bks_tab(band, ik_ibz, spin, Wfd%my_rank) = wfd%s(is)%k(ik)%b(ib)%has_ug 2838 else 2839 Wfd%bks_tab(band, ik_ibz, spin, Wfd%my_rank) = WFD_NOWAVE 2840 end if 2841 end do 2842 end do 2843 end do 2844 2845 ! Gather flags on each node. 2846 nelem = Wfd%mband*Wfd%nkibz*Wfd%nsppol 2847 ABI_MALLOC(my_vtab, (nelem)) 2848 my_vtab(:) = reshape(Wfd%bks_tab(:,:,:,Wfd%my_rank), [nelem]) 2849 2850 ABI_MALLOC(gather_vtabs, (nelem*Wfd%nproc)) 2851 2852 call xmpi_allgather(my_vtab,nelem,gather_vtabs,Wfd%comm,ierr) 2853 2854 Wfd%bks_tab(:,:,:,:) = reshape(gather_vtabs, [Wfd%mband, Wfd%nkibz, Wfd%nsppol, Wfd%nproc]) 2855 ABI_FREE(my_vtab) 2856 ABI_FREE(gather_vtabs) 2857 2858 #if 0 2859 ! This is gonna be slow but if lot of k-points as I cannot assume bands or k-points have been filtered 2860 ! Need to introduce global_filter_ikibz_spin in wfd_init ... 2861 ABI_MALLOC(tab_ranks, (wfd%nproc)) 2862 do spin=1,Wfd%nsppol 2863 do ik_ibz=1,Wfd%nkibz 2864 !if wfd%global_filter_ikibz_spin(ik_ibz, spin) cycle 2865 do band=1,Wfd%nband(ik_ibz, spin) 2866 tab_ranks = .False. 2867 if (len(wfd%bks_ranks(band, ik_ibz, spin) > 0) then 2868 if (any(wfd%my_rank == wfd%bks_ranks(band, ik_ibz, spin)) tab_ranks(wfd%my_rank) = .True. 2869 end if 2870 call xmpi_lor(tab_ranks, wfd%comm) 2871 call bool2index(tab_ranks, wfd%bks_ranks(band, ik_ibz, spin)) 2872 end do 2873 end do 2874 end do 2875 ABI_FREE(tab_ranks) 2876 #endif 2877 2878 if (present(show)) then 2879 if (show >= 0) call wfd%show_bkstab(unit=show) 2880 end if 2881 2882 end subroutine wfdgw_update_bkstab
m_wfd/wfdgw_who_has_ug [ Functions ]
NAME
wfdgw_who_has_ug
FUNCTION
Return the number of processors having a particular (b,k,s) state as well as their MPI rank. Warning: Wfd%bks_tab is supposed to be up-to-date (see wfdgw_update_bkstab).
INPUTS
band=the index of the band. ik_ibz=Index of the k-point in the IBZ spin=spin index
OUTPUT
how_many=The number of nodes owing this ug state. proc_ranks(1:how_many)=Gives the MPI rank of the nodes owing the state.
SOURCE
2737 subroutine wfdgw_who_has_ug(Wfd,band,ik_ibz,spin,how_many,proc_ranks) 2738 2739 !Arguments ------------------------------------ 2740 !scalars 2741 integer,intent(in) :: band,ik_ibz,spin 2742 integer,intent(out) :: how_many 2743 class(wfdgw_t),intent(in) :: Wfd 2744 !arrays 2745 integer,intent(out) :: proc_ranks(Wfd%nproc) 2746 2747 !Local variables ------------------------------ 2748 !scalars 2749 integer :: irank 2750 logical :: bks_select,spin_select,kpt_select 2751 character(len=500) :: msg 2752 2753 !************************************************************************ 2754 2755 bks_select = (band/=0.and.ik_ibz/=0.and.spin/=0) 2756 spin_select = (band==0.and.ik_ibz==0.and.spin/=0) 2757 kpt_select = (band==0.and.ik_ibz/=0.and.spin/=0) 2758 2759 how_many=0; proc_ranks=-1 2760 2761 if (bks_select) then 2762 ! List the proc owining this (b,k,s) state. 2763 do irank=0,Wfd%nproc-1 2764 if (Wfd%bks_tab(band, ik_ibz, spin, irank) == WFD_STORED) then 2765 how_many = how_many +1 2766 proc_ranks(how_many)=irank 2767 end if 2768 end do 2769 2770 else if (spin_select) then 2771 ! List the proc owining at least one state with this spin. 2772 do irank=0,Wfd%nproc-1 2773 if ( ANY(Wfd%bks_tab(:,:,spin,irank)==WFD_STORED) ) then 2774 how_many = how_many +1 2775 proc_ranks(how_many)=irank 2776 end if 2777 end do 2778 2779 else if (kpt_select) then 2780 ! List the proc owining at least one state with this (k-point, spin). 2781 do irank=0,Wfd%nproc-1 2782 if ( ANY(Wfd%bks_tab(:,ik_ibz,spin,irank)==WFD_STORED) ) then 2783 how_many = how_many +1 2784 proc_ranks(how_many)=irank 2785 end if 2786 end do 2787 2788 else 2789 write(msg,'(a,3(i0,1x))')" Wrong value for (b,k,s): ",band,ik_ibz,spin 2790 ABI_ERROR(msg) 2791 end if 2792 2793 end subroutine wfdgw_who_has_ug
m_wfd/wfdgw_write_wfk [ Functions ]
NAME
wfdgw_write_wfk
FUNCTION
This routine writes the wavefunctions to the specified WFK file All the wavefunction are stored on each node, only the spin is distributed.
INPUTS
Wfd<wfd_t>=Initialized wavefunction descritptor. wfk_fname=Name of the WFK file.
OUTPUT
Only writing
SOURCE
4224 subroutine wfdgw_write_wfk(Wfd, Hdr, ebands, wfk_fname, wfknocheck) 4225 4226 !Arguments ------------------------------------ 4227 !scalars 4228 character(len=*),intent(in) :: wfk_fname 4229 class(wfdgw_t),intent(in) :: Wfd 4230 type(Hdr_type),intent(in) :: Hdr 4231 type(ebands_t),intent(in) :: ebands 4232 logical,intent(in),optional :: wfknocheck 4233 4234 !Local variables ------------------------------ 4235 !scalars 4236 integer,parameter :: formeig0=0,master=0 4237 integer :: nprocs,my_rank,iomode,cgsize,npw_k,ik_ibz,spin,nband_k,band,ii 4238 integer :: blk,nblocks,how_many,ierr,how_manyb 4239 real(dp) :: cpu,wall,gflops 4240 logical :: iam_master,nocheck ! MRM 4241 character(len=500) :: msg 4242 type(wfk_t) :: Wfkfile 4243 !arrays 4244 integer :: band_block(2),proc_ranks(Wfd%nproc),my_band_list(Wfd%mband) 4245 integer,allocatable :: blocks(:,:) ! 4246 real(dp),allocatable :: cg_k(:,:) 4247 4248 !************************************************************************ 4249 4250 nocheck=.false. 4251 if(present(wfknocheck)) nocheck=wfknocheck 4252 4253 nprocs = xmpi_comm_size(Wfd%comm); my_rank = xmpi_comm_rank(Wfd%comm) 4254 iam_master = (my_rank == master) 4255 4256 ! Select the IO library from the file extension. 4257 iomode = iomode_from_fname(wfk_fname) 4258 call wrtout(std_out, sjoin('Writing GS WFK file: ',wfk_fname,", with iomode ",iomode2str(iomode))) 4259 4260 if (nprocs > 1 .and. iomode /= IO_MODE_MPI) then 4261 ABI_ERROR("You need MPI-IO to write wavefunctions in parallel") 4262 end if 4263 ! 4264 ! Check consistency between Wfd and Header! 4265 ! The ideal approach would be to generate the header from the Wfd but a lot of info are missing 4266 ABI_CHECK(Wfd%nkibz == Hdr%nkpt,"Different number of k-points") 4267 ABI_CHECK(Wfd%nsppol == Hdr%nsppol,"Different number of spins") 4268 ABI_CHECK(Wfd%nspinor == Hdr%nspinor,"Different number of spinors") 4269 4270 !if (any(Wfd%nband /= reshape(Hdr%nband, [Wfd%nkibz, Wfd%nsppol]))) then 4271 ! ABI_ERROR("Wfd%nband /= Hdr%nband") 4272 !end if 4273 4274 !endif 4275 ! Use bks_tab to decide who will write the data. Remember 4276 ! integer,allocatable :: bks_tab(:,:,:,:) 4277 ! Wfd%bks_tab(mband,nkibz,nsppol,0:nproc-1) 4278 ! Global table used to keep trace of the distribution of the (b,k,s) states on each node inside Wfd%comm. 4279 ! 1 if the node has this state. 0 otherwise. 4280 ! A node owns a wavefunction if the corresponding ug is allocated AND computed. 4281 ! If a node owns ur but not ug, or ug is just allocated then its entry in the table is zero. 4282 ! The main difficulties here are: 4283 ! 4284 ! 1) FFT parallelism (not coded, indeed) 4285 ! 2) Wavefunctions that are replicated, i.e. the same (b,k,s) is treated by more than one node. 4286 4287 ierr = 0 4288 do spin=1,Wfd%nsppol 4289 do ik_ibz=1,Wfd%nkibz 4290 do band=1,Wfd%nband(ik_ibz,spin) 4291 call wfdgw_who_has_ug(Wfd,band,ik_ibz,spin,how_many,proc_ranks) 4292 if (how_many /= 1) then 4293 ierr = ierr + 1 4294 write(msg,'(a,3(i0,1x))')" Found replicated state (b,k,s) ",band,ik_ibz,spin 4295 ABI_WARNING(msg) 4296 end if 4297 end do 4298 end do 4299 end do 4300 4301 if (ierr /= 0) then 4302 ABI_ERROR("Cannot write WFK file when wavefunctions are replicated") 4303 end if 4304 4305 call cwtime(cpu,wall,gflops,"start") 4306 4307 ! Master node opens the file and writes the Abinit header. 4308 if (iam_master) then 4309 do ik_ibz=1,Wfd%nkibz 4310 if (size(Wfd%Kdata(ik_ibz)%kg_k,dim=2)<Hdr%npwarr(ik_ibz)) then 4311 ABI_ERROR("Impossible to continue when the npw in the Hdr is diff. to the npw in the Wfd") 4312 end if 4313 end do 4314 call wfkfile%open_write(Hdr,wfk_fname,formeig0,iomode,get_unit(),xmpi_comm_self,write_hdr=.TRUE.,write_frm=.TRUE.) 4315 end if 4316 4317 ! Other nodes wait here before opening the same file. 4318 call xmpi_barrier(Wfd%comm) 4319 if (.not.iam_master) then 4320 call wfkfile%open_write(Hdr,wfk_fname,formeig0,iomode,get_unit(),xmpi_comm_self,write_hdr=.FALSE.,write_frm=.FALSE.) 4321 end if 4322 4323 do spin=1,Wfd%nsppol 4324 do ik_ibz=1,Wfd%nkibz 4325 ! MRM: Well, we do not check because nocheck is used when Wfd is stored only on the master. So it works for this case! 4326 if(.not.nocheck) then 4327 if (.not. wfd%ihave_ug(band, ik_ibz, spin, how="Stored")) cycle 4328 endif 4329 4330 nband_k = Wfd%nband(ik_ibz,spin) 4331 npw_k = Wfd%npwarr(ik_ibz) 4332 4333 ! Compute my block of bands for this k-point and spin. 4334 call wfd%mybands(ik_ibz, spin, how_manyb, my_band_list, how="Stored") 4335 call list2blocks(my_band_list(1:how_manyb), nblocks, blocks) 4336 4337 !if (proc_distrb_cycle(mpi_enreg%proc_distrb,ik_ibz,1,nband_k,spin,my_rank)) CYCLE 4338 !call mask2blocks(mpi_enreg%proc_distrb(ik_ibz,:,spin)==my_rank, nblocks,blocks) 4339 4340 ABI_CHECK(nblocks==1,"nblocks !=1") 4341 write(msg,"(a,3(i0,2x))")" Will write (ik_ibz, spin, nblocks) ",ik_ibz,spin,nblocks 4342 call wrtout(std_out, msg) 4343 4344 ! Extract the block of wavefunctions from Wfd. 4345 ! Try to allocate all u(g) first, 4346 ! TODO If not enough memory fallback to a blocked algorithm. 4347 cgsize = Wfd%nspinor * npw_k * how_manyb 4348 ABI_MALLOC_OR_DIE(cg_k, (2,cgsize), ierr) 4349 4350 if (size(Wfd%Kdata(ik_ibz)%kg_k,dim=2)<wfkfile%Hdr%npwarr(ik_ibz)) then 4351 ABI_ERROR("Wrong number of npw before printing") 4352 end if 4353 ! Extract the set of u(g) for this (kpoint,spin) 4354 ! This works only if all the bands are on the same node. 4355 !band_block = [1, nband_k] 4356 !call wfd_extract_cgblock(Wfd,[(ii, ii=1,nband_k)],ik_ibz,spin,cg_k) 4357 do blk=1,nblocks 4358 band_block = blocks(:,blk) 4359 call wfd_extract_cgblock(Wfd,[(ii, ii=band_block(1),band_block(2))],ik_ibz,spin,cg_k) ! cg_k extracted from Wfd! 4360 4361 if (band_block(1)==1) then 4362 ! Write also kg_k, eig_k and occ_k 4363 call wfkfile%write_band_block(band_block,ik_ibz,spin,xmpio_single,& 4364 kg_k=Wfd%Kdata(ik_ibz)%kg_k,cg_k=cg_k, & 4365 eig_k=ebands%eig(:,ik_ibz,spin),occ_k=ebands%occ(:,ik_ibz,spin)) ! occs extracted from Bands (i.e. QP_BSt) 4366 ! kg_k obtained from Wfd so OK! It is 4367 ! how Gs are ordered. 4368 else 4369 ABI_ERROR("band_block(1)>1 should not happen in the present version!") 4370 !call wfkfile%write_band_block(band_block,ik_ibz,spin,xmpio_single,cg_k=cg_k(:,1+icg:)) 4371 end if 4372 end do 4373 4374 ABI_FREE(cg_k) 4375 ABI_FREE(blocks) 4376 end do ! k-points 4377 end do ! spin 4378 4379 call xmpi_barrier(Wfd%comm) 4380 4381 ! Close the file. 4382 call wfkfile%close() 4383 4384 call cwtime_report(" write all cg" , cpu, wall, gflops) 4385 4386 end subroutine wfdgw_write_wfk